Vincent Monchiet
email: vincent.monchiet@univ-mlv.fr

Djimedo Kondo
email: djimedo.kondo@upmc.fr

Exact limit-analysis solution of a plastic hollow sphere with a Mises-Schleicher matrix

Keywords: Gurson Model, Limit Analysis, Hollow Sphere, Pressure-Sensitive Plasticity, Mises-Schleicher, Parabolic Type Criterion

In this paper we establish the exact solution for a hollow sphere with a rigid-plastic pressure-sensitive matrix and subjected to hydrostatic traction or compression. The matrix is assumed to obey to a parabolic Mises-Schleicher criterion. The closedform expressions of the velocity field and of the stress field are provided. These exact solutions, expressed by means of the Lambert W function, allow to assess and discuss existing results.

Introduction

The study by Gurson [7] is known to provide an upper bound for the macroscopic yield function of porous plastic materials having von Mises matrix. The approach is based on limit analysis of a hollow sphere made up of a rigid ideal plastic von Mises material and subjected to an arbitrary loading in the context of homogeneous strain rate boundary conditions. The macroscopic yield function has been derived by considering two velocity fields: the first one corresponds to a homogeneous strain rate field, while the second one, which is heterogeneous, accounts for the cavity expansion. For the latter, use has been made of the exact solution of the hollow sphere subjected to a hydrostatic external loading. As a consequence, the obtained macroscopic yield function provides the exact solution for this particular loading case. Improvements of the Gurson criterion have been recently proposed by Monchiet al. [START_REF] Monchiet | An improvement of Gurson-type models of porous materials by using Eshelbylike trial velocity fields[END_REF], [START_REF] Monchiet | A Micromechanics-based modification of the Gurson criterion by using Eshelbylike velocity fields[END_REF] by making use of Eshelby-like trial velocity field. Despite its great interest, the Gurson criterion appears to be inappropriate for a large class of ductile porous materials (for instance some polymers or some geomaterials) having a pressure-dependent matrix. To circumvent this limitation, several studies have been further proposed in order to extend the Gurson approach to porous materials with a Drucker-Prager type plastic matrix [START_REF] Guo | Continuum modeling of a porous solid with pressure sensitive dilatant matrix[END_REF] (see also [START_REF] Jeong | A macroscopic constitutive law for porous solids with pressure-sensitive matrices and its implications to plastic flow localization[END_REF][START_REF] Jeong | A new yield function and a hydrostatic stress-controlled model for porous solids with pressure-sensitive matrices[END_REF]). Still, these authors used a two-field trial velocity which contains the exact solution for the Drucker-Prager hollow sphere subjected to a purely hydrostatic loading (see [START_REF] Thoré | Closed-form solutions for the hollow sphere model with Coulomb and Drucker-Prager materials under isotropic loadings[END_REF][START_REF] Pastor | Limit analysis and numerical modeling of spherically porous solids with Coulomb and Drucker-Prager matrices[END_REF]). It must be emphasized that the macroscopic yield function derived by [START_REF] Guo | Continuum modeling of a porous solid with pressure sensitive dilatant matrix[END_REF] allows to recover the Gurson criterion in the limit case of von Mises matrix. Unfortunately, the above extensions of the Gurson model do not apply to porous media with a pressure-dependent matrix obeying to a parabolic type criterion (instead of the Drucker-Prager linear one), the so-called Mises-Schleicher criterion [START_REF] Lubliner | Plasticity Theory[END_REF] 1 (see also [START_REF] Raghava | The macroscopic yield behavior of polymers[END_REF][START_REF] Raghava | A Macroscopic Yield Criterion for Crystalline Polymers[END_REF]) which exhibits also an asymmetry between tension and compression2 . It must be mentioned that this type of matrix has been considered by Lee and Oung [START_REF] Lee | Yield functions and flow rules for porous pressure-dependent strain-hardening polymeric materials[END_REF] (see also [START_REF] Durban | Plastic response of porous solids with pressure sensitive matrix[END_REF]) to model the ductile failure behavior of glassy polymers. To this end, these authors used a trial velocity field inspired from the one already considered by Gurson and accounting for the matrix pressure sensitivity. A crucial observation, emphasized by the authors themselves, is that their macroscopic yield function does not recover the Gurson criterion in the limit of the von-Mises matrix. In particular, the exact solution established by [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF] for a hydrostatic loading is not retrieved. Based on these observations, the authors proposed a heuristical modification of the original criterion which allows to retrieve the Gurson exact solution. In the case of the pressure-dependent matrix, this modification still misses the exact solution for a hydrostatic external loading. In the present study, we derive the exact stress and velocity fields for the hollow sphere (with a matrix which obeys to a Mises-Schleicher criterion) subjected to an external hydrostatic traction or compression. The closed-form expressions of the velocity field is shown to be expressed in term of the Lambert W function (see for instance [START_REF] Corless | On the Lambert W function[END_REF]). These new results allows to assess the relevance of the Lee and Oung (2000)'s criterion. Their use to derive an upper bound of the macroscopic yield function for arbitrary loadings is addressed.

Statement of the problem

The local constitutive behavior

We consider a homogeneous plastic material obeying to the following pressure sensitive criterion:

f (σ) = σ 2 eq + 3ασ 0 σ m -σ 2 0 = 0 (1)
known as the Mises-Schleicher criterion [17] (see also [START_REF] Lubliner | Plasticity Theory[END_REF] (page 28) [START_REF] Raghava | The macroscopic yield behavior of polymers[END_REF][START_REF] Raghava | A Macroscopic Yield Criterion for Crystalline Polymers[END_REF]) in which α and σ 0 are related to the tensile yield stress, T , and absolute yield stress in compression, C, by:

σ 0 = √ CT , α = C -T √ CT (2)
Note that the Mises-Schleicher criterion predicts asymmetry between tension and compression.

In [START_REF] Barthélémy | Détermination du critère de rupture macroscopique d'un milieu poreux par homognisation non linaire[END_REF], σ eq and σ m denote the equivalent von-Mises stress and the mean stress respectively. σ 0 and α are two material parameters with the assumption that α ≥ 0 which physically means that the yield stress in traction is lower than in compression. The strain rate is derived from the normality flow rule:

d = λ ∂f ∂σ (σ) = λ( 3σ + ασ 0 I ) ( 3
)
where σ is the deviatoric part of the local stress, σ = σ -σ m I and I is the second order identity tensor. In (3), λ is the plastic multiplier. The computation of the von-Mises equivalent strain rate, d eq , and the mean part of the local strain rate, d m , leads to:

d eq = 2 λσ eq , d m = λασ 0 (4)
The plastic multiplier being positive, the means strain rate is then positive, d m ≥ 0. The local dissipation reads:

π(d) = σ : d = σ eq d eq + 3σ m d m = (2σ 2 eq + 3ασ 0 σ m ) λ = λ(σ 2 0 + σ 2 eq ) (5)
Eliminating λ between the two equations in (4), one has:

σ eq = ασ 0 2 d eq d m (6)
Replacing now the equivalent stress in (5) by its expression [START_REF] Guo | Continuum modeling of a porous solid with pressure sensitive dilatant matrix[END_REF] and the plastic multiplier by λ = d m /(ασ 0 ) (see equation ( 4)) one obtains:

π(d) = σ 0 d m α + ασ 0 4 d 2 eq d m (7)

Application of the Gurson limit analysis approach

Let us now consider a hollow sphere, with the internal and external radii a and b respectively, subjected at its outer boundary to homogeneous strain rate conditions:

v(r = b) = D.x (8)
For this homogeneous strain rate boundary conditions, the following inequality holds for all macroscopic stress Σ and macroscopic strain rate D [START_REF] Suquet | Plasticité et homogénéisation[END_REF], [START_REF] Buhan | A fundamental approach to the yield design of reinforced soil structures -Chap. 2: yield design homogenization theory for periodic media[END_REF]:

Σ : D ≤ Π(D) = inf v K.A [ 1 |Ω| ∫ Ω-ω π(d)dV ] (9)
Π(D) represents the macroscopic dissipation, Ω denotes the volume of the unit cell, |Ω| = 4πb 3 /3, whereas ω denotes the volume of the void, |ω| = 4πa 3 /3. In ( 9), the local dissipation π(d) is given by [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF]. The infimum is taken over all kinematically admissible (K.A) velocity fields, v, complying then with (8).

