domain: spi nlin spi.mat

text:
Long Cheng

Yun Jia

Abdelbacet Oueslati

Géry De Saxcé

Djimedo Kondo

A bipotential-based limit analysis and homogenization of ductile porous materials with non-associated Drucker-Prager matrix

Keywords: Bipotential theory, Nonlinear Homogenization, Extended Limit analysis, Ductile porous materials, Non associated plasticity, Drucker-Prager matrix

come

A bipotential-based limit analysis and homogenization of ductile porous materials with non-associated Drucker-Prager matrix

Introduction

In his famous paper, [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth -part I: Yield criteria and flow rules for porous ductile media[END_REF] proposed an upper bound limit analysis approach of a hollow sphere and a hollow cylinder having a von Mises solid matrix. Several extensions of • in soil mechanics, non-associated Drucker-Prager (de Saxcé , 1993;[START_REF] Berga | Elastoplastic Finite Element Analysis of Soil Problems with Implicit Standard Material Constitutive Laws[END_REF]de Saxcé , 1998a;[START_REF] Bousshine | Softening in stress-strain curve for Drucker-Prager nonassociated plasticity[END_REF][START_REF] Hjiaj | A complete stress update algorithm for the non-associated Drucker-Prager model including treatment of the apex[END_REF] and Cam-Clay models [START_REF] De Saxcé | The bipotential method, a new variational and numerical treatment of the dissipative laws of materials[END_REF][START_REF] Zouain | Plastic collapse in non-associated hardening materials with application to Cam-clay[END_REF][START_REF] Zouain | Potentials for the modified Cam-Clay model[END_REF],

• the non linear kinematical hardening rule for cyclic Plasticity [START_REF] De Saxcé | Une généralisation de l'inégalité de Fenchel et ses applications aux lois constitutives[END_REF][START_REF] Bodovillé | Plasticity with non linear kinematic hardening : modelling and shakedown analysis by the bipotential approach[END_REF] and Viscoplasticity [START_REF] Hjiaj | Matériaux viscoplastiques et loi de normalité implicites[END_REF][START_REF] Magnier | Bipotential Versus Return Mapping Algorithms: Implementation of Non-Associated Flow Rules[END_REF][START_REF] Bouby | A comparison between analytical calculations of the shakedown load by the bipotentiel approach and step-by-step computations for elastoplastic materials with nonlinear kinematic hardening[END_REF][START_REF] Bouby | Shakedown analysis: comparison between models with the linear unlimited, linear limited and non linear kinematic hardening[END_REF],

• Lemaitre's coupled plasticity-damage law [START_REF] Bodovillé | On damage and implicit standard materials[END_REF],

• the coaxial laws (de Saxcé , 2002;[START_REF] Vallée | Hill's bipotential[END_REF],

• Coulomb's friction law (de Saxcé , 1998b[START_REF] De Saxcé | Une généralisation de l'inégalité de Fenchel et ses applications aux lois constitutives[END_REF][START_REF] De Saxcé | On the extension of limit analysis theorems to the non-associated flow rules in soils and to the contact with Coulomb's friction[END_REF](de Saxcé , , 1998b,a;,a;[START_REF] Bousshine | Plastic limit load of plane frames with frictional contact supports[END_REF][START_REF] Hjiaj | A variational-inequality based formulation of the frictional contact law with a non-associated sliding rule[END_REF][START_REF] Hjiaj | Three dimensional finite element computations for frictional contact problems with on-associated sliding rule[END_REF]Feng et al., 2006b,a;[START_REF] Fortin | Modlisation numrique des milieux granulaires par lapproche du bipotentiel[END_REF][START_REF] Fortin | An improved discrete element method based on a variational formulation of the frictional contact law[END_REF][START_REF] Laborde | Fixed points strategies for elastostatic frictional contact problems[END_REF],

• the blurred constitutive laws (Buliga et al., 2009b(Buliga et al., , 2010b)).

A complete survey of the bipotential approach can be found in de Saxcé (2002). In the previous works, robust numerical algorithms were proposed to solve structural mechanics problems.

Coming back to the Limit Analysis let us say that a general method to determine the plastic collapse of structures under proportional loading [START_REF] Suquet | Plasticité et homogénéisation[END_REF]Salençon , 1983;[START_REF] Save | Plastic limit analysis of plates, shells and disks[END_REF], even particular in soil mechanics [START_REF] Chen | Limit Analysis and Soil Plasticity[END_REF][START_REF] Chen | Plastic limit state of the hollow sphere model with non-associated Drucker-Prager material under isotropic loading[END_REF], but it is restricted to associated plasticity, then with normality law. The classical presentation of the non-associated plasticity is based on a yield function and a plastic potential. The bipotential offers an alternative formulation which naturally opening to a variational formulation, and then paving the way for an extension of limit analysis techniques to non-associated laws (de Saxcé , 1998a;[START_REF] Bousshine | Softening in stress-strain curve for Drucker-Prager nonassociated plasticity[END_REF][START_REF] Bousshine | Plastic limit load of plane frames with frictional contact supports[END_REF][START_REF] Chaaba | Kinematic Limit Analysis of Nonassociated Perfectly Plastic Material by the Bipotential Approach and Finite Element Method[END_REF][START_REF] Zouain | Plastic collapse in non-associated hardening materials with application to Cam-clay[END_REF]. Extension of limit analysis theory to the repeated variable loading, known as shakedown theory, has been successfully generalized to the ISM1 by the bipotential approach in (de Saxcé , 2002;[START_REF] Bousshine | Softening in stress-strain curve for Drucker-Prager nonassociated plasticity[END_REF][START_REF] Bousshine | A new approach to shakedown analysis for non-standard elastoplastic material by the bipotential[END_REF][START_REF] Bouby | A comparison between analytical calculations of the shakedown load by the bipotentiel approach and step-by-step computations for elastoplastic materials with nonlinear kinematic hardening[END_REF][START_REF] Bouby | Shakedown analysis: comparison between models with the linear unlimited, linear limited and non linear kinematic hardening[END_REF].

The aim of the present study is to formulate a macroscopic model for "ductile porous materials with a non-associated Drucker-Prager"-type matrix, using homogenization techniques combined with the bipotential theory. The paper is organized as follows: the non-associated Drucker-Prager plastic model, for which the yield criterion and plastic potential are respectively defined by two functions, is first summarized in Section 2. Next, we introduce in Section 3, the bipotential theory and its two dual fields (stress and velocity fields) based formulation, which allows us to derive the plastic criterion and the non-associated flow rule. An application of the bipotential theory to the non-associated Drucker-Prager plastic model is particularly discussed in subsection 3.3. Section 4 is devoted to the bipotential-based extended limit analysis approach of non-associated porous media. The proposed formulation provides a fundamental variational theory for the macroscopic modeling of a large class of porous media. In Section 5, the proposed bipotential-based theory is implemented in the case of a hollow sphere having a rigid perfectly plastic matrix obeying to a non-associated Drucker-Prager flow rule. This implementation will be performed by adopting simple trial stress and velocity fields. This allows to derive in subsection 5.3 a closed-form expression of the macroscopic criterion and the non-associated plastic flow rule. Furthermore, some special cases, corresponding to existing models previously proposed in literature, are discussed in subsection 6. Finally, in Section 7, the established macroscopic criterion, flow rule and void evolution are respectively assessed and validated by comparison with Finite Element solutions. Drucker-Prager model (Fig. 1) requires the consideration of a yield criterion in the form:

Brief recall the non-associated Drucker-Prager model

F (σ) = σ e + 3ασ m -σ 0 ≤ 0 , (1)
where σ e is the equivalent stress, σ m the mean stress, σ 0 > 0 the shear cohesion stress of the material and α the pressure sensitivity factor related to the friction angle φ by: tan φ = 3α .

Let us introduce the plastic potential:

G(σ) = σ e + 3βσ m (2)
where β (β ≤ α) depends on the dilatancy angle ψ through:

tan ψ = 3β .

Except for the apex of Drucker-Prager cone (σ e = 0, σ m = σ 0 / 3α ) where σ e is not differentiable, the plastic strain rate is given by the non-associated yielding rule:

d = d eq ∂G ∂σ = d eq 3s 2σ e + β1 , (3)
where σ is Cauchy stress tensor, s the deviatoric stress, 1 the unit tensor. d eq =| 2 3 d : d | 1/2 with d being the deviatoric part of d. The plastic dilatancy reads:

d m = 1 3 tr d = βd eq (4)
This suggests to introduce:

H(d) = βd eq -d m
The plastic yielding rule ( 3) is completed at the apex by the admissibility condition:

H(d) ≤ 0
while, because of (4), H(d) = 0 at the other points of the yielding surface (called regular points). Of course, for the particular event ψ = φ, the normality rule is recovered and the plasticity model is associated. Without loss of generality, we can assume that:

0 ≤ β ≤ α < 1 2 , (5)
or equivalently 0 ≤ ψ ≤ φ < 56 • 18 . In practice, these conditions are fulfilled by the geomaterials and other pressure sensitive dilatant materials. Some examples of experimental data concerning the friction angle can be found for polymers, high strength steels and aluminium in [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF].

Bipotential-based formulation of constitutive models

Rigid perfectly plastic model is usually considered to obtain the analytical solution (plastic criterion and potential) for a large class of materials, simultaneously by adopting the Limit Analysis approach, which is extensively discussed in literature . However, this conventional approach can only be rigorously used for Generalized Standard Materials (GSM), that is materials which obey a normality law. the standard limit analysis framework is then not suitable for materials which obey to a non-associated flow rule (for instance, geomaterials). This question has been discussed in several works (Salençon , 1983;[START_REF] Drucker | Limit analysis of two -and three -dimensional soil mechanics problems[END_REF]Palmer , 1973;[START_REF] Radenkovic | Limit analysis theorems for a Coulomb material with a non standard dilatation[END_REF][START_REF] Telega | Extremum principles for nonpotential and initial-value problems[END_REF][START_REF] Telega | An of Limit Analysis Theorems to Incompressible Material with a Non-Associated Flow Rule[END_REF].

In order to overcome this problem, de Saxcé et al. has proposed in previous papers a new modeling of the non-associated constitutive laws based on the concept of bipotential [START_REF] De Saxcé | New inequality and functional for contact friction: The implicit standard material approach[END_REF][START_REF] De Saxcé | Une généralisation de l'inégalité de Fenchel et ses applications aux lois constitutives[END_REF].

The bipotential in short

First of all, let us recall a basic concept of convex analysis, the subdifferential of a function π in a point x which is the (possibly empty) set:

∂π(x) = {y | ∀x , π(x ) -π(x) ≥ (x -x) : y} . (6)
For more details on convex analysis, the reader is refered for instance to [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Moreau | Fonctionnelles convexes[END_REF][START_REF] Rockafellar | Convex Analysis[END_REF]. Moreover, in mechanics, GSM can represented as a generalized model based on two superpotentials π(x ) and π * (y ), which are depending on a represented strain rate variable x and a stress-like one y. Such a couple of superpotentials satisfies the Fenchel's inequality [START_REF] Fenchel | On conjugate convex functions[END_REF],

∀(x , y ) π(x ) + π * (y ) ≥ x : y (7)
where π(x ) and π * (y ) are convex, lower semicontinuous and conjugate each of the other. The r.h.s. of ( 7) indicates the inner product of x and y . When the equality is achieved, (x, y) is called an extremal couple:

π(x) + π * (y) = x : y
It can be proved that this relation is equivalent to the two following differential inclusions:

y ∈ ∂π(x), x ∈ ∂π * (y).
It is worth to remark that the convexity properties of π and π * are essential in order to state and prove minimum variational principles and use the limit analysis approach2 . When the normality law fails and is replaced by a non-associated flow rule, the classical presentation is based on a yield function (to model the yield criterion) and a plastic potential (to represent the flow rule). Although it is intensively used in the literature, this is, in fact, not very relevant for the variational methods. On the ground of this observation, de Saxcé and collaborators proposed in [START_REF] De Saxcé | New inequality and functional for contact friction: The implicit standard material approach[END_REF][START_REF] De Saxcé | Une généralisation de l'inégalité de Fenchel et ses applications aux lois constitutives[END_REF]) a suitable modeling based on more general generating functions called bipotentials and defined by the following properties:

y ∈ ∂b(•, y)(x) ⇐⇒ x ∈ ∂b(x, •)(y) ⇐⇒ b(x, y) = x : y (9)
In a mechanical point of view, the bipotential represents the plastic dissipation power (by volume unit) and ( 9) is the constitutive law. The couples (x, y) for which ones equivalence (9) holds are called extremal couples. The cornerstone inequality (8) clearly generalizes Fenchel's one (7).

Variational framework of bipotential-based formulations for constitutive laws

Let us now replace the above notations x and y respectively by the strain rate tensor d and the stress tensor σ. In a mechanical point of view, the corresponding bipotential represents the plastic dissipation power (by volume unit) and from Eq.( 9) the constitutive law can be obtained. Accounting for the definition (6) of the subdifferential and the cornerstone inequality (8), the constitutive law reads (Buliga et al., 2009a(Buliga et al., , 2010a;;[START_REF] Laborde | Fixed points strategies for elastostatic frictional contact problems[END_REF]Buliga et al., 2010b):

min d (b(d , σ) -d : σ) = min σ (b(d, σ ) -d : σ ) = 0 . ( 10
)
It is worth remarking that, with respect to the previous minimization problems, the bipotential has the required convexity properties.

Next, let us show how to recover simply the plastic yielding condition F (σ) = 0 by the bipotential formalism. To this end, the first minimization problem in (10) becomes: min

H(d)≤0 (b 0 (d, σ) -d : σ) = 0 , ( 11
)
where b 0 is the finite part of the bipotential when the extremal value is taken. Relaxing the kinematical condition H(d) ≤ 0 by use of Lagrange's multiplier λ, this constrained minimization problem is transformed into an equivalent saddle-point problem

max λ≥0 min d (L(d, σ, λ) = b 0 (d, σ) -d : σ + λH(d)) = 0 , (12)
where L(d, σ, λ) is the lagrangian function. Its stationarity with respect to d:

∂L ∂d = ∂b 0 ∂d (σ) -σ + λ ∂H(d) ∂d = 0
Eliminating the lagrangian multiplier λ in above system of equations, the resultant functional depends only on stress tensor σ. Let us denote it F ; it follows the yield criterion

F (σ) = 0
In a similar way, it is possible to recover the plastic flow rule (3) at a regular point. The second minimization problem in (10) becomes: min

F (σ)≤0 (b 0 (d, σ) -d : σ) = 0 .
Relaxing the plastic yielding condition F (σ) ≤ 0 by use of Lagrange's multiplier λ * , this problem is transformed into an equivalent saddle-point problem max

λ * ≥0 min σ (L * (d, σ, λ * ) = b 0 (d, σ) -d : σ + λ * F (σ)) = 0 , (13)
By calculating the stationarity of the lagrangian L * (d, σ, λ * ) with respect to σ

∂L ∂σ = ∂b 0 ∂σ (d) -d + λ * ∂F (σ) ∂σ = 0
and eliminating λ * , the resultant functional (denoted H) depends only on the strain rate tensor d. Hence, the flow rule can be obtained

H(d) = 0

Case of the non-associated Drucker-Prager materials

The finite value bipotential of non-associated Drucker-Prager model (see Section 2) takes the form [START_REF] Hjiaj | A complete stress update algorithm for the non-associated Drucker-Prager model including treatment of the apex[END_REF]:

b(d, σ) = σ 0 α d m + (β -α) 3σ m -σ 0 α d eq if F (σ) ≤ 0 and H(d) ≤ 0 +∞ otherwise , (14)
In view of what will be done in subsections 5.2 and 5.3, for the homogenization problem, it is convenient to indicate how the derivation of the non-associated yield criterion can be done from ( 12) and ( 14). The lagrangian function reads:

L(d, σ, λ) = σ 0 α d m + (β -α) 3σ m - σ 0 α d eq -(σ e d eq + 3d m σ m ) + λ(βd eq -d m ) .
Its stationnarity with respect to d eq and d m gives:

σ e = (β -α) 3σ m - σ 0 α + βλ , 3σ m = σ 0 α -λ .
Eliminating λ between these relations leads to the plastic criterion :

F (σ) = σ e + 3ασ m -σ 0 = 0 .
Simultaneously, from the second minimization problem of (10), we have the corresponding lagrangian by introducing the multiplier λ *

L * (d, σ, λ * ) = σ 0 α d m + (β -α) 3σ m - σ 0 α d eq -(d : s + 3d m σ m ) + λ * (σ e + 3ασ m -σ 0 ) .
In the same way, the stationnarity with respect to s and σ m reads,

d = λ * 3s 2σ e , ( 15
) (β -α)d eq -d m + αλ * = 0 . (16)
From ( 15) one obtains λ * = d eq . Eliminating λ * in (16) leads to the kinematical condition at the regular points upon the yield surface:

H(d) = βd eq -d m = 0 , (17)
That allows to recover the non-associated yielding rule (3):

d = d + d m 1 = d eq 3s 2σ e + β1 .
For the treatment of the apex, the reader is refered to [START_REF] Hjiaj | A complete stress update algorithm for the non-associated Drucker-Prager model including treatment of the apex[END_REF], which is specifically devoted to this question.

Extended limit analysis of porous materials with a non-associated matrix

Unlike the classical presentation of the non-associated constitutive laws by means of the yield function and the plastic potential, the bipotential formulation naturally opens into a variational formulation; this is crucial for an extension of limit analysis techniques to the context of non-associated laws. We present here the main elements of this variational framework in the context of porous media.

Determination of the macroscopic bifunctional and its variational properties

This presentation is directly done in the framework of homogenization of porous material, considering a reference cell Ω composed of a void ω and a matrix Ω M = Ωω made of an Implicit Standard Material. The macro-cell Ω is enclosed by surface ∂Ω and the void ω by ∂ω. The external boundary of the cell is subjected to a uniform strain rate: v = D • x, x being the position vector at the boundary. The macroscopic stress Σ and strain rate D are then classically defined as volume averages of their microscopic counterpart σ and d:

Σ = 1 | Ω | Ω σ dV , D = 1 | Ω | Ω d dV . ( 18
)
Note that the set of kinematical admissible velocity fields is defined in the following sense:

K a = {v s.t. v(x) = D.x on ∂Ω} . (19)
and the associated strain rate field is given by d(v) = grad s v = 1 2 gradv + grad T v . The set of statically admissible stress fields is:

S a = {σ s.t. div σ = 0 in Ω M , σ • n = 0 on ∂ω, σ = 0 in ω} . ( 20
)
The set of admissible couples is the product A = K a × S a and the set of extremal ones is defined by:

E = {(v, σ) s.t. (d(v), σ) is extremal in Ω M } .
The homogenization problem consists in determining the set A×E of admissible and extremal fields. Owing to the non linearity of the problem, no exact solution can be found in general.

Due to this difficulty, we present an equivalent variational formulation, more appropriate for simple approximations, thanks to relevant choice of trial fields and minimization procedure. By Hill's lemma, any admissible couple (v, σ) complies with:

D : Σ = 1 | Ω | Ω d(v) : σ dV = 1 | Ω | Ω M d(v) : σ dV , (21)
This suggests introducing the following two field macroscopic bifunctional:

B(v , σ ) = 1 | Ω | Ω M b(d(v ), σ ) dV -D : Σ ,
As previously indicated, we are interested for homogenization purpose in finding the admissible and extremal couples (v, σ). In fact, they are solutions of the following simultaneous minimization problems:

B(v, σ) = min v ∈Ka B(v , σ) = min σ ∈Sa B(v, σ ) = 0 . (22)
Indeed, if (v , σ ) is admissible, relation ( 21) and (8) entail:

B(v , σ ) = 1 | Ω | Ω M (b(d(v ), σ ) -d(v ) : σ ) dV ≥ 0 .
In particular, this occurs for admissible couples (v , σ), (v, σ ), (v, σ). Moreover, for the latter, owing to (9):

B(v, σ) = 0 .
In short, one has for all admissible fields v ∈ K a and σ ∈ S a :

B(v , σ) ≥ B(v, σ) = 0 and B(v, σ ) ≥ B(v, σ) = 0 ,
which prove ( 22). Now, let us discuss some relevant aspects of the variational principles for a rigid perfectly plastic matrix such as the one described in the previous sections. The set of plastically admissible velocity and stress fields are respectively defined as:

K p = {v s.t. H(d(v)) ≤ 0 in Ω M } S p = {σ s.t. F (σ) ≤ 0 in Ω M } . ( 23
)
while the sets of licit velocity and stress fields are respectively K l = K a ∩K p and S l = S a ∩S p . We considered the finite valued bifunctional:

B 0 (v , σ ) = 1 | Ω | Ω M b 0 (d(v ), σ ) dV -D : Σ . ( 24
)
the finite valued bipotential b 0 being introduced in (11).

Hence, the bipotential-based variational homogenization problem becomes:

B 0 (v, σ) = min v ∈K l B 0 (v , σ) = min σ ∈S l B 0 (v, σ ) = 0 . ( 25
)
Note that the determination of the above macroscopic bifunctional can be done by means of any of the two minimization principles, providing that the exact stress field or exact velocity field is given.

Application of the variational principle to the plastic porous material

For a rigid perfectly plastic matrix, since b 0 is positively homogeneous of order one in d, there is a trivial kinematical solution to the previous problem (equation ( 25) together with ( 24)) where v and D vanish. The limit analysis approach consists in finding non trivial solutions qualified as ruin mechanisms. It is expected that these non trivial solutions exist only under an equality condition on Σ that can be interpreted as the equation of the macroscopic yielding surface in the model.

It is worth noting that if both D and Σ are chosen arbitrarily, there is in general no solution to the problem (25). In a practical point of view, it is more convenient for instance to fix only Σ and to find D and v satisfying the first minimization problem in (25). Introducing Lagrange's multiplier field x → Λ(x), this constrained minimization problem is transformed into an equivalent saddle-point problem

max Λ≥0 min v∈Ka L(v, σ, Λ) = B 0 (v, σ) + 1 | Ω | Ω M ΛH(d) dV .
We perform a first approximation by imposing Lagrange's multiplier field to be uniform in

Ω M : max Λ≥0 min v∈Ka L(v, σ, Λ) = B 0 (v, σ) + Λ 1 | Ω | Ω M H(d) dV , (26)
that is equivalent to minimize the bifunctional B 0 under the relaxed kinematical condition:

1 | Ω | Ω M H(d) dV = 0 . ( 27
)
Satisfying the kinematical condition only in an average sense but not locally anywhere in Ω M is a strong approximation but leading to easier calculations. As consequence of this approximation, it is crucial to remark that the minimum of B 0 may not be expected to be zero. Nevertheless, in the spirit of Ladevèze's method of the error on the constitutive law [START_REF] Ladevèze | Comparaisons de modèles de milieux continus[END_REF][START_REF] Ladevèze | Accuracy of elastoplastic and dynamic analysis[END_REF][START_REF] Ladevèze | Error estimation and mesh optimization for classical finite element[END_REF][START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF][START_REF] Ladevèze | La maîtrise du calcul en mécanique linéaire et non linéaire[END_REF]Ladevèze et al. , , 2006a,b),b), its value for the minimizer can be used as a variational error estimator [START_REF] Fortin | Modlisation numrique des milieux granulaires par lapproche du bipotentiel[END_REF]. The minimum principle allows obtaining the "better"solution within the framework imposed by the approximations.

Introducing

Y (v, σ, Λ) = 1 | Ω | Ω M b 0 (v, σ)dV + Λ 1 | Ω | Ω M H(d) dV
and considering (24), the Lagrangian function can be recast into

L(v, σ, Λ) = Y (v, σ(Σ), Λ) -D : Σ (28)
from which, as it will be shown by introducing the trial stress and velocity fields, one obtains the macroscopic criterion and flow rule.

For now, let us just indicate that the ultimate step is to solve the Saddle-point problem by computing its subdifferentials with respect to parameters D:

∂L ∂D (Λ, Σ) = 0 (29)
Eliminating the Lagrangian multiplier Λ in the system of functionals ( 29), one obtains

F(Σ(φ, ψ, f )) = 0 (30)
A priori, the above macroscopic criterion depends not only on the porosity f and the friction angle φ, but also on the dilatancy angle ψ of the matrix.

For completeness, the macroscopic non-associated flow rule, with the boundary conditions v = D • x, can be directly obtained from the stationnarity of the lagrangian function ( 28) with respect to the multiplier Λ:

G = 1 | Ω | Ω M H(d) dV = 0 . ( 31
)

The hollow sphere model with a non-associated Drucker-Prager matrix

The major objective of this section is to apply the above bipotential-based variational approach and limit analysis technique to the hollow sphere model, which is made up of a spherical void embedded in a homothetic cell of a rigid-plastic isotropic and homogeneous matrix, the latter being described by a non-associated Drucker-Prager model. The inner and outer radii of the hollow sphere are respectively denoted a and b, giving the void volume fraction f = (a/b) 3 < 1. The hollow sphere is subjected at its exterior boundary to a uniform strain rate tensor D (see Fig. 2). Primarily, we aim at deriving a macroscopic criterion for the non-associated porous material and the corresponding flow rule.

a b r Ω M ω v = D • x

Proposed trial stress and velocity fields

As mentioned in Section 4, in order to derive the macroscopic model, it is indispensable to propose a couple of trial stress and velocity fields. In order to limit the errors due to approximations, we will consider trial fields for which the macroscopic model is exact at least for pure hydrostatic loadings. In Cheng et al. (2012), we obtain closed analytical formula for the limit hydrostatic stresses. For this case, the stress field and limit load do not depend on the dilatancy angle ψ and they are identical to the ones of the associated case with same friction angle φ, previously obtained in [START_REF] Thoré | Closed-form solutions for the hollow sphere model with Coulomb and Drucker-Prager materials under isotropic loadings[END_REF]. Only the collapse mechanism is dilatancy angle dependent. This insensitivity of the hydrostatic limit load to the dilatancy angle agrees with the model of [START_REF] Maghous | Micromechanical approach to the strength propoerties of frictional geomaterials[END_REF] already mentioned. For this reason, and taking into account the symmetry of the hollow sphere model, the trial stress field is considered as the sum of the two following fields:

• A heterogeneous part corresponding to the exact field under pure hydrostatic loadings [START_REF] Thoré | Closed-form solutions for the hollow sphere model with Coulomb and Drucker-Prager materials under isotropic loadings[END_REF]; it reads, in spherical coordinates with orthonormal frame {e r , e φ , e θ }:

σ (1) = A 0 b r 3γ e r ⊗ e r + 1 - 3γ 2 (e θ ⊗ e θ + e φ ⊗ e φ ) (32)
where A 0 is a constant to be determined, s = 1 + 2 α and γ = 1 -2 α 1+2 α , with a loading parameter = ±1, which will be interpreted later.

• A homogeneous part in the cylindrical coordinates with orthonormal frame {e ρ , e φ , e z }:

σ (2) = A 1 (e ρ ⊗ e ρ + e φ ⊗ e φ ) + A 2 e z ⊗ e z ( 33
)
where A 1 and A 2 are also constant parameters.

It should be noted that σ (2) allows to capture the macroscopic shear effect.

The resultant three parameters based trial stress field is defined in the matrix Ω M as:

σ = σ (1) + σ (2) (34)
Note that a vanishing stress field is considered in the void ω.

It is worth to remark that with a stress field is in internal equilibrium, one has:

Σ void = 1 | Ω | ω σdV = 1 | Ω | ∂ω (σn) ⊗ xdS Hence: Σ void m = 1 | Ω | ω σ m dV = (3V ) -1 ∂ω x • (σn)dS
As the continuity condition:

(σn) -+ (σn) + = 0 on ∂ω (35)
is difficult to satisfy by the very simple chosen trial field, we relax it as follows:

(3V ) -1 ∂ω x • ((σn) -+ (σn) + )dS = 0 (36)
which can be equivalently written as:

f -γ A 0 + 2A 1 + A 2 3 = 0 , (37)
On the other hand, the macroscopic stress field is:

Σ = A 0 1 -f 1-γ 1 + (1 -f ) [A 1 (e x ⊗ e x + e y ⊗ e y ) + A 2 e z ⊗ e z ] , (38)
Taking into account (37), the macroscopic mean stress is:

Σ m = A 0 (1 -f -γ ) , (39)
while the macroscopic deviatoric stress reads:

Σ e = (1 -f ) | A 1 -A 2 | . ( 40
)
These two last relations allow to express the stress parameters in terms of the macroscopic stress:

A 0 = Σ m 1 -f -γ , | A 1 -A 2 |= Σ e 1 -f . (41)
Next, following [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF] (see also [START_REF] Thoré | Closed-form solutions for the hollow sphere model with Coulomb and Drucker-Prager materials under isotropic loadings[END_REF])), we adopt, in cylindrical coordinates, the following trial velocity field which depend on the dilatancy angle ψ, not on the friction angle φ,

v = C 0 b r 3/s (ρe ρ + ze z ) + C 1 ρe ρ + C 2 ze z , (42)
with r = ρ 2 + z 2 , s = 1 + 2 β, where is the sign of C 0 . The first term is the exact solution for the hydrostatic case (Cheng et al., 2012). As in Gurson's model [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth -part I: Yield criteria and flow rules for porous ductile media[END_REF] and in its extension to pressure sensitive dilatant materials [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF], this term is completed by two linear terms in order to capture the shear effects. D being the applied macroscopic strain rate, the trial velocity field (42) must comply with the boundary conditions:

v = D • x
In the case of axisymmetric macroscopic strain rate (D = D xx (e x ⊗ e x + e y ⊗ e y ) + D zz e z ⊗ e z ), considered in the present study, C 0 , C 1 and C 2 are such that:

D m = 1 3 trD = C 0 + 1 3 (2C 1 + C 2 ) D zz -D xx = 2 3 (C 1 -C 2 ) (43)
from which it follows:

D e = 2 3 D : D = 2 3 | C 1 -C 2 | (44)
D being the deviatoric part of D.

Closed-form expression of the macroscopic bifunctional

In order to derive the non-associated macroscopic model by solving the saddle point problem (26), we aim now at parametrically expressing the macroscopic bifunctional (24) thanks to the proposed trial stress (34) and trial velocity fields (42). It should be pay attention that the bipotential ( 14) depends on the microscopic mean stress σ m , mean strain rate d m and equivalent strain rate d eq . From ( 34) and ( 42), and considering (41), these quantities can be respectively calculated as

σ m (r) = Σ m 1 -f γ 1 - f γ s b r 3γ , ( 45
)
d m (r) = 1 3 trd = 1 - 1 s C 0 b r 3/s + 1 3 (2C 1 + C 2 ) (46) d eq (r) = 2 3 (C 1 -C 2 ) 2 + (C 1 -C 2 ) 3C 0 s b r 3/s (3 cos 2 θ -1) + 3C 0 s 2 b r 6/s (47)
for which one must have in mind the relations:

d m (r) = D m + C 0 (1 - 1 s )( b r ) 3/s -1 (48) d eq (r) = 2 3 D 2 e + sign(C 1 -C 2 )D e 3C 0 s b r 3/s (3 cos 2 θ -1) + 3C 0 s 2 b r 6/s (49)
For simplicity, let us introduce the normalized stress tensor

T = Σ σ 0 , (50)
and the contribution of the void to average strain rate

D(a) = C 0 f γ 1 + f [C 1 (e x ⊗ e x + e y ⊗ e y ) + C 2 e z ⊗ e z ] , (51)
with γ = 1s-1 . Combining ( 14), ( 22) and ( 50), the normalized macroscopic bifunctional can be obtained as follows:

B0 (v, σ) = B 0 (v, σ) σ 0 = 1 3α tr (D -D(a)) + 1 - β α Π(v, σ) -D : T (52) with Π(v, σ) = 1 | Ω | Ω M d eq dV - 1 | Ω | Ω M 3α σ m σ 0 d eq dV (53)
Hence, it is convenient to introduce

Π(v) = 1 | Ω | Ω M d eq dV ,
which by considering ( 47) and ( 43), can be reduced into

Π(C 0 , D e ) = D e 1 f K(ξ) 1 + τ 2 x -2/s dx , (54)
with

τ = 2C 0 sD e , x = ( r b ) 3 , ξ = 2τ x -1/s 1 + τ 2 x -2/s sign(C 1 -C 2 ), | ξ |≤ 1 .
(55)

K(ξ) = 1 2 π 0 1 + 1 2 (3 cos 2 θ -1)ξ sin θ dθ , (56)
Putting ( 45) into ( 53), Π can be recast into the following parametric form

Π(C 0 , D e , T m ) = 1 - 3αT m 1 -f γ Π(C 0 , D e ) + 3αT m 1 -f γ f γ s I(D e ) (57)
with

I(D e ) = D e 1 f x -γ K(ξ) 1 + τ 2 x -2/s dx . (58)
Finally, the closed form expression of the macroscopic bifunctional reads:

B0 (C 0 , D m , D e , T m , T e ) = 1 α (1 -f ) D m -f γ -f C 0 + 1 - β α Π(C 0 , D e , T m ) -(D e T e + 3D m T m ) (59)
This constitute one of the key practical respect of the study.

Determination of the macroscopic criterion and of the macroscopic flow rule

Having in hand the bifunctional B 0 , we are already now to determine the kinematical admissibility condition (27), written as

βΠ(C 0 , D e ) -(1 -f ) D m -f γ -f C 0 = 0 (60)
which plays the role of the macroscopic flow rule.

Concerning the determination of the macroscopic criterion, let us first introduce a normalized multiplier Λ = Λ/σ 0 , the normalized lagrangian can be written as

L(C 0 , D m , D e , T m , T e , Λ) = L(C 0 , D m , D e , T m , Λ) σ 0 = 1 α -Λ (1 -f ) D m -f γ -f + ΛβΠ(C 0 , D e ) + 1 - β α Π(C 0 , D e , T m ) -(D e T e + 3D m T m ) (61)
For the minimization of the macroscopic bifunctional ( 59

T e (τ ) = ΛβΠ ,De (τ ) + 1 - β α Π,De (τ ) , ( 62
) 3T m (τ ) = 1 α -Λ(τ ) (1 -f ) , (63)
ΛβΠ ,C 0 (τ ) + 1 - β α Π,C 0 (τ ) - 1 α -Λ(τ ) (f γ -f ) = 0 . ( 64
)
from which, we deduce the expression of the normalized multiplier:

Λ(τ ) = 1 α (f γ -f ) + β α -1 Π,C 0 (τ ) f γ -f + βΠ ,C 0 (τ )
.

Eliminating Λ in ( 62) and ( 63), delivers the closed-form macroscopic criterion in the form:

T e = (f γ -f ) β α Π ,De + 1 -β α Π,De + 1 -β α β Π ,C 0 Π,De -Π ,De Π,C 0 f γ -f + βΠ ,C 0 , 3T m = (1 -f ) β α Π ,C 0 + 1 -β α Π,C 0 f γ -f + βΠ ,C 0 . ( 65
)
in which Π is given by ( 57), and it is recalled that γ = 2 β 1+2 β . The explicit expressions for Π ,C 0 , Π ,De , Π,C 0 , Π,De and their antiderivatives Π, Π are calculated and detailed in Appendix A. It is worthy to noted that the above macroscopic criterion (65) is established in a parametric form which depends on the strain rate ratio τ (55) corresponding to the velocity imposed condtion v = D • x. More precisely, the points of plastic limit stress curve can be obtained from the parametric macroscopic criterion (65) for different fixed values of τ . This also allows to deduce the normalized triaxiality T = T m /T e .

Finally, the stationarity of the normalized lagrangian (61) with respect to Λ gives directly the macroscopic flow rule (see ( 60)).

Note also that, macroscopic associated flow rule is obtained by β = α, non-associated one otherwise. Let us recall from (55) that

C 0 = s 2 τ D e (66)
and introduce

P = Π D e (67)
which only depends on τ . Inserting ( 66) and ( 67) into (60) provides a new form of the admissibility condition (60) (flow rule):

D m D e = 1 1 -f βP(τ ) + (f γ -f ) s 2 τ (68)
Hence, the plastic flow direction Υ can be obtained once the value of τ is pre-proposed,

Υ = acot 1 1 -f βP(τ ) + (f γ -f ) s 2 τ (69)
Finally, owing to the matrix compressibility, the void growth rate readily reads,

ḟ = 3(1 -f )D m - 1 | Ω | Ω M trd dV (70)
which, by considering ( 42), ( 46) and ( 66), takes the final expression,

ḟ = 3(f γ -f )C 0 = 3 2 (f γ -f )sτ D e (71)
Note that ( 65), ( 68) and ( 71) are probably some of the most important and practical results of the study. They are the basic blocks of the non-associated constitutive law of the porous material having a non associated Drucker-Prager matrix.

Examination of some special cases

We analyze in this subsection the predictions obtained for some special cases for which results are available in literature. Let us first note that the general problem involves three constants defining the velocity field (42). These constants are linked by the three relations (43), ( 44) and (60). These equations can be easily explicit in the particular cases examined below.

• Hydrostatic case: D e = 0 and C 0 = 0 From ( 43) and ( 47), the microscopic equivalent strain rate in this case reads

d eq (r) = 2 | C 0 | s b r 3/s
It follows from ( 54) and ( 58) that

Π = C 0 β 1 -f γ , I = 2 | C 0 | s 1 -f γ-γ γ -γ in which τ = C 0
De has been also considered. Hence, taking also into account (57), the macroscopic stress is given from (65) in the form:

T e = Σ e σ 0 = 0 , T m = Σ m σ 0 = 1 3α (1 -f γ ) , (72)
which is the plastic limit state for the pure hydrostatic loading Σ e = 0 (traction and compression). This result corresponds to the exact solution in non-associated case derived by Cheng et al. (2012), which have also verified that for the hydrostatically loaded hollow sphere, the limit loads for the non-associated case are the same as for the corresponding associated one given by [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF] (see also [START_REF] Thoré | Closed-form solutions for the hollow sphere model with Coulomb and Drucker-Prager materials under isotropic loadings[END_REF]). Finally, due to D e = 0, (60) readily implies C 0 = D m and then, by ( 71), ḟ = 3(f γf )D m .

• C 0 = 0 and D e = 0 In this case, from ( 47), ( 54) and ( 58) we have

d eq = 2 3 | C 1 -C 2 |= 2 3 D e Π = D e (1 -f ) , I = D e 1 -f 1-γ 1 -γ
The macroscopic limit stresses can be obtained from ( 65)

T m = Σ m σ 0 = 0, T e = Σ e σ 0 = 1 -f (73)
Similarly, the plastic flow direction can be derived from ( 68) or ( 69):

D m D e = β, or Υ = acot β (74)
In this case, the macroscopic admissibility condition appears as the exact couterpart of the microscopic one. For completeness to this case, the void growth rate can be immediately obtained from (71), that is ḟ = 0. All of the above results obtained from C 0 = 0 reflects and proves that the hollow sphere model is under a pure shear loading (73). They provide the same result as the solution obtained by [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth -part I: Yield criteria and flow rules for porous ductile media[END_REF] and [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF]. Obviously, it leads to the conclusion that no matter the matrix of the porous media take a normality rule or not, solution for pure shear loading depends only on the value of porosity, neither on the friction angle nor the dilatancy one. This is a limitation which comes from the simplicity of the trial fields. Equally important, the plastic flow follows a constant and regularly direction (74), which is independent on the friction angle φ. It can be calculated when and only when the dilatancy angle ψ is fixed.

• Case of associated matrix When the matrix complies with an associated flow rule, ψ = φ and β = α in Eqs. ( 1) and (2). It is worthy to indicate that for the pressure-sensitive matrix there is such as β = α = 0. Consequently, one gets from (65) the following macroscopic criterion

T e = f γ -f f γ -f + α ∂Π ∂C 0 ∂Π ∂D e 3T m = 1 -f f γ -f + α ∂Π ∂C 0 ∂Π ∂C 0 (75)
and from ( 69) and ( 71) that the plastic flow direction and the void growth rate, respectively

Υ = acot 1 1 -f αP(τ ) + (f γ -f ) s 2 τ ḟ = 3(f γ -f )C 0 = 3 2 (f γ -f )sτ D e
which is precisely the so called Upper Bound Model (UBM)3 proposed in [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF]. For completeness, we present in Appendix B the derivation of our methodology in the case of associated matrix.

• von Mises matrix For the porous material with an incompressible matrix, for instance the von Mises yield criterion at microscopic level, the pressure-sensitive parameter φ and ψ both vanish, or in another word α = β = 0. In this case, we have s = s = 1, γ = γ = 0 and owing to (60), one obtains D m = C 0 . As a result, (47) then takes the form

d eq = 2 3 D 2 e + 3D m D e (3 cos 2 θ)( b r ) 3 + (3D m ) 2 ( b r ) 6
and (65) reduced to

T e = Σ e σ 0 = ∂Π ∂D e 3T m = 3 Σ m σ 0 = ∂Π ∂D m (76)
Considering ( 54) and ( 71), we get the following macroscopic criterion

T 2 e + 2f cosh( 3 2 T m ) -(1 + f 2 ) = 0 , (77)
and the void growth equation

ḟ = 3(1 -f )D m . ( 78
)
which are the well-known results of [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth -part I: Yield criteria and flow rules for porous ductile media[END_REF] obtained by a kinematical limit analysis approach.

Illustration and numerical validation of the established criterion

In this section, the predictions of established macroscopic criterion (65) in non-associated cases will be firstly compared with the associated one (75) (see also the so called UBM of [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF]) in subsection 7.1. The expected influence of the non-associated feature is clearly illustrated. Next, Finite Element Method (FEM) based limit analysis computations are performed in subsection 7.2 and their results allow to assess the obtained theoretical criterion. For completeness, due to the fact that the plastic flow rule has been formulated in an implicit form except for the pure shear loading, we will perform in subsection 7.3 the illustration of the analytical plastic flow direction (68) for such a particular case, which will be validated from the corresponding FEM solutions.

Preliminary illustration of the established criterion

We aim now at illustrating the macroscopic criterion (65) established in subsection 5.3 both in associated and non-associated cases. As mentioned before, the matrix pressure sensitivity is characterized by the friction angle φ and the dilatancy one ψ for the Drucker-Prager model (see Eqs.( 1) and ( 2)). These two angles must satisfy the condition 0 ≤ ψ ≤ φ < 56 • 18 (see also Eq.( 5)). It is convenient to note that porosity values in geomaterials are relatively bigger comparatively to porous metals or polymers, etc.. Accordingly, in this subsection, we provide illustration of the established criterion for a porosity f = 0.2 and friction angle φ = 30 • . The corresponding associated case ψ = φ = 30 • is denoted AC, while two nonassociated cases are considered; they are respectively defined by dilatancy angles ψ = 15 • (denoted NAC1) and ψ = 5 • (NAC2). As already mentioned, for hydrostatic loadings (traction and compression), the non-associated cases provide the same predictions as that of the associated one (see [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF]; [START_REF] Thoré | Closed-form solutions for the hollow sphere model with Coulomb and Drucker-Prager materials under isotropic loadings[END_REF]. Note again that this observation is in full agreement with the theoretical and numerical results already established in Cheng et al. (2012). Furthermore, unlike the previous work of [START_REF] Maghous | Micromechanical approach to the strength propoerties of frictional geomaterials[END_REF], as mentioned in sections 4 and 5.3, the non-associated cases show in general different yield loci with respect to the associated one: as expected, the yield surface for a non associated case is lower than for the associated one. Note that a decrease of the dilatancy angle leads to a weaker strength, the difference between the cases ψ = 15 • and ψ = 5 • being slight. These results will be investigated in subsection 7.2 by means of numerical results.

To simplify the presentation here, additional results are provided in Appendix C.1 (see Figs. C.10 and C.11); they allow to illustrate the effects of porosity (f = 0.15, f = 0.25) on the macroscopic yield surfaces in the context of the non-associated plastic matrix. In the same appendix, effects of the friction angle φ are also provided for a porosity f = 0.2.

-3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 T m T e AC (ψ=30)°N AC1 (ψ=15°) NAC2 (ψ=5°) -2.4 -2.3 -2.2 0.95 1 f=0.2, φ=30°F
igure 3: Comparison of yield surfaces between the associated case (denoted AC) with dilatancy angle ψ = 30 • and two non-associated cases (65) (denoted NAC1 and NAC2) with ψ = 15 • and 5 • , respectively. Porosity: f = 0.2; friction angle φ = 30 • 7.2. Numerical investigations of the macroscopic yield surface and plastic flow rule in the context of the Druker-Prager non associated matrix In this subsection, the predictions of the established macroscopic criterion will be compared with the Finite Element Method (FEM) solutions. For the FEM analysis, we consider an axisymmetric model of the spherical shell. Hence, owing to the geometrical symmetry, only a quarter of this model is considered by adopting 1500 quadratic axisymmetric elements (see Fig. 4). Moreover, the numerical analysis is carried out in the context of non-associated elastoplasticity and small deformations. The computations are performed by means of ABAQUS/Standard software and a user subroutine MPC (Multi-Points Constraints). The main reason for which we need to enforce MPC conditions in the code is that we have to impose the velocity field v from v = D•x (on the external boundary of the hollow sphere) such that the constraint of constant macroscopic stress triaxiality (T = Σ m /Σ e ) be fulfilled. In practice, as in [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF], this is done by applying a constant macroscopic stress ratio Σ ρ /Σ z corresponding to the desired Σ m /Σ e . Note that the implementation of this procedure has been already described by Cheng and Guo (2007) for their study of voids interaction and coalescence in an associated Drucker-Prager matrix. Coming now to the results, an excellent agreement between the DA and SFA is noted (see Fig. 5), the plastic flow direction (DA) at each FEM point being normal to the yield surface (SFA). Concerning the plastic flow direction (DNA1 and DNA2) for the cases of non-associated matrix, a lack of normality to the corresponding yield surfaces (SFNA1 and SFNA2) is noted. These FEM results proves the non-associated character of the macroscopic flow rule in the case of a non-associated matrix. It must be noted that the lack of normality is more pronounced when the dilatancy angle ψ is small, that is a material with a pronounced non associated matrix.

For completeness, additional results and validations are provided on Figs. 12(a), 13(a), 14(a) and 15(a) in Appendix C.2. This allows to illustrate the effects of the porosity and friction angle.

Validations of the established criterion and of the corresponding flow rule

For validation purpose, the analytical yield surfaces for the associated case (AC) as well as the non-associated ones (NAC1 and NAC2) are compared with the corresponding numerical limit stresses on Fig. 6. The FEM results confirm that the limit stresses of non-associated cases (denoted FNAC1 and FNAC2) and the associated one (denote FAC) are very close in the vicinity of traction dominant region T m > 0 (or Σ m > 0). In contrary, a slight difference is observed in a part of the compression dominate region T m < 0 (or Σ m < 0). The above numerical results validate the predictive capabilities of the analytical criterion.

Let us recall that other validating results showing also the influences of the porosity and friction angle are provided on Fig. 12(b), 13(b), 14(b) and 15(b) in Appendix C.2.

-3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 T m T e AC (ψ=30°) NAC1 (ψ=15°) NAC2 (ψ=5°) FAC (ψ=30°) FNAC1(ψ=15°) FNAC2(ψ=5°) -2.4 -2.2 -2 1 1.1 f=0.2, φ=30°F
igure 6: Comparison between the yield surfaces obtained from the established criteria ( 65) and FEM results (denoted FAC for associated case, FAC1 and FAC2 for non-associated ones).

Concerning the assessment of the obtained macroscopic flow rule (see Eq. ( 68)). We display on Fig. 7 the plastic direction obtained from the analytical function ( 69) and that numerically computed from the FEM computation. Both of these results are illustrated with respect to the different values of macroscopic stress triaxiality T m /T e . The illustrations have been realized for f = 0.2 and φ = 30 • , ψ = 30 • , ψ = 15 • and ψ = 5 • . Noticeable difference between the associated case and the non associated ones is obtained, both for the analytical results and for the numerical data. Moreover, for any fixed value of triaxiality, the obtained value of acot(D m /D e ) representing the plastic flow direction is smaller when the dilatancy angle ψ diminishes. Finally, very good qualitative agreement is observed between the theoretical predictions and the corresponding numerical data. At the difference of porous media with an incompressible plastic matrix [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth -part I: Yield criteria and flow rules for porous ductile media[END_REF], void growth (see Eq. ( 71)) occurs under macroscopic pure shear loading in the case of a pressure-sensitive matrix (as considered in the present work). For this particular shear loading, Eq.( 74) indicates that the plastic flow is only influenced by the dilatancy parameter β, but not by the friction one α. Still for the pure shear, we perform on Fig. 8 a comparison of the plastic flow direction given by the analytical solution (74) and the FEM computations. The following values of fixed material parameters, f = 0.2 and φ = 30 • , are considered. An excellent agreement between the analytic and FEM results is obtained. Finally, Fig. 9 illustrates the variation of porosity with the macroscopic stress triaxiality for a initial porosity f = 0.2, a fixed value of friction angle φ = 30 • and for three dilatancy angles ψ = 30 • , ψ = 15 • and ψ = 5 • . Noticeable difference between the associated case and non-associated ones is observed, particularly for high stress triaxialities.

Conclusion

In this study, a bipotential-based variational framework of ductile porous media has been proposed. It allowed to extend classical limit analysis of porous media to the context of a matrix obeying to a non-associated flow rule. This is generally the case of various porous geomaterials or porous polymers displaying also pressure-sensitivity of the matrix. The proposed variational formulation combines the bipotential theory earlier introduced by de Saxcé et al. (1991) with homogenization techniques. It delivers closed-form expression of the macroscopic criterion for the porous materials, as well as the non-associated flow rule.

As detailed in Section 4, application of the proposed approach to porous materials with non associated plastic matrix (characterized by a friction angle and a dilatancy angle) has led to the formulation of a macroscopic bifunctional. Minimization procedures of the bifunctional have been proposed with respect to a trial kinematically admissible velocity fields (respectively the statically admissible stress fields). This can deliver for the considered class of porous materials an upper bound (respectively lower bound) provided that the exact stress field (respectively the exact velocity field) has been adopted.

Practical implementation of the approach has been done by considering a hollow sphere subjected to uniform strain rate boundary conditions. To this end, a choice of a simple stress field together with a class of trial velocity fields has been made. The whole procedure has allowed to establish a closed form expression of the macroscopic yield function (see equations ( 65)), as well as the macroscopic flow rule (see equations ( 68)), in a parametric form. The non-associated character of the matrix (dilatancy angle) affects not only the macroscopic yield surface, but also the macroscopic flow rule which is shown to be non associated. Moreover, due to the suitable choice of the trial fields, the model also preserves the exact solution established in Cheng et al. (2012) for the hollow sphere with the non associated matrix, under pure hydrostatic loadings. It is also worth noticing that the obtained results allow to retrieve (as a particular case) the kinematically-based model proposed by [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF] for the porous material with an associated Drucker-Prager matrix. This automatically includes the Gurson model for a von Mises matrix [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth -part I: Yield criteria and flow rules for porous ductile media[END_REF]. The predictions of the general model are fully assessed, both for the macroscopic yield surfaces and for the flow rule. To this end, numerous Finite Elements computations in non associated plasticity have been carried out on the hollow sphere. A good agreement has been observed between theoretical results and numerical data for various configurations of porosity, friction angle and dilatancy angle.

Finally, it should be noted that some improvements of the basic model are possible, in particular by searching more refined trial velocity and stress fields. The consideration of voids saturation by an internal pressure will also constitute a challenging extension which will pave the way for various applications in geomechanics including non associated poroplasticity.

The regular solution is classically written in the form of the following power series

2 F 1 (a, b; c; z) = ∞ n=0 (a) n (b) n (c) n z n n! , (A.3)
where (a) n is the Pochhammer symbol defined by

(a) n = Γ(a + n) Γ(a)
= a(a + 1)...(a + n -1)

.

Let us introduce that ι = τ f 1/s Π and I (Eqs. (A.1) and (A.2)) can be indirectly expressed as follows:

P = Π D e = 2 F 1 - 1 2 , - s 2 ; 1 - s 2 ; -τ 2 -f • 2 F 1 - 1 2 , - s 2 ; 1 - s 2 ; -ι 2 (A.4) Q = I D e =s -1 2 F 1 - 1 2 , - ss -1 2 ; 1 -- ss -1 2 ; -τ 2 -f s -1 • 2 F 1 - 1 2 , - ss -1 2 ; 1 -- ss -1 2 ; -ι 2 (A.5)
Finally, considering (A.1), (A.2) and ( 57), the derivatives Π ,C 0 , Π ,De , Π,C 0 and Π,De can be also computed as follows:

Π ,C 0 = τ s 2 -1 2 F 1 1 2 , 1 - s 2 ; 2 - s 2 ; -τ 2 -f 1-2/s • 2 F 1 1 2 , 1 - s 2 ; 2 - s 2 ; -ι 2 (A.6) Π ,De = 2 F 1 1 2 , - s 2 ; 1 - s 2 ; -τ 2 -f • 2 F 1 1 2 , - s 2 ; 1 - s 2 ; -ι 2 , (A.7) Π,C 0 = Π ,C 0 + 3αT m 1 -f γ • f γ s • I ,C 0 -Π ,C 0 (A.8) Π,De = Π ,De + 3αT m 1 -f γ • f γ s • I ,De -Π ,De (A.9)
where I ,C 0 and I ,De have the following expressions

I ,C 0 = 2τ s 1 s -2
Next, by introducing (B.6) and (B.4) into ( 24), the macroscopic bifunctional in the case of an associated matrix can be written as:

B 0 (v, σ) = 1 | Ω | Ω M σ 0 α d m dV -D : Σ . (B.9)
with the relaxed admissibility constraint (B.8): Ω M d m dV = Ω M αd eq dV .

For the first minimization problem of (25), as for the non associated case, it is suggested to introduce the normalized Lagrangian L by omitting σ 0 (or taking it equal to 1), with Λ the Lagrange multiplier which is assumed constant:

L = 1 | Ω | Ωm
L = α ΛΠ + ( 1 α -Λ)[(1 -f )D m -(f γ -f )C 0 ] -(3T m D m + T e D e ) (B.11)
The following minimization relations are obtained:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ L ∂C 0 = α Λ ∂Π ∂C 0 -( 1 α -Λ)(f γ -f ) = 0 ∂ L ∂Dm = ( 1 α -Λ)(1 -f ) -3T m = 0 ∂ L ∂De = α Λ ∂Π ∂De -T e = 0 (B.12)
The first equation of (B.12) delivers the optimal expression of Λ

Λ = 1 α (f γ -f ) f γ -f +α ∂Π ∂C 0 (B.13)
which, when reported in the two last ones, leads to the expression of the macroscopic admissible stress components, the parametric expression of the macroscopic criterion then reads:

3T m = (1 -f ) 1 α - 1 α (f γ -f ) f γ -f + α ∂Π ∂C 0 = (1 -f ) ∂Π ∂C 0 f γ -f + α ∂Π ∂C 0 (B.14) T e = (f γ -f ) ∂Π ∂De f γ -f + α ∂Π ∂C 0 (B.15)
which corresponds to the result obtained by [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF]. It is worth noticing that, due to the relaxation of the admissibility condition, (B.14) and (B.15) do not guarantee the upper bound character of the result.

Appendix C. Complementary results concerning the effects of the porosity, friction angle and dilatancy angle on the macroscopic criterion

In this section, we will illustrate the influences of the porosity f and the friction angle φ, together with that of the dilatancy angle ψ. The analytical results described in Appendix C.1 will be validated from the FEM computations shown in Appendix C.2. andNAC2). Moreover, the same values of material parameters as described in subsection 7.1 (for the case f = 0.2) are respectively adopted for the three cases. As observed through this figure, the yield loci of non-associated cases are lower than the corresponding associated one. It is interesting to point out that the difference between the yield surfaces of the associated case and non-associated one becomes smaller when the porosity is bigger. Next, in order to estimate the influence of the friction angle φ on the proposed macroscopic criterion (65), two groups of comparisons with a fixed value of porosity f = 0.2 but different friction angles φ = 20

Appendix C.2. Assessment of analytical results with FEM computations

The influence of the porosity f on the plastic limit states (directions of plastic flow and limit stresses) are reported on Figs.C.12 and C.13. It is observed that the difference between the non-associated case and the associated one is more significant for f = 0.15 in the compression zone (T m = Σ m /σ 0 < 0) than for the case f = 0.25. This fact numerically confirms that, as in Appendix C.1, the difference of the yield locus between the associated case and the non associated one is smaller with the decrease of the friction angle φ.

Finally, the influence of friction angle φ is assessed through , it is validated that the difference between non-associated case with a fixed dilatancy angle ψ and corresponding associated one is negligible when the friction angle φ is adequately smaller. It is important to point out that for the case φ = 40 • , the FEM yield surfaces are relatively higher than the analytical one; this probable results from the proposed stress fields, for which the internal boundary condition is relaxed because of its homogeneous part.

Figure 1 :

 1 Figure 1: Drucker-Prager model: yield criterion and non-associated flow rule

  (a) b is convex and lower semicontinuous in each argument. (b) For any x and y we have b(x , y ) ≥ x : y(8)(c) For x and y we have the equivalences:

Figure 2 :

 2 Figure 2: Hollow sphere model

  ), one needs to calculate the partial derivatives of the normalized lagrangian (61). Considering τ = τ (C 0 , D e ), Π(C 0 , D e ) = Π(τ, D e ), Π(C 0 , D e , T m ) = Π(τ, D e , T m ) its stationarity with respect to D e , D m and C 0 gives for the normalized stress tensor T = Σ σ 0 :

Figure 4 :

 4 Figure 4: Hollow sphere model: Geometry of the elementary cell and boundary conditions.

Fig. 5

 5 Fig.5displays the FEM results not only for the macroscopic limit stress, but also for the direction of plastic flow. As in subsection 7.1, the values of porosity f = 0.2 and friction angle φ = 30 • are considered here. Also, the direction of plastic flow for the associated case and two non-associated cases are denoted DA(ψ = 30 • ), DNA1 (ψ = 15 • ) and DNA2 (ψ = 5 • ), respectively. Moreover, the numerical yield surfaces, obtained by connecting each FEM point of plastic limit state, are indicated by SFA(ψ = 30 • ), SFNA1(ψ = 15 • ) and SFNA2(ψ = 5 • ), respectively. Note that each FEM point has been obtained by performing computation at fixed stress triaxialities Σ m /Σ e (equivalently at fixed T m /T e ).Coming now to the results, an excellent agreement between the DA and SFA is noted (see Fig.5), the plastic flow direction (DA) at each FEM point being normal to the yield surface (SFA). Concerning the plastic flow direction (DNA1 and DNA2) for the cases of non-associated matrix, a lack of normality to the corresponding yield surfaces (SFNA1 and SFNA2) is noted. These FEM results proves the non-associated character of the macroscopic flow rule in the case of a non-associated matrix. It must be noted that the lack of normality is more pronounced when the dilatancy angle ψ is small, that is a material with a pronounced non associated matrix.For completeness, additional results and validations are provided on Figs. 12(a), 13(a), 14(a) and 15(a) in Appendix C.2. This allows to illustrate the effects of the porosity and friction angle.

  Illustration of FEM results: plastic flow directions (denoted DA for associated case, DNA1 and DNA2 for non-associated ones) and yield surface (denoted SFA for associated case, SFNA1 and SFNA2 for non-associated ones).

Figure 7 :

 7 Figure 7: Analytical result and FEM solutions for plastic flow direction with the fixed values of porosity f = 0.2 and friction angle φ = 30 • .

  Comparison of plastic flow direction for pure shear loading case obtained from Eq.(68) and FEM solutions with the fixed values of porosity f = 0.2 and friction angle φ = 30 • .

Figure 9 :

 9 Figure 9: Analytical result of the plastic void growth with a initial value of porosity f = 0.2. Friction angle φ = 30 • .

  eqd m )dV ] (B.10) Putting (B.6) and (B.7) into (B.10) leads to:

  Appendix C.1. Analytical results Fig.C.10 displays the yield surfaces with a relatively smaller porosity f = 0.15 and a bigger one f = 0.25, both for the associated case (AC) and two non-associated cases (NAC1

  Figure C.10: Comparison of yield surfaces illustrated from (65) between the associated case (denoted AC) with dilatancy angle ψ = 30 • and two non-associated cases (denoted NAC1 and NAC2) with ψ = 15 • and 5 • , respectively.

  Figure C.11: Comparison of plastic flow direction for pure shear loading case obtained from Eq.(69) and FEM solutions.

  Figs.C.14 and C.15 by adopting a fixed value of porosity f = 0.2 and two values of friction angle φ = 20 • and 40 • , respectively. From the illustrations on Figs. 14(b) and 15(b) 5

  of FEM results: plastic flow directions (denoted DA for associated case, DNA1 and DNA2 for non-associated ones) and yield surface (denoted SFA for associated case, SFNA1 and SFNA2 for nonassociated ones).

  between the yield surfaces obtained from the established criteria (65) and FEM results (denoted FAC for associated case, FAC1 and FAC2 for non-associated ones).

Figure

  Figure C.12: Numerical validation for the established model with fixed porosity f = 0.15 and friction angle φ = 30 • .

  between the yield surfaces obtained from the established criteria (65) and FEM results (denoted FAC for associated case, FAC1 and FAC2 for non-associated ones).

Figure C. 13 :

 13 Figure C.13: Numerical validation for the established model with fixed porosity f = 0.25 and friction angle φ = 30 • .

  between the yield surfaces obtained from the established criteria (65) and FEM results (denoted FAC for associated case, FAC1 and FAC2 for non-associated ones).

Figure

  Figure C.14: Numerical validation for the established model with fixed porosity f = 0.2 and friction angle φ = 20 • .

  between the yield surfaces obtained from the established criteria (65) and FEM results (denoted FAC for associated case, FAC1 and FAC2 for non-associated ones).

Figure

  Figure C.15: Numerical validation for the established model with fixed porosity f = 0.2 and friction angle φ = 40 • .

  Figure D.16: Analytical result and FEM solutions for plastic flow direction with the fixed values of porosity f = 0.15 and friction angle φ = 30 • .

See also the use of Shakedown theory by[START_REF] Boulbibane | Application of shakedown theory to soils with non-associated flow rules[END_REF] for non-associated soils.

Lower and upper bound solutions obtained from the Hill's principle and Markov's one.

It should be emphasized that the macroscopic model proposed by[START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF] cannot be seen as an upper bound. For the corresponding demonstration, readers can be referred to[START_REF] Cheng | Homogenization of porous media with plastic matrix and non-associated flow rule by variational methods[END_REF].

For the same reductions of dilatancy angle ψ with respect to the associated one.

For the non-associated cases with material parameters φ = 40 • , ψ = 30 • and 20 • , there are not enough FEM results to represent the plastic flow and to construct the yield surface due to the so strong nonlinearity of the non-associated model that the corresponding FEM plastic limit states cannot be obtained in this paper.

Appendix A. Explicit expressions of Π, Π and the derivatives Π ,C 0 , Π ,De , Π,C 0 , Π,De

In order to explicitly express the closed-form macroscopic criterion (65), we provide as follows the expressions of Π and Π and their partial derivatives with respect to C 0 and D e .

Eqs. (54) for Π(v) and (58) for I(γ) simultaneously contain the term K(ξ), which is smooth over the compactly supported domain with extreme values K max = K(0) = 1 and K min = K(-1) = 0.962. Following [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth -part I: Yield criteria and flow rules for porous ductile media[END_REF], this function is taken to be unity for simplicity of calculation; that reduces (54) and (58) to:

Additionally, Eqs.(A.1), (A.2) and the derivatives Π ,C 0 , Π ,De , Π,C 0 and Π,De can not be calculated into simple forms. Fortunately, they can be expressed by means of the Gauss hypergeometric function (see for example [START_REF] Fine | Basic Hypergeometric Series and Applications[END_REF]) defined by:

This function is a solution of the hypergeometric differential equation

Appendix B. Macroscopic criterion of ductile porous media with an associated Drucker-Prager matrix

We apply here the general procedure proposed in sections 4 and 5 to the particular case of the associated Drucker-Prager matrix (ψ = φ = 0 and β = α = 0). This is the case already studied by [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF] by using the classical kinematical limit analysis approach which is retrieve here from the proposed variational formulation. In this case, the associated flow rule reads:

The volumetric plastic strain is such that:

2)

The plastic flow rule (B.1) is completed at the apex by the condition:

The finite valued bipotential is reduced into:

It should be emphasized that the bipotential in this case takes the same expression as the support function for the associated Drucker-Prager model (see for instance by Salençon (1983)). Moreover, considering the velocity field (42) in the case of associated matrix (see [START_REF] Guo | Continuum modeling of a porous solid with pressure-sensitive dilatant matrix[END_REF])

the microscopic mean strain rate (46) and the equivalent one (47) can be recast into:

and

which obviously do not comply with (B.3) everywhere in the hollow sphere. Due to this difficulty, the idea is to relax the admissibility condition in an average sense by imposing