As classically, the limit stress states at the macroscopic scale are shown to be of the form:

Σ = ∂Π ∂D ( 10
)

Lee and Oung's macroscopic yield surface [10]

An upper bound for the macroscopic yield function has been derived by Lee and Oung [START_REF] Lee | Yield functions and flow rules for porous pressure-dependent strain-hardening polymeric materials[END_REF]. These authors considered a trial velocity field inspired of the one used by Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF] in the case of a von-Mises matrix:

v = (D + a 0 I).x + (D m -a 0 ) b 3 r 2 e r ( 11
)
This field is compatible with boundary conditions ( 8) and coincide with Gurson's trial velocity field when a 0 = 0. The constant a 0 is introduced in order to takes into account the compressibility of the matrix. Due to the presence of this constant, it must be emphasized that the velocity field given by ( 11) is not entirely determined by the boundary conditions. Its determination has required the minimization of the macroscopic dissipation (9) according to a 0 which is then given by (see [START_REF] Lee | Yield functions and flow rules for porous pressure-dependent strain-hardening polymeric materials[END_REF]):

a 0 = α 2 ( 4D 2 m + f 2 D 2 eq f + α 2 ) 1/2 (12)
The macroscopic yield function, derived from (10) reads:

Σ 2 eq + 9f 4 Σ 2 m + 3(1 -f )ασ 0 Σ m -(1 -f ) 2 σ 2 0 = 0 ( 13
)
Note that the above criterion constitutes an upper bound of the exact yield surface. However, the most important criticism which can be formulated about ( 13) is that it does not allow to recover the Gurson criterion [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF] when the limit α → 0 is taken 3 . Moreover, under purely hydrostatic loadings (traction or compression) and for α = 0, it provides:

|Σ m | = 2 3 σ 0 1 -f f ( 15
)
Which differs from the well-known exact solution:

|Σ m | = - 2 3 σ 0 ln(f ) (16)
Since the trial velocity field [START_REF] Lubliner | Plasticity Theory[END_REF] contains the exact solution for purely hydrostatic loading and α = 0 (then a 0 = 0), this observation is questionable. The explanation may be found in the fact that, although the velocity field [START_REF] Lubliner | Plasticity Theory[END_REF] coincides with the one used by Gurson for α = 0, it may not be valid for the 3 To correct this shortcoming of the upper bound (13), Lee and Oung [START_REF] Lee | Yield functions and flow rules for porous pressure-dependent strain-hardening polymeric materials[END_REF] has proposed a heuristical modification of their result in the form:

Σ 2 eq + 2f σ 2 0 cosh ( 3Σ m 2σ 0 ) + 3(1 -f )ασ 0 Σ m -(1 -f ) 2 σ 2 0 = 0 (14)
dissipation [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF]. Indeed, when the hollow sphere is subjected to an hydrostatic traction or compression (D = 0 and D = D m I), the exact expression of the local dissipation is:

π(d) = 2σ 0 |D m | b 3 r 3 (17)
whereas with the trial velocity field [START_REF] Lubliner | Plasticity Theory[END_REF], π(d) takes the following form when α → 0:

π(d) = σ 0 |D m | √ f ( 1 + f b 6 r 6 ) ( 18
)
3 Derivation of the local stress fields and the yield strength

The local stress field

Let us first introduce the spherical basis (e r , e θ , e φ ) and the associated coordinates system (r, θ, φ), with θ ∈ [0, 2π] and φ ∈ [0, π]. The macroscopic loading is assumed to be isotropic (the hollow sphere is subjected at its external boundary to hydrostatic traction or compression). Since the problem is invariant by any rotation with an angle φ or θ, the stress tensor has the form:

σ ij =        σ rr 0 0 0 σ θθ 0 0 0 σ φφ        (e r ,e θ ,e φ ) (19)
where σ rr and σ θθ = σ φφ are both functions of the coordinate r. The equivalent stress as well as the mean stress, computed from [START_REF] Thoré | Closed-form solutions for the hollow sphere model with Coulomb and Drucker-Prager materials under isotropic loadings[END_REF], read then:

σ eq = ϵ(σ θθ -σ rr ), σ m = 1 3 (σ rr + 2σ θθ ) ( 20
)
in which ϵ = sign(σ θθ -σ rr ). Let us introduce the positive function of the radial coordinate, G(r), such that

σ eq = σ 0 G(r) (21)
Due to the condition f (σ) = 0, the mean stress reads, for α ̸ = 0:

σ m = σ 0 3α (1 -G 2 (r)) (22)
The components of the stress tensor are:

σ rr = σ m - 2ϵ 3 σ eq = σ 0 3α [ 1 -G 2 (r) ] - 2ϵ 3 σ 0 G(r) σ θθ = σ m + ϵ 3 σ eq = σ 0 3α [ 1 -G 2 (r) ] + ϵ 3 σ 0 G(r) (23)
As classically, the local equilibrium is expressed as:

dσ rr dr + 2 r (σ rr -σ θθ ) = 0 (24)
Introducing ( 23) into (24) leads to the following non linear differential equation for G(r):

G(r)G ′ (r) + ϵαG ′ (r) + 3ϵα r G(r) = 0 (25)
The integration of the above equation allows to obtain the solution in the following implicit form:

ln(G(r)) + G(r) ϵα + ln(r 3 ) = K ( 26
)
Where K is a constant. The above solution can be put into the equivalent form:

G(r) ϵα exp ( G(r) ϵα ) = p a 3 r 3 ( 27
)
where the notation p = exp(K)/(ϵαa 3 ) has been adopted. The explicit solution reads then:

G(r) = ϵαW ( p a 3 r 3 ) ( 28
)
where W denotes the "Lambert W" function which is defined to be the inverse of the function x → x exp(x). Function W is not injective and is double-valued on the interval [-1/e, 0]. Indeed, if x > 0 the equation x exp(x) = y has only one solution for x but has two solutions when -1/e ≤ x > 0 (see figure 1). For this reason, two branches of function W are distinguished. The upper branch, corresponding to W (x) ≥ -1, is denoted W -1 (x) whereas W 0 (x) denotes the lower branch, corresponding to W (x) ≤ -1. On figure 1 are represented the two branches of the Lambert W function. Function G must be positive since G(r) = σ eq /σ 0 . However W (x) has the sign of x; it follows that, in (28), ϵ = sign(p). Finally, reporting (28) into (23) gives:

σ rr = σ 0 3α [ 1 -2α 2 W ( p a 3 r 3 ) -α 2 W 2 ( p a 3 r 3 )] σ θθ = σ 0 3α [ 1 + α 2 W ( p a 3 r 3 ) -α 2 W 2 ( p a 3 r 3 )] (29)
Coefficient p is determined with condition σ rr (r = a) = 0 (traction free condition on the void boundary). This identification leads to :

                 p = p + = z + exp(z + ), z + = -α + √ α 2 + 1 α or p = p -= z -exp(z -), z -= -α - √ α 2 + 1 α (30)
The variations of p + and p -according to α are represented on figure 2. It can be observed that coefficient p + is positive and tends to zero when α → +∞ whereas p -is negative and tends to the finite value z -= -2 exp(-2) when α → +∞.

The yield strength

Under uniform strain rate boundary conditions, the macroscopic stress is defined as the average of the local stress field:

Σ =< σ > Ω (31)
Here, the macroscopic strain rate is assumed to be isotropic, D = D m I. Therefore, the macroscopic stress is isotropic and takes the form Σ = Σ m I with:

Σ m =< σ m > Ω = b 3|Ω| ∫ ∂Ω σ rr (r = b)dS (32)
for which use has been done of the divergence theorem. Since σ rr is only function of the radial coordinate, r, one has:

Σ m = σ rr (r = b) (33)
(28) together with (30) lead to four expressions for each component of the stress field, depending on whether one considers the upper or lower branch of the Lambert W function but also whether one considers expressions p + or p -for coefficient p. Note that W -1 (x) takes complex values for x < 0 and consequently the solution W -1 (p -a 3 /r 3 ) must be discarded. The solution which uses the branch W 0 with coefficient p -can be also discarded since it leads to the solution Σ m = 0 when α = 0. The solution for traction is given by (29) in which the upper branch W 0 together with p + are used:

Σ m = σ 0 3α [ 1 -α 2 W 2 0 (f p + ) -2α 2 W 0 (f p + ) ] (34)
where f = a 3 /b 3 is the porosity. Remark 1: The Taylor expansion of W 0 (f p + ) at α = 0 gives:

W 0 (f p + ) = 1 α + ln(f ) -1 + o(α)
. It can then be shown that the expression of Σ m , given by (34), takes the finite limit Σ m = -2σ 0 3 ln(f ) when α → 0 which coincides with the result of Gurson (equation ( 16)). The solution for compression use the lower branch W -1 together with p -:

Σ m = σ 0 3α [ 1 -α 2 W 2 -1 (p -f ) -2α 2 W -1 (p -f ) ] ( 35
)
Remark 2: The Taylor expansion of

W -1 (f p -) at α = 0 is: W -1 (f p -) = -1 α + ln(f ) -1 + o(α).
It can then be shown that the expression of Σ m , given by (35), has a finite limit when α → 0. This limit is Σ m = 2σ 0 3 ln(f ) which coincides also with the result of Gurson (equation ( 16)).

Note that in [START_REF] Lee | Yield functions and flow rules for porous pressure-dependent strain-hardening polymeric materials[END_REF] the authors derived the macroscopic flow stress for purely hydrostatic loading by a direct method, i.e not requiring the determination of the local stress field. This solution can be recovered from ( 34) and ( 35) by reexpressing W (f p) as function of Σ m and taking thereafter the exponential. Doing that, the two solutions (34) and ( 35) can be reexpressed in the following implicit form:

• if Σ m ≥ 0 ln   -α + √ 1 + α 2 -3αΣ m /σ 0 -α + √ 1 + α 2   + 1 α √ 1 + α 2 -3αΣ m /σ 0 - 1 α √ 1 + α 2 = ln(f ) (36) • if Σ m < 0 ln   α + √ 1 + α 2 -3αΣ m /σ 0 α + √ 1 + α 2   - 1 α √ 1 + α 2 -3αΣ m /σ 0 + 1 α √ 1 + α 2 = ln(f ) (37)
On figure 3 is represented the yield stress Σ m /σ 0 as function of the parameter α for a porosity f = 0.1. The upper bound [START_REF] Monchiet | A Micromechanics-based modification of the Gurson criterion by using Eshelbylike velocity fields[END_REF] and the modified criterion [START_REF] Pastor | Limit analysis and numerical modeling of spherically porous solids with Coulomb and Drucker-Prager matrices[END_REF] both proposed by Lee and Oung [START_REF] Lee | Yield functions and flow rules for porous pressure-dependent strain-hardening polymeric materials[END_REF] are compared to the exact solution (34) for a macroscopic hydrostatic tension. It is observed that the upper bound (13) differs from the exact solution (34) when low values of the coefficient α are considered. However for larger values of α the upper bound (13) coincide with (34). The case of a hydrostatic compression is illustrated on figure 4 in which the upper bound ( 13) and the modified criterion ( 14) are compared to the exact solution (35). It appears that the upper bound differs from the exact solution, whatever the value of α. The modified criterion [START_REF] Pastor | Limit analysis and numerical modeling of spherically porous solids with Coulomb and Drucker-Prager matrices[END_REF] allows to retrieve the exact solution provided by Gurson criterion (α = 0). However, for strictly positive values of α, the approximate solution gives an estimation of the yield stress which is different from the exact solution.

Derivation of the velocity field

We now propose to derive the expression of the velocity field, solution of the rigid-plastic hollow sphere obeying to the Mises-Schleicher criterion. When the sphere is subjected to a hydrostatic loading, the velocity field is in the form:

v = F (r)e r ( 38
)
and the local strain rate field reads:

d ij =        d rr 0 0 0 d θθ 0 0 0 d φφ        (e r ,e θ ,e φ ) (39)
with:

d rr = F ′ (r), d θθ = d φφ = 1 r F (r) (40)
From relation (3), one has ϵ = sign(σ θθ -σ rr ) = sign(d θθ -d rr ). Consequently the equivalent and mean strain rate read:

d eq = 2 3 |d θθ -d rr | = 2ϵ 3 ( 1 r F (r) -F ′ (r) ) d m = 1 3 (d rr + 2d θθ ) = 1 3 ( F ′ (r) + 2 r F (r) ) (41)
Recalling that σ eq = σ 0 G(r) (see ( 21)) where the expression of G(r) is given in (28), and taking into account relation ( 6) with (41), it comes:

W ( p a 3 r 3 ) = F (r) -rF ′ (r) 2F (r) + rF ′ (r) (42)
Introducing in the above equation F (r) = exp(X(u)), with the change of variable u = a 3 /r 3 , one obtains:

X ′ (u) = W (pu) u(1 + W (pu)) - 1 3u (43)
Due to the property:

dW (x) dx = W (x) x(1 + W (x)) (44)
the integration of X(u) from (43) gives:

X(u) = 1 3 W (pu) - 1 3 ln(u) + K ′ (45)
It follows that:

F (r) = C r 2 W (pa 3 /r 3 ) (46)
Where C = a 2 exp(K ′ ) is a constant which has to be determined from boundary conditions.

C = b 3 D m W (pf ) (47)
for which it is recalled that the porosity reads f = a 3 /b 3 . The velocity field is then given by:

v r = D m b 3 r 2 W (pf ) W (pa 3 /r 3 ) (48)
In the case of a hydrostatic tension, the above solution is used with the upper branch W 0 and coefficient p + :

v r = D m b 3 r 2 W 0 (p + f ) W 0 (p + a 3 /r 3 ) (49)
whereas in the case of compression:

v r = D m b 3 r 2 W -1 (p -f ) W -1 (p -a 3 /r 3 ) (50)
d eq = 2|D m | b 3 r 3 + o(α); d m = |D m | b 3 r 3 α + o(α 3 ) (51)
whereas, from the trial velocity field (11), we get:

d eq = 2|D m | b 3 r 3 + o(α); d m = |D m | √ f α + o(α 3 ) (52)
It is observed that the terms of order 0 in the series of d eq and d m are the same for the exact and approximate solutions. However, the terms of order 1 in the series of d m differs. We aim now at comparing the exact and approximate expressions of the local dissipation. The sphere being still submitted to an external hydrostatic strain rate, the exact expression of the local dissipation is then obtained by replacing in [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF] the expressions of d eq and d m obtained with (38) with (48). The resulting expression reads:

π(d) = σ 0 |D m | b 3 r 3 W (pf )[1 + α 2 W (pa 3 /r 3 )] αW (pr 3 /a 3 )[1 + W (pr 3 /a 3 )] (53)
whereas the approximate expression of the local dissipation, obtained by Lee and Oung [START_REF] Lee | Yield functions and flow rules for porous pressure-dependent strain-hardening polymeric materials[END_REF], is:

π(d) = σ 0 |D m | r 6 /b 6 + f + 2α 2 -2α √ f + α 2 r 6 /b 6 √ f + α 2 (54)
On figure 7 are compared the two expressions (53) and (54) of the local dissipation for α = 0 and a porosity f = 0.01. It is observed that, although the approximate velocity field coincides with the exact solution for α = 0 (see figure 5), the local dissipation significantly differs. This is due to the fact that in the two solutions (see ( 51) and ( 52)) the ratio d m /α which directly enters in the expression of the local dissipation does not deliver the same limit when α → 0.

Conclusion

In this paper we derive the exact solution for a pressure-sensitive plastic hollow sphere subjected to a purely hydrostatic traction and compression. The hollow sphere is made up of a rigid-plastic material which obeys to the Mises-Schleicher criterion [17]. The closed-form expressions for the velocity and strain field, as well as for the stress field are shown to be expressed in term of the Lambert W function. These exact solutions are then used to assess the accuracy of approximate expressions in tension and compression established by [START_REF] Lee | Yield functions and flow rules for porous pressure-dependent strain-hardening polymeric materials[END_REF]. The exact velocity field found in the present study can be used to construct a trial velocity field for the limit analysis of the plastic hollow sphere under arbitrary loadings.

Fig. 1 .

 1 Fig. 1. At the left: function y = x exp(x). At the right: the two branches of the W Lambert function, y = W 0 (x) and y = W -1 (x).

Fig. 2 .

 2 Fig. 2. Representation of the coefficients p + and p -as function α

Fig. 3 .

 3 Fig. 3. Normalized yield stress Σ m /σ 0 as function of α for a porosity f = 0.1. Case of purely hydrostatic traction.

Fig. 4 .

 4 Fig. 4. Normalized yield stress Σ m /σ 0 as function of α for a porosity f = 0.1. Case of purely hydrostatic compression.

  The hollow sphere being subjected to the uniform strain boundary conditions v(r = b) = D.x with D = D m I, F (r) must satisfy F (r = b) = bD m . It comes, for the coefficient C:

Figures 5 and 6

 6 Figures5 and 6show the variations of v r /(bD m ) with r/a for traction and compression, respectively. The porosity value is f = 0.01 and different values of α are considered. On figures 5 and 6 are also reported the variations of the trial velocity field (in fact the unique radial component) used by Lee and Oung given by[START_REF] Lubliner | Plasticity Theory[END_REF] together with[START_REF] Monchiet | An improvement of Gurson-type models of porous materials by using Eshelbylike trial velocity fields[END_REF] in which has been put D = 0. As already mentioned, this velocity field coincides with the exact solution for α = 0. For completeness we propose to compare the values of d eq and d m , which enters into the definition of the local dissipation, for low values of the coefficient α. The Taylor expansions of d eq and d m at α = 0 reads, for the exact solution:

Fig. 5 .

 5 Fig. 5. Normalized radial component of the velocity field, v r /(bD m ) as function of r/a for f = 0.01 and various values of coefficient α in the case of traction.

Fig. 6 .

 6 Fig. 6. Normalized radial component of the velocity field, v r /(bD m ) as function of r/a for f = 0.01 and various values of coefficient α in the case of compression.

Fig. 7 .

 7 Fig. 7. Normalized dissipation, π(d)/(σ 0 D m ) as function of r/a for f = 0.01 and α = 0 in the case of a hydrostatic traction.

Figure 8

 8 Figure 8 corresponds to the comparison for α = 0.2.

Fig. 8 .

 8 Fig. 8. Normalized dissipation, π(d)/(σ 0 D m ) as function of r/a for f = 0.01 and α = 0.2 in the case of a hydrostatic traction.

The original publication of this criterion can be found in[17]

Note that recent investigations on ductile porous metallic materials with matrix having asymmetry between tension and compression has been carried out by[START_REF] Cazacu | Plastic potentials for porous aggregates with the matrix exhibiting tensioncompression asymmetry[END_REF].


Vincent Monchiet
email: vincent.monchiet@univ-mlv.fr

Eric Charkaluk
email: eric.charkaluk@univ-lille1.fr

Djimedo Kondo
email: djimedo.kondo@upmc.fr

A Micromechanics-based modification of the Gurson criterion by using Eshelby-like velocity fields

Keywords: ductile porous materials, Gurson type model, spherical cavities, Eshelby-like trial velocity fields, micromechanics, limit analysis

In this study, we propose a micromechanics-based modification of the Gurson criterion for porous media subjected to arbitrary loadings. The proposed formulation, derived in the framework of limit analysis, consists in the consideration of Eshelbylike exterior point trial velocity fields for the determination of the macroscopic dissipation. This approach is implemented for perfectly plastic rigid von Mises matrix containing spherical voids. After the minimization procedure required by the use of the Eshelby-like trial velocity fields, we derive a two-field estimate of the macroscopic yield function. It is shown that the obtained closed-form estimate provides a significant modification of the Gurson criterion, particularly in the domain of low stress triaxialities. This estimate is first compared with existing criteria. Moreover, its accuracy is assessed through comparison with results derived from numerical exact two-field criterion and with recently available numerical bounds.

Introduction

During the four past decades, the field of ductile damage has been the subject of important researches in non linear solids mechanics. Among the first to study the mechanical behavior of porous materials, McClintock [START_REF] Michel | The constitutive law of nonlinear viscous and porous materials[END_REF] and Rice and Tracey [START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF] derived, from a variational procedure, a cavity growth law for spherical and cylindrical voids in a von Mises matrix. These pioneering works were developed in the domain of dilute concentration of voids, which justifies an uncoupling between the plastic behavior of the matrix and the void growth. For finite porosities, Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF] developed a limit analysis approach of a hollow sphere made up of a rigid ideal plastic von Mises material and subjected to an arbitrary loading. In the case of homogeneous strain rate boundary conditions, Gurson derived a macroscopic yield function by considering two velocity fields: the first one corresponds to a homogeneous strain rate, while the second one takes into account the cavity expansion. For the latter, Gurson makes use of the exact solution of the hollow sphere subjected to a hydrostatic external loading. Consequently, the obtained macroscopic yield surface provides the exact solution in this particular loading case. Moreover, it was shown (see for instance [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF]) that the Gurson yield surface constitutes an upper bound for porous media, at least for Hashin composite-sphere assemblage. Later, Tvergaard [31] 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 [32]. Such extension, known as the GTN model, introduces three parameters, q 1 , q 2 and q 3 , which have to be determined 1 .

Alternatively, making use of variational techniques in non linear homogenization framework, Ponte-Castañeda [START_REF] Rice | On a ductile enlargement of voids in triaxial stress fields[END_REF] 2 obtained rigorous Hashin-Shtrikman type upper bounds for a von Mises plastic materials containing spherical voids. An important observation is that the Gurson model violates the non linear homogenization-based upper bound for low values of the stress triaxiality T = Σ m /Σ eq . However, the Gurson model predicts values more realistic than the Hashin-Shtrikman bound which yields unreasonably high values. Note that Danas and Ponte-Castañeda [START_REF] Danas | A homogenization-based constitutive model for isotropic viscoplastic porous media[END_REF] recently derived accurate yield surfaces for porous materials by using a second order homogenization method. A notable advantage of this second order approach is to exhibits an effect of the third invariant for isotropic microstructures. Note also that this effect has also be investigated in recent numerical studies by Thore et al. (2011). However, since one major advantage of the limit analysis approach over the latter is that it provides closed-form expression of the yield function, it would be interesting to improve the predictions of the original Gurson model. A possible way to do this consists to consider refined trial velocity fields still in the limit analysis framework. For instance, this was done by [START_REF] Gãrãjeu | Contribution à l'étude du comportement non linéaire de milieux poreux avec ou sans renfort[END_REF][START_REF] Gãrãjeu | A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids[END_REF] who considered the exact solution of the elastic hollow sphere subjected to an arbitrary loading. However this criterion still violates the Hashin-Shtrikman bounds.

The main objective of the present paper is to develop a limit analysis approach based on Eshelby-like velocity fields and to derive new expression of the yield function in the case of spherical voids. The paper is organized as follows. In section 2, we briefly recall the Gurson approach and state the trial velocity fields used in order to improve it (see subsection 2.2). Then, we discuss the question of boundary conditions associated to the consideration of the Eshelby-type velocity field and point out the requirement of a minimization procedure to derive the macroscopic yield surface. Section 3 is devoted to the determination of a closed-form expression of the macroscopic yield function. To this end, some approximations are introduced (see subsection 3.1). The minimization procedure is described in subsection 3.2. In order to underline the improvement of the new obtained criterion, some comparisons with other criteria (Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF] and Hashin-Shtrikman [START_REF] Rice | On a ductile enlargement of voids in triaxial stress fields[END_REF]) are provided in section 4.

Limit analysis approach and Eshelby-based velocity field

Statement of the Gurson approach

Since substantial modifications of the Gurson limit analysis will be proposed, we briefly recall here the main steps of the methodology which will be followed. In the following we consider the spherical basis (e r , e θ , e φ ) and the associated coordinates system (r, θ, φ), with θ ∈ [0, 2π] and φ ∈ [0, π]. Let us consider a hollow sphere with an external radius b and a void of radius a. This hollow sphere is subjected at its outer boundary to a homogeneous strain rate conditions:

v(r = b) = D.x (1)
The determination of an approximate macroscopic criterion, in the framework of limit analysis requires the choice of the velocity field into the solid matrix. As classically, and in agreement with [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF], we adopt the velocity in the following general form:

v = A.x + v E (2)
in which the term A.x is associated to a uniform strain rate A and v E is a heterogeneous field which corresponds to the expansion of the cavity. The matrix incompressibility implies that tr(A) = 0 and div(v E ) = 0. For v E , the choice made by Gurson corresponds to the solution of the plastic hollow sphere subjected to an hydrostatic loading; this field takes the form:

v E = C r 2 e r ( 3
)
in which C is a constant. The consideration of the boundary conditions leads to:

v = D.x + D m b 3 r 2 e r ( 4
)
where D m = 1 3 tr(D) and D represents the deviatoric part of D. It is assumed that the matrix material obeys to a von Mises criterion, with σ 0 the uniaxial tensile yield stress (material constant). Accordingly, the macroscopic plastic dissipation Π(D) is defined as:

Π(D) = σ 0 |Ω| ∫ Ω-ω d eq dV (5)
where d eq is the equivalent strain rate, d eq = √ 2 3 d : d. Tensor d is the deviatoric part of the strain rate tensor d deduced from the velocity field (4). The following inequality then holds for all macroscopic stresses, Σ and strain rate, D (Hill-Mandel lemma):

Σ : D ≤ Π(D) (6)
As a consequence, the macroscopic yield locus of the plastic porous medium is given by:

Σ = ∂Π ∂D (7)
In the case of von Mises matrix, the computation performed by Gurson, using (4), leads to:

Σ 2 eq σ 2 0 + 2f cosh ( 3 2 Σ m σ 0 ) -1 -f 2 ≤ 0 (8)
where Σ eq is the macroscopic von Mises equivalent stress, Σ eq = √ 3 2 Σ : Σ, Σ m = 1 3 tr(Σ) and f is the material porosity. As indicated before and proposed by [START_REF] Gãrãjeu | Contribution à l'étude du comportement non linéaire de milieux poreux avec ou sans renfort[END_REF], an improvement of the above yield function can be expected by investigating other velocity fields than the one used by Gurson.

Eshelby-based velocity field

Our proposal is to adopt, for v E , a field inspired from the solution of the Eshelby inclusion problem [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problem[END_REF], [START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF]. To this end, consider an infinite linear viscous medium in which a finite subdomain ω has a spherical shape with radius a. This subdomain is subjected to a uniform eigenstrain rate d * . From Eshelby's well-known results, the strain rate field is homogeneous inside the inclusion, while for an exterior point (outside the inclusion), it is heterogeneous. In the particular case of a spherical inclusion, the exterior point solution reads [START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF] (see also [START_REF] Pardoen | An extended model for void growth and coalescence[END_REF] or [START_REF] Leblond | Mécanique de la rupture fragile et ductile[END_REF]):

v E = a 5 5r 4 [ 5d * m 1 + 2 d * ] .e r + a 3 r 2 [ 1 - a 2 r 2 ]
d * rr e r [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF] in which d * m = 1 3 tr(d * ) and d * represents the deviatoric part of d * . In spherical coordinates, v E is expressed:

v E r = a 3 r 2 [( 1 - 3a 2 5r
It must be emphasized that the Eshelby-based velocity field allows to recover the one considered by Rice and Tracey [START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF] and by Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF] as a particular case associated to d * = 0 and d * m = C/a 2 in (9). This corresponds to the case of an inclusion subjected to uniform dilatational eigenstrain rate. Note that the deviatoric part of d * induces the inclusion shape change. Moreover, the resulting field, (2), corresponding to the Eshelby-like velocity field introduces 11 unknown kinematic parameters (the components of d * and A). The strain rate field d derived from the velocity field (2) with (9) reads:

d = A + d E ( 13
)
with d E = ∇ s v E , the symmetric part of the velocity gradient associated to v E . The components of the d E are given by:

d E rr = - 2a 3 r 3 d * rr + 4 5 a 5 r 5 [ 2d * rr -d * θθ -d * φφ ] d E φφ = a 3 r 3 d * rr + 1 5 a 5 r 5 [ -4d * rr + d * θθ + 3 * d * φφ ] d E θθ = a 3 r 3 d * rr + 1 5 a 5 r 5 [ -4d * rr + d * θθ + 3d * φφ ] d E rφ = a 3 r 3 [ 1 - 8 5
a 2 r 2 ] d * rφ ; d E rθ = a 3 r 3 [ 1 - 8 5
a 2 r 2 ] d * rθ ; d E φθ = 2 5 a 5 r 5 d * θφ ( 14
)
and can be put on the form:

d E = D : d * ( 15
)
where D is a fourth order tensor, which has the minor symmetries D ijkl = D jikl = D ijlk but not the major symmetry D ijkl ̸ = D klij . Components of tensor D are given by:

D rrrr = -2D θθrr = -2D φφrr = -2 a 3 r 3 [ 1 - 4 5 a 2 r 2 ] D θθθθ = D φφφφ = 3 5 a 5 r 5 D rrθθ = D rrφφ = - 4 5 a 5 r 5 D θθφφ = D φφθθ = D θφθφ = 1 5 a 5 r 5 D rθrθ = D rφrφ = 1 2 a 3 r 3 [ 1 - 8 5 a 2 r 2 ] (16)
Note that the components of the fourth order tensor D are such that D rrrr + D θθrr + D φφrr = 0, which corresponds to the matrix incompressibility.

Analysis of the associated boundary conditions

The verification of homogeneous strain rate boundary conditions (1), for the velocity field ( 2) with ( 9), implies:

A rr + f [( 1 - 3f 2/3 5 ) d * rr + 3f 2/3 5 d * m ] = D rr A rθ + 2f 5/3 5 d * rθ = D rθ ; A rφ + 2f 5/3 5 d * rφ = D rφ (17)
Introducing the deviatoric and mean part of A, d * and D into (17

) (i.e. A = Ā, d * = d * m 1 + d * , D = D m 1 + D) gives: f d * m = D m A rr + f ( 1 - 3f 2/3 5 ) d * rr = Drr A rθ + 2f 5/3 5 d * rθ = Drθ ; A rφ + 2f 5/3 5 d * rφ = Drφ (18)
Solution of ( 18) is:

d * m = D m f ; A = D; d * = 0 ( 19
)
It appears that the considered class of Eshelby-based velocity fields ((2) with ( 9)) complies with the homogeneous strain rate boundary conditions if and only if d * = 0. This is equivalent to d * rr = constant. The velocity field (2) with (9) reduces then to (4) used by Gurson. It follows that all supplementary velocity fields (associated to d * ̸ = 0) cannot comply with the macroscopic strain rate tensor D. In fact, the consideration of the Eshelby-based velocity field can be rigorously justified in the context of limit analysis only for homogeneous stress boundary conditions, σ.n = Σ.n. Then, and in agreement with Hill-Mandel lemma, the macroscopic strain rate D is related to the local strain rate, d, by the average rule:

D = 1 |Ω| ∫ Ω ddV ( 20
)
Note that, the use of the kinematic approach (developed in section 3) leads to an upper bound of the macroscopic yield surface in the context of the homogeneous stress boundary conditions. Introducing now ( 13) with ( 15) into (20) leads to:

D = 1 |Ω| ∫ Ω ddV = A + f [ d * m + 2 5 d * ] (21)
which links the macroscopic strain rate D to A and d * . Consideration of the deviatoric and the mean part of D gives:

         D = A + 2f 5 d * D m = f d * m ( 22
)
In these relations, one can observe that the mean part of the tensor d * is known and explicitly given in terms of the macroscopic strain rate, D m . Moreover [START_REF] Ponte-Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF] provides 6 scalar relations for the 11 unknown parameters (the components of A and d * ). It still remains 5 parameters, namely the components of d * , whose determination will imply a minimization procedure.

3 The two-field estimate of the macroscopic yield criterion

The approximate macroscopic dissipation

This section is devoted to the determination of an approximate criterion by considering the Eshelby-based trial velocity field. Although the general methodology has been depicted in the section 2.1, it must be emphasized that our approach presents some substantial differences with the Gurson one. Indeed, as mentioned before, the macroscopic plastic dissipation will be obtained after a minimization procedure on the remaining unknown parameters:

Π(D) = min d * [ Γ(D, d * ) ] (23)
with:

Γ(D, d * ) = σ 0 |Ω| ∫ Ω-ω d eq dV ( 24
)
in which

d 2 eq = A 2 eq + 4 3 A : D : d * + 2 3 d * : D T : D : d * ( 25
)
We now propose to derive an explicit expression of Γ(D, d) which requires the integration of the equivalent strain rate d eq as given by [START_REF] Suquet | On bounds for the overall potential of power law materials containing voids with an arbitrary shape[END_REF]. Let us introduce the first approximation used to perform the integration of the local dissipation:

A1 : The mean value of d eq over the unit sphere, i.e. < d eq > S(r) , is replaced by

√ < d 2 eq > S(r)
, such as:

< d 2 eq > S(r) = A 2 eq + 4 3 A :< D > S(r) : d * + 2 3 d * :< D T : D > S(r) : d * (26)
Note that this approximation has been used by [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I. -Yield criterion and flow rules for porous ductile media[END_REF] and has the advantage to preserve the upper bound character of the derived macroscopic yield function.

As the mean value of D over the unit sphere is zero, < D > S(r) = 0, the crossed term is null. Concerning the term < D T : D > S(r) , its computation gives:

< D T : D > S(r) = { 2J + 6 25 G(u)K } u 2 (27)
with J = 1 3 1 ⊗ 1, 1 being the second order unit tensor. Tensor K is defined by K = I -J, with I the fourth order symmetric identity tensor. G(u) is a function of variable u = a 3 /r 3 , defined by:

G(u) = 1 3 ( 15 -40u 2/3 + 28u 4/3 ) (28)
To summarize, < d 2 eq > S(r) is expressed as:

< d 2 eq > S(r) = A 2 eq + [ 4d * m 2 + 6G(u) 25 d * eq 2 ] u 2 (29)
With approximation A1, Γ(D, d) reads:

Γ(D, d) = σ 0 f ∫ u=1 u=f A 2 eq + ( 4d * m 2 + 6G(u) 25 d * eq 2 ) u 2 du u 2 (30)
We aim now at deriving a closed-form expression of Γ(D, d). Due to the presence of G(u) in (30), a second approximation is needed, the objective being to put Γ(D, d) in the form σ 0 f

∫ 1 f √ A 2 + B 2 u 2 du u 2 (
which can be easily computed) where A and B are two constants. This second approximation is expressed as:

A2 : function G(u) is replaced by its mean value, denoted thereafter g(f ), over the interval [u 1 , u 2 ] = [1, f ]
The computation of g(f ) leads to:

g(f ) = 1 1 -f ∫ 1 f G(u)du = 1 -4f (1 -f 2/3 ) 2 1 -f (31)
from which:

Γ(D, d * ) = σ 0 f ∫ 1 f { A 2 eq + B 2 u 2 } 1/2 du u 2 = σ 0 f   B arcsinh { uB A eq } - √ A 2 eq + u 2 B 2 u   1 f ( 32
)
with

B 2 = 4d * m 2 + 6g(f ) 25 d * eq 2 = 4 f 2 D 2 m + 6g(f ) 25 d * eq 2 (33)
Note that, in this form, Γ(D, d * ) is equivalent to the one obtained by Gurson in the case d * eq = 0 and A eq = D eq . It is important to emphasize that, since any inequality can be established, A2 introduces an "uncontrolled" approximation in the computation of Γ(D, d * ). Consequently, at the difference of Gurson result, the upper bound character of the approximate criterion which will be established in the following is not ensured. Note however that such upper bound can be numerically determined starting from (34) in which Γ(D, d * ) is obtained from a numerical integration.

Minimization procedure and determination of the macroscopic criterion

The main objective of the present section is to determine an expression of the macroscopic yield function from the procedure depicted in section (3.1). This requires first the determination of the unknown parameters d * from the minimization of Γ(D, d * ). Once this minimization is performed, it still remains to determine the macroscopic yield surface from [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids -case of axisymmetric prolate ellipsoidal cavities[END_REF]. However, it was possible to perform these two steps at the same time as3 :

Σ = ∂Γ(D, d * ) ∂D with: ∂Γ(D, d * ) ∂ d * = 0 (34)
A complete resolution of this non linear system can be done by introducing into (34) the following change of variables: Γ(D, d * ) = Γ(A eq , B). In one hand, the macroscopic stress tensor Σ is given by:

Σ = ∂Γ ∂A eq ∂A eq ∂D + ∂Γ ∂B ∂B ∂D (35)
In another hand, in (34), the second relation reads:

∂Γ ∂A eq ∂A eq ∂ d * + ∂Γ ∂B ∂B ∂ d * = 0 ( 36
)
Derivatives of A eq and B according to D and d * provide:

∂A eq ∂D = 2 3 A A eq ; ∂A eq ∂ d * = -f 4 15 A A eq ∂B ∂D = 4D m 3f 2 B 1; ∂B ∂ d * = 4g(f ) 25B d * (37)
Introducing expressions in (37) into ( 35) and (36) gives the following expressions for the mean and equivalent macroscopic stress, as well as the minimum condition:

Σ m = 4 3 ∂Γ ∂B D m f 2 B Σ eq = ∂Γ ∂A eq ∂Γ ∂A eq = 3g(f ) 5f d * eq B ∂Γ ∂B (38)
The elimination of B between the first and the last equation in (38) leads to:

∂Γ ∂A eq = Σ eq ∂Γ ∂B = f √ 9 4 Σ 2 m + 2 3g(f ) Σ 2 eq ( 39
)
The computation of ∂Γ ∂A eq and ∂Γ ∂B is performed by considering (32):

∂Γ ∂A eq = σ 0 [ √ 1 + f 2 ξ 2 -f √ 1 + ξ 2 ] ∂Γ ∂B = σ 0 f [ arcsinh(ξ) -arcsinh(f ξ) ] (40)
with ξ = B/A eq , B being given by (33).

The macroscopic criterion is obtained by replacing the derivatives of Γ in (40) by (39) and eliminating ξ in (40):

F (Σ, f ) = Σ 2 eq σ 2 0 + 2f cosh   9 4 Σ 2 m σ 2 0 + 2 3g(f ) Σ 2 eq σ 2 0   -1 -f 2 ≤ 0 (41)
which is the most important practical result of the present paper. In the case of low values of the porosity, one can consider for simplicity that g(f ) ≃ 1 which leads to:

F (Σ, f ) ≃ Σ 2 eq σ 2 0 + 2f cosh   9 4 Σ 2 m σ 2 0 + 2 3 Σ 2 eq σ 2 0   -1 -f 2 ≤ 0 ( 42
)
Remark: it must be emphasized that the closed form expression of the macroscopic criterion depends on the different approximations considered for the derivation of Γ(D, d). Indeed, consideration of alternative approximation instead of A2 (see section 3.1) may lead to a criterion different of (41). This has been done and reported in appendix B.

Evaluation and validation of the established criterion

We aim now at comparing the new macroscopic yield function (41) with existing criteria. Moreover, the accuracy of this criterion will be evaluated by its comparison with recent numerical bounds.

Comparison with existing criteria

The main difference of the new criterion (41) with ( 8) lies in the presence of the term introduced by the equivalent stress and the porosity into the hyperbolic cosine. However, the new criterion still preserves, as Gurson's one, the exact solution of the hollow sphere subjected to an hydrostatic loading, Σ m = -2σ 0 3 ln(f ). This is a consequence of the fact that the velocity field considered in the present study contains the one already used by Gurson. Moreover, it is interesting to compare the obtained results to the one deduced from non linear variational homogenization techniques proposed by Ponte Castañeda [START_REF] Rice | On a ductile enlargement of voids in triaxial stress fields[END_REF], [START_REF] Suquet | On bounds for the overall potential of power law materials containing voids with an arbitrary shape[END_REF] and Suquet [29]. Indeed, for a representative elementary volume of a porous medium with spherical voids, these authors established a macroscopic elliptic criterion which reads [START_REF] Rice | On a ductile enlargement of voids in triaxial stress fields[END_REF]:

( 1 + 2f 3 ) Σ 2 eq σ 2 0 + 9f 4 Σ 2 m σ 2 0 -(1 -f ) 2 ≤ 0 (43)
and constitutes an upper bound for the yield surface of the porous medium.

It is convenient to recall that the Gurson criterion violates this upper bound for low values of the stress triaxiality T = Σ m /Σ eq . This has motivated modifications such as the ones proposed by Tvergaard [31] or [START_REF] Monchiet | An improvement of Gurson-type models of porous materials by using Eshelbylike trial velocity fields[END_REF] and widely used in literature for numerical applications. However, due to the heuristic nature of these modifications, they are not considered in the following comparative analysis.

The results are shown on figure 1 to 3 for three values of porosity, f = 0.01, f = 0.1 and f = 0.2. On these figures are also reported data numerically obtained by computing the exact value of Γ(D, d * ) given by ( 24) without any approximation. For this "exact" solution, the limit analysis of the hollow sphere is performed by the consideration of the Eshelby type velocity field. So, one has to compute Γ(D, d * ) by integrating the local dissipation [START_REF] Rice | On a ductile enlargement of voids in triaxial stress fields[END_REF] and to minimize Γ(D, d * ) according to the unknown kinematical parameters, d * . In a second time, the macroscopic criterion is obtained from [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids -case of axisymmetric prolate ellipsoidal cavities[END_REF]. These numerical results validate the approximations made for the derivation of the closed-form criterion (41).

We compare now the new criterion to the one obtained from variational homogenization techniques and to the Gurson one. Note first that for low values of porosity (see figure 1), the new criterion seems to coincide with the Gurson one. The comparison also shows that for higher porosities, the established criterion provides a significant modification of the Gurson criterion for low triaxialities. Moreover, it appears that the new criterion is close to the one predicted by the variational technique for low stress triaxialities while it is interestingly in agreement with the Gurson criterion at high triaxialities. The predictions of the different criteria for pure deviatoric macroscopic stress states (triaxiality equal to 0) are illustrated on figure [START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF] for various values of porosities. For completeness, it must be emphasized that the very recent works by Danas et al. [START_REF] Danas | A homogenization-based constitutive model for isotropic viscoplastic porous media[END_REF], using a second order method, lead to significant improvements of the predictions provided by non linear homogenization techniques. Note however that the latter do not deliver a closed-form expression of the macroscopic yield function.

Comparison to numerical limit analysis results

As mentioned in section 2.3, boundary conditions used for the derivation of the closed form criterion can be interpreted as homogeneous stress boundary conditions. The determination of the macroscopic criterion has been made by using a kinematical approach and so leads to an upper bound of the macroscopic yield surface for these stress boundary conditions. However to obtain the closed form expression (41), two approximations have been done. The first one, A1, already used by Gurson, has the advantage to preserve the upper bound character. But, the second approximation (A2) is "uncontrolled". In order to assess the accuracy of this approximation, we now propose, in table 1, 2 and 3, some comparisons with numerical solutions recently reported by Trillat and Pastor [30]. These comparisons are made for the following values of the porosity: f = 0.01, f = 0.064 and f = 0.1. The numerical data given in table 1, corresponding to the porosity f = 0.1, are reported on figure 5. In each table, are given the ratio between the macroscopic equivalent stress and the yield shear stress, namely Σ eq /τ 0 (with σ 0 = τ 0 √ 3) for a number of fixed values of the ratio Σ m /τ 0 . The last line in each table provides the numerical value of the ratio Σ max m /τ 0 corresponding to Σ m /τ 0 associated to Σ eq = 0 (purely hydrostatic macroscopic stress). In the first column, are given the prediction of the new approximate criterion (41) whose accuracy is first compared to the "exact"numerical solution without any approximation. These numerical data are reported in the second column of each table. It can be first observed that, for the three porosities, and whatever the value of the stress triaxiality T = Σ m /Σ eq , the approximated criterion (41) is close to the exact two-field criterion. Moreover, it appears that, despite the approximations required for its derivation, (41) preserves (for the problem with homogeneous stress boundary conditions) the upper bound character of the approach. To evaluate the accuracy of the above results, we now propose some comparisons with data provided by Trillat and Pastor [30] and numerically obtained by means of Finite Element-based limit analysis of the hollow sphere combined with convex optimization. In [30], the unit cell is subjected either to a homogeneous strain rate conditions (v = D.x) or, as in the present study, to homogeneous stress boundary conditions (σ.n = Σ.n). For each boundary conditions, the problem is solved by means of a static approach (leading to a lower bound for the macroscopic criterion) and a kinematical approach (giving an upper bound). Tables 1, 2 and 3 give the numerical results corresponding to homogeneous stress boundary conditions (BC) and strain rate boundary conditions. It is generally observed that, comparatively to homogeneous stress boundary conditions, the values obtained in the present study are larger than the upper bound provided by the Finite Element bounds. This can also be clearly observed on figure 5 For sake of completeness, it can be useful to compare in table 1 to 3 our results to the ones given by the Finite Element bounds in the case of homogeneous strain rate boundary conditions. Interestingly, it appears that our results compare well to these numerical bounds.

Conclusion

A micromechanics-based modification of the classical Gurson criterion for porous media has been derived. The approach has been developed for a hollow sphere in the framework of limit analysis. The originality of the present study lies in the consideration of new trial velocity fields inspired from the Eshelby exterior point solution to inclusions problem in which the eigenstrains are unknown. Due to the presence of these unknowns, the new estimate of the macroscopic yield criterion has been obtained by implementing a minimization procedure. It has been shown that the derived criterion exhibits some remarkable properties:

• the new criterion is of Gurson type (i.e. non elliptic), but has the significant difference to contain the equivalent stress, Σ eq , into the hyperbolic cosine while preserving the accuracy of the Gurson model for high stress triaxialities. In particular, as for the Gurson model, the new criterion allows to retrieve the exact solution for hydrostatic loadings. • the most important modification of the Gurson criterion appears for low stress triaxialities, where the derived yield surface is close to the Hashin-Shtrikman bound established by Ponte-castañeda [START_REF] Rice | On a ductile enlargement of voids in triaxial stress fields[END_REF] in the framework of the non-linear homogenization applied to a representative elementary volume of isotropic porous media.

In addition to these general remarks, the good accuracy of the various approximations used to obtain the closed-form expression of the yield function is shown by means of comparison with the numerical exact two-field criterion. Moreover, complete evaluation of the criterion is provided by comparisons with recent numerical upper and lower bounds reported by [30]. Both homogeneous strain rate boundary conditions and homogeneous stress boundary conditions are considered in these numerical computations. It is observed that the obtained estimate of the criterion provides yield surfaces close to the numerical results with homogeneous strain boundary conditions. It must be recalled that the Eshelby-based velocity fields do not comply with homogeneous strain rate boundary conditions for the hollow sphere problem. Consequently, the obtained result must be rigorously interpreted in terms of homogeneous stress boundary conditions. A first step of generalization of the above results in order to formulate a constitutive model of isotropic porous media has been proposed in [START_REF] Mcclintock | A criterion for ductile fracture by the growth of holes[END_REF] in which the derived model (and its predictions) is compared to the standard Gurson model and the GTN model. Moreover, a generalization of the present work to anisotropic microstructures has been already introduced in Monchiet et al. [START_REF] Monchiet | Macroscopic yield criteria for plastic anisotropic materials containing spheroidal voids[END_REF] by taking into account voids shape effects. Investigation of ductile porous media made up of a von Mises solid matrix and spheroidal (prolate or oblate) cavities is the subject of a forthcoming paper [START_REF] Mura | Micromechanics of Defects in Solids[END_REF].

Σ eq A A eq = 3g(f ) 5 Σ b B d * (A.1)
which can also be put on the form:

A A eq = d * d * eq (A.2)
where d * eq is given by:

d * eq B = 5B 3g(f ) Σ eq Σ b (A.3)
In the last relation B is also function of d * eq (see equation ( 33)), the computation of d * eq gives:

d * eq = 20 9
Σ eq Σ m D m (A.

B Alternative expression of the macroscopic yield function

We aim at presenting here an alternative form of the macroscopic yield function of the porous medium. To this end, instead of approximation A2 (see section 3.1), we use in (30) giving Γ(D, d) the following alternative one:

A2 * : function 6G(u)u 2 /25 is replaced by a constant denoted h(f ). Γ(D, d) is then approximated by: Γ(D, d) = σ 0 f ∫ u=1 u=f √ A 2 eq + h(f )d * eq 2 + 4d * m 2 u 2 du u 2 (B.1)
The approximation consisting to replace 6G(u)u 2 /25 by its average value over the interval [u1, u2] = [f, 1] appears to be poor here. A better approximation has be found by introducing the change of variable w = f /u in (B.1) :

Γ(D, d) = σ 0 ∫ w=1 w=f √ A 2 eq + H(w)d * eq 2 + 4D 2 m w 2 dw (B.2)
with:

H(w) = 2 25 ( 15 -40(f /w) 2/3 + 28(f /w) 4/3 ) (f /w) 2 (B.3)
and by replacing H(w) by its average over the interval [w1, w2] = [f, 1]. This leads to the following expression for h(f ):

h(f ) = 1 1 -f ∫ 1 f H(w)dw = 6f g(f ) 25 (B.4)
where g(f ) is defined by (31). Now, the integration of (B.2) yields: which appears to be more simple than (A.6). It follows that A is given by: Putting (B.9) in (B.5) yields the macroscopic dissipation Π(D). Recalling that the macroscopic dissipation is the one derived by Gurson with A eq replaced by A, the macroscopic criterion is obtained by a change of variable:

Γ(D, d * ) = σ 0 ∫ 1 f √ A 2 + 4D 2 m w 2 dw = σ 0 [ B arcsinh ( 2D m Aw ) - √ A 2 w 2 + B 2
A 2 = h(f )
(

1 + 2 3 f g(f ) ) Σ 2 eq σ 2 0 + 2f cosh ( 3 2 Σ m σ 0 ) -1 -f 2 ≤ 0 (B.10)
It is worth noticing that, for low porosities (f ≃ 0):

1 + 2 3 f g(f ) ≃ 1 + 2 3 f + o(f 2 ) (B.11)
which leads to the following approximate expression of the macroscopic yield function for low values of porosities:

(

1 + 2 3 f ) Σ 2 eq σ 2 0 + 2f cosh ( 3 2 Σ m σ 0 ) -1 -f 2 ≤ 0 (B.12)
Remarkably, this expression is exactly the one heuristically proposed by [START_REF] Lin | Numerical implementation of a recent improved Gurson-type model and application to ductile fracture[END_REF].

It must be emphasized that the alternative criterion (B.10), still preserve the exact solution for purely hydrostatic loading, i.e Σ m = -2σ 0 3 ln(f ). In the case of purely deviatoric macroscopic stress states (Σ m = 0), we represent on figure B.1 the value of the equivalent stress as function of the porosity for the two approximate criterion (41) and (B.10) and the numerical two-field criterion. It can be observed that the predictions of the two approximate criteria (41) and (B.10) are equivalent. Also for comparison purpose, the predictions of the alternative criterion, (B.10), and the previously proposed one, (42), are plotted on figure B.2 for a porosity f = 0.1. On this figure are also reported the data numerically obtained for the exact two-field solution. Despite the high proximity of the results, it can be observed that the first approximate form (42) is slightly more close to the numerical points than the second form (B.2). expression (42) of the criterion is then preferred.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Comparison between the new criterion (41), Gurson's model (8), the Hashin-Shtrikman bound (43) and the numerical exact solution (obtained by using Eshelby-type velocity fields). The porosity has been chosen equal to f = 0.01.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison between the new criterion (41), Gurson's model (8), the Hashin-Shtrikman bound (43) and the numerical exact solution (obtained by using Eshelby-type velocity fields). The porosity has been chosen equal to f = 0.2.

Fig. 5 .

 5 Fig. 5. Comparison between the new criterion (41), numerical results provided by Trillat and Pastor [30], and exact solution, obtained numerically, of the problem depicted in section 3.1 for the porosity f = 0.1.

4 )

 4 Let us recall that D = A + 2f 5 d * , so from (A.2) it appears that the directions of A, d * and D are the same:

2

 2 eq + h(f )(d * eq )2 (B.6)At this stage, it is interesting to notice that the macroscopic dissipation is exactly the one derived by Gurson, A eq being replaced here by A whose expression has to be established now. As Γ(D, d * ) only depends on d * through the scalar A, the minimum of Γ(D, d * ) according to d * is obtained for: see equation (22)), it can be shown that the minimum is attained for:

Fig. B. 1 .

 1 Fig. B.1. Comparison between the approximate criteria (41), (B.10) and numerical data for the exact two field criterion. Value of the equivalent stress as function of the porosity.

  Fig. B.2. Comparison between the approximate criteria (41), (B.10) and numerical data for the exact two field criterion. The porosity has been chosen equal to f = 0.1.

Table 1 :

 1 . numerical values of Σ eq /τ 0 versus Σ m /τ 0 corresponding to f = 0.01 Σ m /σ 0

	Σ m /σ 0	Present approach	Trillat and Pastor's results Trillat and Pastor's results:
					stress BC		strain rate BC
		approx.	exact	static	kinematic	static	kinematic
		criterion	2 fields
			criterion
	0	1.7084		1.6937	1.7048	1.7	1.7086
	1.6	1.6867	1.6817	1.6732	1.6835	1.6784	1.6894
	3.2	1.5711	1.5631	1.5016	1.5511	1.5435	1.5805
	4	1.4058	1.3893	1.2745	1.3713	1.3452	1.4271
	4.4	1.2564	1.2290	1.0705	1.2214	1.1568	1.2888
	4.8	1.0141	0.9674	0.7698	1.0009	0.845	1.067
	5	0.8245	0.7636	0.5472	0.8467	0.5467	0.9013
	5.1	0.6954	0.6150	0.2001	0.7509	0.3856	0.7966
	5.2	0.5211	0.4358	-	0.6355	-	0.6698
	5.3	0.2055	0.1649	-	0.4872	-	0.5084
	Σ max m /τ 0 5.3176	5.3175	5.1791	5.4385	5.1746	5.4388

Table 2 :

 2 numerical values of Σ eq /τ 0 versus Σ m /τ 0 corresponding to f = 0.064 Σ m /σ 0

		Present approach	Trillat and Pastor's results Trillat and Pastor's results:
					stress BC	strain rate BC
		approx.	exact	static	kinematic	static	kinematic
		criterion	2 fields
			criterion
	0	1.4882	1.4736	1.3913	1.4213	1.5020	1.5206
				(1.3774)	(1.4150)	(1.4805)	(1.5206)
	0.8	1.4368	1.4113	1.3621	1.3853	1.4598	1.4774
				(1.3430)	(1.3849)	(1.4624)	(1.4774)
	1.6	1.2473	1.2297	1.1557	1.1875	1.2699	1.2869
				(1.1550)	(1.1872)	(1.2669)	(1.2868)
	2	1.0562	1.0451	0.9490	0.9913	1.0373	1.0820
				(0.9436)	(0.9912)	(1.0373)	(1.0820)
	2.2	0.9136	0.9038	0.8039	0.8511	0.8729	0.9269
				(0.7918)	(0.8511)	(0.8598)	(0.9269)
	2.4	0.7113	0.6994	0.5789	0.6608	0.6347	0.7108
				(0.5973)	(0.6608)	(0.6315)	(0.7106)
	2.5	0.5674	0.5525	0.4504	0.5305	0.4300	0.5632
				(0.4318)	(0.5305)	(0.4621)	(0.5601)
	Σ max m /τ 0 2.6587	2.6587	2.6239	2.6705	2.6316	2.6707
				(2.6229)	(2.6704)	(2.5650)	(2.6316)

Table 3 :

 3 numerical values of Σ eq /τ 0 versus Σ m /τ 0 corresponding to f = 0.1.

	The values within parenthesis correspond to the initial results in [30]
	(Numerical results without parenthesis are recent results kindly
	communicated by the J. Pastor).

For completeness, mention must be made of[START_REF] Lin | Numerical implementation of a recent improved Gurson-type model and application to ductile fracture[END_REF] and of several extensions of the Gurson model accounting for voids shape effects (see for instance[START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids -case of axisymmetric prolate ellipsoidal cavities[END_REF],[START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids -case of axisymmetric oblate ellipsoidal cavities[END_REF],[START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF],[START_REF] Gãrãjeu | A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids[END_REF],[START_REF] Ponte-Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF],[START_REF] Benzerga | Anisotropic ductile fracture. Part II: Theory[END_REF],[START_REF] Monchiet | Macroscopic criteria for ductile materials containing spheroidal voids: an Eshelby-like velocity fields approach[END_REF])

A similar derivation has been done bySuquet [29] for cylindrical voids.

Note that the computation of d * has been done and the result is reported in appendix A.