Limit analysis and Gurson's model

European Journal of Mechanics A/Solids 24 (2005) 800–819 Limit analysis and Gurson’s model Malorie Trillat, Joseph Pastor ∗ Laboratoire Optimisation de la Conception et Ingénierie de l’Environnement (LOCIE), École Supérieure d’Ingénieurs de Chambéry (ESIGEC), Université de Savoie, 73376 Le Bourget du Lac, France Received 28 January 2005; accepted 15 June 2005 Available online 2 August 2005 Abstract The yield criterion of a porous material using Gurson’s model [Gurson, A.L., 1977. Continuum theory of ductile rupture by void nucleation and growth – Part I: Yield criteria and flow rules for porous ductile media. ASME J. Engrg. Mater. Technol. 99, 2–15] is investigated herein. Both methods of Limit Analysis are applied using linear and conic programming codes for solving resulting non-linear optimization problems. First, the results obtained for a porous media with cylindrical cavities [Francescato, P., Pastor, J., Riveill-Reydet, B., 2004. Ductile failure of cylindrically porous materials. Part 1: Plane stress problem and experimental results. Eur. J. Mech. A Solids 23, 181–190; Pastor, J., Francescato, P., Trillat, M., Loute, E., Rousselier, G., 2004. Ductile failure of cylindrically porous materials. Part 2: Other cases of symmetry. Eur. J. Mech. A Solids 23, 191–201] are summarized, showing that the Gurson expression is too restrictive in this case. Then the hollow sphere problem is investigated, in the axisymmetrical and in the three-dimensional (3D) cases. A plane mesh of discontinuous triangular elements is used to model the hollow sphere as RVE in the axisymmetrical example. This first model does not provide a very precise yield criterion. Then a full 3D model is applied (using discontinuous tetrahedral elements), thus solving nearly exactly the general three-dimensional problem. Several examples of loadings are investigated in order to test the final criterion in a variety of situations. As a result, the Gurson approach is slightly improved and, for the first time, it is validated by our rigorous static and kinematic approaches.  2005 Elsevier SAS. All rights reserved. Keywords: Porous media; Gurson’s model; Limit analysis; Lower and upper approach; Linear and conic programming 1. Introduction Concerning the ductile failure of porous materials, Gurson’s criterion (Gurson, 1977) is the most widely accepted because it is based on a homogenization method and on the kinematic approach of limit analysis. The plastic domain is approached from the outside by a semi-analytical approach, proved to be an upper bound approach by Leblond (2003). Gurson’s model treats a hollow sphere with macroscopic strain imposed on the boundary. The criterion that he proposed for a rigid plastic isotropic matrix around a spherical cavity is expressed as follows:
√  2 Σeqv 3 Σm + 2f cosh (1) =1+f2 2k 3k 2 * Corresponding author. E-mail address: Joseph.Pastor@univ-savoie.fr (J. Pastor). 0997-7538/$ – see front matter  2005 Elsevier SAS. All rights reserved. doi:10.1016/j.euromechsol.2005.06.003 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 801 where f is the porosity rate of the material studied and k the flow stress in shear or cohesion. Σeqv is the macroscopic equivalent stress and Σm the macroscopic mean stress. By definition, this model does not take into account interactions between cavities. The original Gurson criterion for porous materials has undergone various modifications to improve its adequacy with experimental and numerical results. Tvergaard (1981) uses empirical parameters to take into account cavity growth and the cavities coalescence. The Richmond and Smelser (1985) criterion introduces a parameter that can take into account shear bands, which could dominate the plastification process. The Sun and Wang (1989) criterion is the only one to be based on an inner approach of the Gurson model under a macroscopic stress on the hollow sphere boundary. The model of Gologanu et al. (1997) extended the Gurson model to the spheroidal void incorporating void shape effects. The Rousselier (2001) criterion is based on a macroscopic thermodynamic approach of nonreversible phenomena under the local equilibrium assumption, i.e. the state of a system can be defined locally using the same variables as at the equilibrium state. The parameter values of this criterion must be determined experimentally. Numerous models have been developed for metal processing, such as Kim and Kwon (1992). More recently, Bigoni and Piccolroaz (2003) has proposed a macroscopic yield function to agree with a variety of experimental data relative to soil, concrete, rock, metallic and composite powders, metallic foams, porous materials and polymers. Their yield function is presented as a generalization of several classical criteria, also giving an approximation of the Gurson criterion. In Francescato et al. (2004) and Pastor et al. (2004), both approaches of limit analysis make it possible to determine the yield criteria of a porous material with cylindrical cavities. In this case of porous material, we have shown that a Gurson-like formulation cannot properly represent the general solution. A realistic formulation must be specified in terms of at least three loading parameters, the Gurson expression being too restrictive. In the present paper, the same tools of homogenization, limit analysis and optimization are used to obtain the plastic domain in the case of a porous material with spherical cavities. The Gurson model idealizes the porous material as a single cavity in a homothetic cell composed of a rigid plastic Mises material, called the Representative Volume Element (RVE) in the following. The homogenization method links macroscopic stress and strain rates Σij and Eij with microscopic stress and strain rates σij and vij by the averaging relations: Σij = 1 V  V σij dV , Eij = 1 V  vij dV . (2) V This RVE is submitted either to a an average strain rate on the boundary of RVE (ui = Eij xj , uniform strain rate loading) or to an average stress (σij nj = Σij nj , uniform stress loading ) leading to different mechanical problems. Here, we essentially analyse the first case, i.e. ui = Eij xj on the RVE boundary, for testing the original Gurson criterion. Then we use the limit analysis static (or lower bound), and kinematic (or upper bound), approaches. They lead to non-linear optimization problems (because of the plasticity criterion), second order conic programing problems (SOCP) in fact. They are solved either by XA (after a pre-linearization of the non-linearities) or directly by MOSEK (an SOCP code), both optimization codes based on so-called Interior Point methods. These optimizers are the result of intensive algorithmic research since the pioneering work of Karmarkar in 1984. XA is a linear programming code: a pre-linearization is done using the very efficient Ben-Tal and Nemirovski (2001) method detailed in Pastor et al. (2004). The drawback of this linearization is the amount of memory required for the resulting large-scale problems. On the contrary, MOSEK solves the problem directly and therefore can deal with thinner RVE mesh. After recalling the main results for porous materials with cylindrical cavities in Section 2, the admissible loading domain of the Gurson problem is investigated using the hollow sphere as RVE. In Section 3, the problem is considered an axisymmetric one (subject to an axisymmetric loading), using a plane mesh (and triangular elements) to discretize the RVE. This first model does not provide a sufficiently accurate yield criterion. Therefore, in Section 4, a three-dimensional mesh (and tetrahedral elements) is used. These models give a quasi-exact yield criterion, improving the Gurson criterion slightly. They do not confirm the corner found in Francescato et al. (2004) and Pastor et al. (2004) on the average stress axis in a porous material with cylindrical cavities. More precisely, the kinematic 3D result is very close to the kinematic axisymmetric bound, which slightly improves the Gurson approach. This proves that the kinematic axisymmetric approach is excellent. The static 3D approach improves the static axisymmetric approach and becomes very close to the kinematic approach. Then the Gurson criterion, for the first time, is validated as a good approximation with both our kinematic and static approaches. Moreover, the kinematic axisymmetric approach is reinforced as a good prediction tool of the macroscopic isotropic criterion. Furthermore, we confirm that, at least in the case studied, Gurson’s criterion does not depend on the third invariant of the stress tensor. Finally, using the 3D model, macroscopic plane strain and plane stress loadings are also investigated. 802 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 Fig. 1. Yield criteria in plane stress (left), plane strain (right) and Gurson’s results. Porosity f = 0.16, 672 element-mesh. 2. Previous results for porous media with cylindrical cavities In a recent paper (Francescato et al., 2004; Pastor et al., 2004) we tested the validity of the Gurson criterion for a porous material with cylindrical cavities. The yield criteria were determined within the framework of limit analysis in plane strain, plane stress, generalized plane strain and 3D cases. The main results are recalled and complemented in this section, without the description of the numerical methods. 2.1. Plane stress and plane strain cases To study a porous material with cylindrical cavities, the RVE is a hollow truncated cylinder (outer radius of 5 length units, cohesion k of 1 stress unit) with specific boundary conditions; it is made up of a homogeneous Mises material whose radius is such that the cavity rate equals the porosity rate of the real material. Both static and kinematic methods of limit analysis lead to a non-linear optimization problem that is solved using the Ben-Tal and Nemirovski (2001) linearization algorithm. In plane stress, the macroscopic stresses are defined by: √ 3 1 Σh = (Σx + Σy ), (Σx − Σy ) Σcp = 2 2  2 + Σ 2 ), which is why Σ and the macroscopic equivalent stress by Σeqv = (Σcp cp is plot versus Σh . In plane strain, the h macroscopic stresses are defined by the same relationships, using Σps instead of Σcp . In addition, Σeqv = |Σps | = Σps (being Σps non-negative in the present case), then Σeqv is a true loading parameter. The curves on Fig. 1 show the numerical results in plane stress (left), plane strain (right) and Gurson’s results. The numerical static and kinematic yield curves are not matched by any previous analytical yield criteria, particularly in the vicinity of the Σh axis, except the Rousselier (2001) criterion. Pellegrini (2002) gives a heuristic expression to fit our numerical curves in plane stress for the porosity rates of 16, 10, 5, 1%: 
sinh (Σps /k) Σ2 sinh (Σm /2k) 2 + 2m + =1 (3) sinh (Yps /k) sinh (Ym /2k) Yeqv with: √ 3k(1 − f ) exp(−af 2/3 ), a = 0, 501, √ Ym = k ln f, Yeqv = 3k(1 − f ), Yps = In plane strain, we recently found that this expression also fits our numerical criterion by taking a = 1/3 and Yeqv = 10k. It is represented on Fig. 2 for the porosity rates of 16 and 10 % (1152 elements), 5 and 1 % (3600 elements). 2.2. Generalized plane strain case In the case of generalized plane strain, the macroscopic stresses are defined by: √ 3 1 1 (Σx − Σy ), Σps = Σgps = (Σx + Σy ) − Σz . Σh = (Σx + Σy ), 2 2 2 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 803 Fig. 2. Modified Pellegrini’s heuristic expression and our numerical results in plane strain. Fig. 3. Yield criterion in generalized plane strain: Σeqv = f (Σh ) numerical curve and Gurson’s results. Porosity f = 0.16, 672-element mesh. Fig. 4. Generalized plane strain: Σps = f (Σh ), Σgps fixed. Porosity f = 0.16, 672-element mesh. For clarity, only the static results are plotted on Figs. 3–5, keeping in mind that static and kinematic results are very close. Fig. 3 shows that the yield criterion in the (Σeqv , Σh ) plane verifies nearly exactly Gurson’s values. But in the (Σps , Σh ) plane in function of Σgps , the corner appearing for Σgps = 0 (plane strain) vanishes as Σgps increases (see Fig. 4). This fact is 804 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 Fig. 5. Yield criterion in generalized plane strain case in 3D-axis representation: (Σps , Σh , Σgps ). Porosity f = 0.16, 672-element mesh. Fig. 6. Discontinuous mesh with 672 elements and the RVE. certainly related to the anisotropy (i.e. transverse isotropy here) of the homogenized material. It demonstrates that a Gurson-like expression of the criterion is too restrictive. Fig. 5 is a 3D representation in the (Σps , Σh , Σgps ) axis of the yield criterion in the generalized plane strain case. It can be seen that the criterion must be expressed using at least three loading parameters. To investigate the influence of the average shear strain, the mechanical problem must be studied from the 3D point of view in the (x, y) plane. By meshing a half-ring with the appropriate symmetry conditions, we obtained the criterion under shear loading for the porous material with cylindrical cavities. As a result, the Gurson criterion overestimates the true criterion: the influence of the axial shear stress needs specific terms in the macroscopic criterion. To conclude for porous materials with cylindrical cavities, the Gurson criterion appears to be insufficient, even erroneous in some cases. We have proved these facts using both upper and lower approaches of limit analysis, giving very precise determination of the true criterion. In the generalized plane strain case, an analytical expression of this true criterion cannot be a Gurson-like criterion, but it must take the form of a function of the loading parameters, at least in a three-dimensional representation. In the following, the case of porous media with spherical cavities is investigated. 3. Axisymmetric approach 3.1. Description of the problem From the axisymmetry assumption, the problem only needs a plane finite element mesh. To model the RVE, from the isotropy of the homogenized material, a quarter of a ring (outer radius of 5 length units, cohesion k of 1 stress unit in all the following) is discretized in triangular elements in the reference frame (ORθZ). This frame rotates around the Z axis and a half hollow sphere is thus obtained (Fig. 6). Indeed, the loading is also axisymmetric. Under the boundary conditions ui = Eij Xj , the total dissipated power rate P tot can be written as follows with the actual loading parameters: P tot /V tot = Σm Em + Σps Eps + Σgps Egps (4) where the macroscopic stresses and the macroscopic strain velocities are defined by: √ 3 ΣX + ΣY 1 − ΣZ , Σps = (ΣX − ΣY ) Σm = (ΣX + ΣY + ΣZ ), Σgps = 3 2 2
 2 EX + EY 1 − EZ , Eps = √ (EX − EY ). Em = EX + EY + EZ , Egps = 3 2 3 From axisymmetry, EY Z = EZX = 0 and EX = EY (Fig. 6); from the isotropy of the homogenized material, EXY = 0 is assumed, without loss of generality. Then the following relationships hold: ΣY Z = ΣZX = ΣXY = 0 and ΣX = ΣY , i.e. Σps = 0. M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 805 In the present case, the total dissipated power rate is reduced to the expression: P tot /V tot = Σm Em + Σgps Egps . (5) The Mises macroscopic equivalent stress Σeqv is linked to the stresses Σgps and Σps by the relation: 2 = Σ2 + Σ2 = Σ2 . Σeqv gps ps gps (6) The macroscopic strain rate is imposed via the following boundary conditions in the reference frame (OXY Z): uX = EX X, uY = EX Y , uZ = EZ Z. In the frame (ORθZ), these boundary conditions become: uR = EX R and uZ = EZ Z. 3.2. Static method Owing to the isotropy of the Mises material, in the reference frame (ORθZ) the Mises criterion is written: f(σ ) =

2 2 σR + σZ + (σR − σZ )2 + (2τRZ )2 − (2k)2 . − σθ √ 2 3 (7) The need for a linear variation field to verify both the criterion and the equilibrium conditions throughout the element leads to the following stress field in each triangular element (Pastor and Turgeman, 1982):  σR = aσ + bσ Z + cσ R,       σθ = aσ + bσ Z + 2cσ R, (8)  σZ = dσ + eσ Z + fσ R,      τRZ = −eσ R/2 with aσ , bσ , cσ , dσ , eσ , fσ reals. The union of these real constants of all elements form the variables of the final optimization problem. The following conditions must be imposed to make the whole stress field admissible: • Definitions of the macroscopic stresses ΣX , ΣY , ΣZ and loading parameters Σm , Σgps . • Continuity conditions: the stress vector components σn and τnt are continuous across any common side, of normal n, of the elements. • Boundary conditions: the stress vector is equal to zero on the cavity boundary. • Symmetry conditions: the shear stress τRZ is zero on the vertical and horizontal sides (Fig. 6). • σ is plastically admissible i.e. f ( σ )
0. The Mises criterion is linearized to solve the problem with the linear pro- 2 2 2 2 gramming √ code XA. Indeed, the yield criterion (7) in each element apex can be written as: x + y + z
α with x = ((2/ 3 )((σR + σZ /2) − σθ )), y = (σR − σZ ), z = (2τRZ ) and α  0. It is linearized using the Ben Tal and Nemirovski algorithm twice (see Pastor et al., 2004, for a presentation of this algorithm, called the BTN algorithm in the following). Finally Σgps (= Σeqv here) is maximized for fixed positive values of Σm . After the optimization, a post-process checks that the microscopic stress field is statically and plastically admissible. But as shown in Fig. 8, the discontinuous stress field defined in (8) does not supply a satisfactory result. For example, for Σeqv = 0 the static criterion is remote from the Gurson criterion. To improve this inner approach for high triaxilities, a global continuous field, derived from the solution of the hollow sphere under inner pressure given in Chakrabarty (1987), is added to the locally defined stress field. This continuous field is the solution of the problem of a spherical cavity (inner radius of a, outer radius of b) subject to an outer isotropic stress field p. This field is written in the spherical frame (Fig. 7):  √ b  σr = −2σ0 ln + p; σ0 = k 3,    r     σθ = σ0 + σr , (9)  σφ = σθ ,        σ (a) = 0 ↔ p = 2σ ln b . r 0 a 806 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 Fig. 7. Spherical frame. In the (ORθZ) frame of the study, we obtain:  σR = σr + σ0 sin2 θ,       σZ = σr + σ0 cos2 θ, (10)   σRZ = −σ0 cos θ sin θ,     σθ = σr + σ0 . From the spherical symmetry of the continuous problem, we can verify:  √ b tr[σ ]XY Z  tr[σ ]rθϕ  1 σX dV = ΣX  = = = 2 3k ln = 3.1741k. V 3 3 a (11) V After optimization, a post-correction is needed because the addition of a continuous field introduces a small amount of nonlinearity on each triangular element. As can be seen on Fig. 8, this refinement is pertinent: the addition of the continuous stress field gives results that are closer to the kinematic curve. For example, on the horizontal axis the static value now equals the analytical and exact value Σm = 3.1741k, as expected. But in spite of the addition of this analytical stress field, the static results still are not good enough, specially for low Σm values. 3.3. Kinematical method To obtain linear strain rates, the components uR and uZ of the displacement velocity field are defined as: uR = R(au + bu R + cu Z), uZ = du + eu R + fu Z + hu RZ + iu R 2 + ju Z 2 (12) with au , bu , cu , du , eu , fu , hu , iu , ju reals (Pastor et al., 2002). From the axisymmetry, uθ = 0. Indeed, this velocity displacement field only takes into account the RVE shape change. As in Gurson (1977), a global incompressible radial field in gu /r 2 (gu real) is added to provide a global volume change. The displacement velocity field is quadratic, and, as in the static program, any side of the triangle is a potential discontinuity line. The following conditions must be imposed to make the displacement velocity field admissible: • From the linearity of the strain rate v, the condition tr v = 0 is imposed at each element apex. • Across the discontinuity lines: the jump of the normal component [un ] of the displacement velocity is null across a discontinuity line. The jump of the tangential component [ut ] equals the difference of the two plastic multipliers: [ut ] = ξ1 − ξ2 with ξ1 , ξ2  0. As the displacement field is quadratic, this equation is written at the ends and the middle of each element side. • Symmetry conditions: uZ = 0 on the horizontal boundary. As above, this equation is written at the middle of each element side and at its two ends. • Uniform strain boundary conditions: uR = EX R and uZ = EZ Z are imposed on the RVE boundary. This condition is written at four equidistant points on the circle of the same extremities as the triangle side along the boundary. • Definition of the macroscopic strain rates Em , Egps . • Definition of the functional. The total power rate P tot is the sum of the dissipated power rate P disc along the discontinuity lines and the volume power rate P vol : P tot = P vol + P disc (13) M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 807 with P disc =  k|[ut ]| dS (14) Sd where [ut ] is the tangential velocity jump across a discontinuity and Sd is the set of discontinuity surfaces in the mesh. By assuming a linear variation of the plastic multipliers ξ1 , ξ2 along the half-side, the integration of k(ξ1 + ξ2 ) gives an upper bound Pξ of the dissipated power rate. The volume power rate P vol is defined by:  (15) P vol = Π( v ) dV V with Π( v ) is the unit power rate here defined by: Π( v ) = √  2 + v 2 + v 2 + 2v 2 . 2k vR θ Z RZ (16) Using the incompressibility condition, Π( v ) becomes: Π( v ) = 2k
√ 2
2 1 3 2 . (vR + vZ ) + (vR − vZ ) + vRZ 2 2 This expression is linearized using twice the BTN algorithm at each element apex, and is written as: √  3 1 (vR + vZ ), y = (vR − vZ ), z = vRZ , x 2 + y 2 + z2
β, x = 2 2 (17) (18) with two BTN cones and two additional variables, β being the second one. The convexity of Π( v ) is used to upper-bound P vol by integrating 2kβ, β taken as linearly varying upon the triangle:  (19) P vol
Pβ = 2kβ dV V Finally, the optimal value of Σgps is obtained as: opt Σgps = Min 0E (Pβ + Pξ )/V tot − Σm m 0 Egps (20) 0 and Σ 0 are fixed. where Egps m After optimization, using the solution velocity field, all imposed conditions are verified in order to check the admissible character of the solution. The total power rate is calculated in adding the volume and discontinuity power rates. Since the dissipated power rate is non-linear, the following numerical integration method is used to compute the volume power rate, by dividing the nt triangular element in nsub sub-triangles. According to (18), the volume power rate is calculated by: nt nsub P vol = 2π 2k j =1 i=1  2 + y 2 + z2 Rij Sij xGij Gij Gij (21) where Gij is the position of the centre of gravity (of horizontal coordinate Rij ) and Sij the area of each triangular subdivision. The power rate dissipated on the discontinuities is calculated in an analogous manner by dividing the discontinuity segment into a sufficient number of sub-segments. The admissible character of the displacement velocity field is also checked. As expected, the post-processing of the total power rate becomes almost useless with the linearization of the criterion used, because the correction is about 0.1% of the values of Σm and Σgps . For a 672-element problem and k = 8 for the BTN coefficient (i.e. 28 linearization planes for each cone), the problem has 108,849 rows, 80,489 columns and 459,999 non-zero terms. The optimization is done in 30-90 s with the XA code on a Apple G5-2GHz using 2G of RAM. 808 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 Fig. 8. Yield criterion of a spherical cavity material: axisymmetric method, porosity 6.4%, 672-triangle mesh. 3.4. Results Fig. 8 shows the macroscopic equivalent stress Σeqv versus the macroscopic mean stress Σm for a 672-element mesh drawn on Fig. 6 for a porosity rate of 6.4%. We can see that the kinematic curve is slightly exterior, close to the Gurson criterion, but the static curve is distant for low triaxialities, not close enough to qualify the obtained yield criterion accurately. Leblond (2003) points out that Gurson’s limit analysis involves an approximation resulting in an overbounding the plastic dissipation. Indeed, our rigorous kinematic analysis is interior to the Gurson yield surface. This is an indication that the mesh used contained enough elements for the calculation being accurate. Fig. 9 shows the evolution with Σm of the displacement velocity field. For Σm  2.5k, the field becomes higher on the cavity boundary and the volume cavity expansion offsets the outer volume expansion. To go further we need to study the RVE with three-dimensional geometry, displacement velocity and stress fields. 4. 3D approach 4.1. Mesh From the isotropy of the homogenized material, without generality loss, the loading E is taken as principal, i.e. Eij = 0 for i = j . Therefore an eighth of the hollow sphere is meshed by using triangular basis prisms as shown on Fig. 10. On this figure, each triangle on the outer surface is the top basis of a prism (and so on going to the cavity). Fig. 11 shows how each prism is divided into two tetrahedrons and three pyramids. Each pyramid is also divided into four tetrahedrons. Hence, each prism is meshed using 14 tetrahedral elements. Each triangular surface between two tetrahedrons is a potential surface of discontinuity. This discontinuity surface is characterized on Fig. 12 by its normal n of components (nx , ny , nz ). Fig. 13 shows the outer surface of the 14,000-element mesh used. Here, the outer surface and the homothetic surfaces are made up of 100 triangles. 4.2. Kinematical method The displacement velocity field is chosen linear in x, y, z in each tetrahedral element, and any triangular surface common to two tetrahedrons is a potential surface of velocity discontinuity. The following conditions must be checked so that a displacement solution field is admissible: • Incompressibility condition: tr v = 0 in each element. • Symmetry conditions: the normal velocity is zero on each coordinate plane. M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 809 Fig. 9. Displacement velocity field versus Σm . • Boundary conditions: uij = Eij xj (with Eij = 0 for i = j ), so the macroscopic strain rates Eij are the average of the microscopic ones vij . • Continuity of the normal displacement, i.e. [un ] = 0 at each apex of a triangular discontinuity surface. • Definition of the functional. 810 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 (a) (b) Fig. 10. General view and meridian plane of a 672-3D-element mesh. (a) General view of an eighth of tetrahedron meshed hollow sphere. (b) Mesh view in y = 0 plane. Fig. 11. Decomposition of the prism into two tetrahedrons and three pyramids. Then each pyramid is cut into four tetrahedrons. Fig. 12. A discontinuity surface characterized by the pairs (l1 , l2 ), (l3 , l4 ) and (l5 , l6 ) and by its normal n. Fig. 13. 14,000-element mesh, outer view. In the present case, the total dissipated power rate is expressed as: P tot /V tot = Σm Em + Σgps Egps + Σps Eps . Moreover, P tot is defined in (13) and the volume power rate is defined in (15) with: (22) M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 √  2 + 2v 2 + 2v 2 Π( v ) = k 2 vx2 + vy2 + vz2 + 2vyz xz xy
√ 2
2 3 1 2 + v2 + v2 . (vx + vy ) + (vx − vy ) + vyz = 2k xz xy 2 2 The power rate dissipated on the discontinuities is given by:    P disc = k|[ut ]| dS = k [ut1 ]2 + [ut2 ]2 dS Sd 811 (23) (24) Sd Sd is the set of discontinuity surfaces. [ut1 ] and [ut2 ] are the tangential displacement velocities in the discontinuity side interelements. In order to solve this problem, we use the second order conic programming optimization code MOSEK (2002). Conic optimization is a generalization of linear optimization, which allows the formulation of a class of non-linear convex optimization problems. The main idea of conic optimization is to include constraints of the form x t ∈ C in the optimization problem where x t consist of a subset of the variables and C is a convex cone. Must be remember that C is a convex cone if and only if x ∈ C ⇒ αx ∈ C for all α  0. MOSEK can handle conic constraints (called a quadratic cone) of the type:  n C = x ∈ Rn+1 :  xj2
xn+1 (25) j =1 Compared to XA, the linear optimization code used to solve the axisymmetric problem, this code solves larger problems. The expressions (23) of Π( v ) and (24) of P disc are directly written in the form (25) for each tetrahedron (Π( v ) is constant in the tetrahedron) and for each apex of the discontinuity triangle, respectively, using additional variables as in axisymmetry. After bounding P disc over the discontinuity triangles in the same manner than in axisymmetry, the final objective results in a linear expression, as in (20), adequate fixed components of the macroscopic stress being given by the static approach. After optimization, the admissible character of the solution field is checked. To take into account the non-linearity of the discontinuity power rate, each discontinuity surface is divided in 64 sub-triangles and a numerical integration is performed: the corresponding correction on the total power rate results is very small. The incompressibility condition and the maximum jump across a discontinuity surface is about 10−7 for k = 1 (stress unit) and the outer radius equals 5.0 (lenght units), i.e. is highly precise when using an interior point method. For a 1750-element mesh, MOSEK is much quicker than using XA with the BTN linearization. For the kinematic problem with 14,000 tetrahedrons (334,204 rows, 497,706 columns, 2,497,681 non-zero terms), MOSEK only requires roughly 2.5 hours on a Transtec PC 2.2GHz using 1.5G of RAM. Keeping in mind that the axisymmetric value is also a 3D kinematic value, the tests did not decrease these upper bounds (see Fig. 14). This result confirms the pertinence of our above kinematic axisymmetric. Consequently, we have to implement a true 3D-static method to improve the axisymmetric static results. Fig. 14. Yield criterion. Comparison with the axisymmetric method. Porosity 6,4%, 14,000-tetrahedron mesh. 812 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 4.3. Static method As in the kinematic method, the microscopic stress field is chosen as linearly varying in x, y, z in each tetrahedral element and it can be discontinuous through any element boundary, as in Pastor (1983) and Pastor et al. (1990). As above, we want to find the equivalent macroscopic stress Σeqv defined by: 2 = Σ2 + Σ2 Σeqv gps ps (26) In the 3D-general case, the Mises criterion is written: 2

2 σx + σx − σz + (σx − σy )2 + (2τyz )2 + (2τxz )2 + (2τxy )2
(2k)2 f(σ ) = √ 2 3 (27) To have a microscopic stress solution field statically and plastically admissible, the imposed conditions are implemented, briefly: • In each element, the equilibrium equations σij,j = 0 in the Cartesian frame. • Continuity conditions: the stress vector components are continuous across a discontinuity surface of normal n: for each discontinuity surface, the continuity of the stress vector Ti = σij nj is written between the three opposite pairs (see Fig. 12) of the apices of the discontinuity triangle. • Boundary conditions: the stress vector Ti = σij nj is null at each apex of the element sides of the boundary cavity. • Definitions of macroscopic stresses Σij as the average of the microscopic stresses σij . In the case of average strain boundary conditions, the Σij are free variables. In the case of average stress boundary conditions, the Σij are fixed: at the three apices of the element sides of the boundary we write: σij nj = Σij nj . • Symmetry conditions: the microscopic tangential stresses are null on the three coordinate planes. • Stress field plastically admissible: the criterion (27) is linearized from the inside, using four BTN cones as in Pastor et al. (2004) to use XA, or directly written in the form of (25) to solve the final problem with MOSEK. We wish to determine Σeqv , the macroscopic equivalent stress defined by (26). In the case of an axisymmetric loading, i.e. Σx = Σy , Σps is zero and Σgps is optimized, noting that the result is also a static value for the general 3D case. As in the static axisymmetric method, a global continuous field was added to the locally defined stress field to improve the lower bound approach for high triaxialities. This extra field slightly improves the static curve, much less than in the previous axisymmetric situation. Hence this field is not used to simplify the post-processing. As usual, the admissible character of the solution stress field is verified a posteriori. After reading the macroscopic and microscopic stresses, (27) is calculated in each element apex and the functional optimal value is corrected, in the case of a BTN pre-linearization. All other conditions are also verified. 4.4. Results under axisymmetric loading On Fig. 14, the static approach is very close to the kinematic approach obtained with the axisymmetric method. The static axisymmetric curve is far below the static 3D curve for low Σm /k. Finally, using the full 3D model, we have seen that the axisymmetric kinematic approach is excellent and could be used as the macroscopic criterion. The static approach confirms the kinematic approach and shows that the Gurson’s criterion is a good approximation for spherical cavity media. The Gurson criterion was improved at low triaxialities by Garajeu (1995). Its estimate is also an upper bound for the plastic dissipation which for f = 6.4% and Σm /k = 0 gives Σeqv /k = 1.5995 compare to 1.6212 for Gurson. This value is in good agreement with our numerical kinematic result (1.5911). The next table compares XA and MOSEK, the two optimizers, for a 1750-tetrahedron mesh solved on a Transtec PC 2.2GHz with 1.5G RAM. Optimizer Number of rows Number of columns Number of terms Number of iterations Time (s) XA (linear) MOSEK (conic) 485,210 79,210 343,010 84,011 1,514,606 418,706 ≈ 45 ≈9 ≈ 1500 ≈ 70 In general, the size of a problem processed with MOSEK is roughly five times smaller than the same problem solved with XA. Hence, for the same RAM memory, thinner mesh problems are solved with MOSEK, the conic optimizer. With 1.5G of RAM, XA is limited to roughly 8000 elements (with 64 linearization planes), whereas MOSEK is limited to roughly 18,000 elements. With MOSEK, the 14,000-element static problem (631,810 rows, 672,011 columns, 3,409,939 non-zero terms) is solved in 2 hours. M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 813 Fig. 15. Yield criteria with 3D general loading. Homogeneous boundary strain and stress rate. 14,000-element mesh. f = 6.4%. We have seen that in the 3D case with an axisymmetric loading, Gurson’s criterion is slightly improved. If our numerical approaches and this analytic criterion are correct, we must find the same plastic multiplier λ for each macroscopic stress state. The plastic multipliers are defined in the expression: Ei = λ ∂fGurson , ∂Σi λ  0, i = x, y, z. (28) Actually, using our numerical values in (28) , we find the three values of the plastic multiplier λ very close (0.1%) for Σm not too high. For Σm in the vicinity of its maximum, the difference between the two multipliers increases because two points on the same vertical are compared. This comparison also validates the elasto-plastic calculations using Gurson’s criterion. 4.5. Results with other loadings 4.5.1. General 3D loading For a general 3D loading, a 14,000-element mesh and a porosity rate of 6.4%, Fig. 15 shows the numerical criteria obtained under average strain boundary conditions (E-values, as in Gurson’s model) and average stress boundary conditions (quoted S-values). As coalescence implies a concentration of the strain rate in thin layers linking voids, conditions of homogeneous boundary stress are a more pertinent choice than conditions of homogeneous boundary strain to deal with cases where coalescence can be expected to occur. In the appendix, Table 1 gives numerical values of Σeqv versus Σm corresponding to Fig. 15. It is worth noting that the little difference between our results of Figs. 15 and 14 comes from the fact that, in the present 3D case, there is not additional continuous field. With the purpose of having the numerical values of the yield criterion for more current porosity rates, Figs. 16–19 represent our numerical criteria for f = 1% and 10% with the criteria of Gurson, Bigoni, Tvergaard and Richmond (with their standard parameter values) under 3D loading with average boundary strain and stress rates for the 14,000-element mesh. Tables 2 and 3 of the appendix provide numerical values of Σeqv as a function of Σm for these two porosity rates. For f = 10% (Figs. 16 and 17), static values are quite close to kinematic values. For f = 1% (Figs. 18 and 19) static and kinematic values are more distant since the meshing is not thin enough for this low porosity rate. A thinner mesh is not yet possible since MOSEK cannot process more than 2 RAM-Gbytes with a 32-bit system. All the numerical curves are located between the Gurson and Tvergaard criteria. For these two porosity rates, Tvergaard and Richmond criteria are too distant from the real criterion. Finally, the Gurson criterion is the closest analytic expression among all those quoted, making it possible to approximate the real criterion in this case. Fig. 20 represents the displacement velocity field in the z = 0 plane in homogeneous boundary strain rate for f = 6.4% and 14,000 elements (or 400 along the z = 0 plane) for two Σm values. We again find the displacement velocity field found previously in the axisymmetric method. Fig. 21 represents the displacement velocity field in the x = 0 plane in the homogeneous boundary strain rate for f = 6.4% and 14,000 elements (or 400 in the x = 0 plane) for two Σm values. As expected, the symmetry conditions can be verified on these figures. 4.5.2. Invariants in general 3D loading Let (I1 , I2 , I3 ) be the three invariants of the stress tensor and (J1 , J2 , J3 ) the three invariants of the deviatoric stress tensor. 814 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 Fig. 16. Yield criteria with 3D general loading. Homogeneous boundary strain rate. 14,000-element mesh. f = 10%. Fig. 17. Yield criteria with 3D general loading. Homogeneous boundary stress rate. 14,000-element mesh. f = 10%. Fig. 18. Yield criteria with 3D general loading. Homogeneous boundary strain rate. 14,000-element mesh. f = 1%. M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 815 Fig. 19. Yield criteria with 3D general loading. Homogeneous boundary stress rate. 14,000-element mesh. f = 1%. Fig. 20. Displacement velocity field in z = 0 plane. Homogeneous boundary strain rate, f = 6.4%, 3D loading. 14,000 elements. Field scale = 3. Fig. 21. Displacement velocity field in x = 0 plane. Homogeneous boundary strain rate, f = 6.4%, 3D loading. 14,000 elements. Field scale = 3. 816 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 Gurson’s criterion is written as a function of three invariants: 
I J fGurson = 22 + 2f cosh √1 −1−f2 k 2 3k (29) with:  I1 = tr Σ = Σx + Σy + Σz = 3Σm ,       2  1  Σeqv 1 J2 = − tr S 2 − (tr S )2 = (Σx − Σy )2 + (Σx − Σz )2 + (Σy − Σz )2 = ,  2 6 3      I = det Σ = Σ Σ Σ . x y z 3 (30) Since the Gurson’s criterion is isotropic, it depends only on Σm (linked to I1 ), Σeqv (linked to J2 ) and det Σ (I3 ). Leblond (2003) notes that the Gurson’s criterion is obtained by assuming it does not depend on det Σ . We want to check the dependence of the criterion in det Σ . The difficult point is that we cannot directly optimize I3 , nonlinear and non-quadratic or conic. To bypass this, Σgps is maximized with I1 fixed and for various values of Σx ; after solving, J2 and I3 are calculated from the optimal values of Σx , Σy , Σz . As a result, we see that J2 (and Σeqv ) are quasi-constant whereas I3 varies a great deal. Hence, at least in the case studied here, the criterion really seems independent of the third invariant I3 , and consequently, the expression of Gurson using Σeqv and Σm is valid. It was not the case for a media with cylindrical cavities: to fit our numerical criterion in shear loading, it was necessary to use all five invariants of a transverse isotropic media. 4.5.3. Macroscopic plane stress and strain For the purpose of analysing other features of the real criterion, we now examine the cases of macroscopic plane stress and strain. 2 = Σ 2 + Σ 2 with: In plane stress (i.e. Σz = 0 here), Σeqv is defined by Σeqv cp h 1 Σh = (Σx + Σy ), 2 Σcp = √ 3 (Σx − Σy ), 2 Eh = (Ex + Ey ), 1 Ecp = √ (Ex − Ey ). 3 For Σh fixed Σcp is optimized. The Gurson criterion for the spherical cavities here becomes: 

Σh 2 = 3 1 + f 2 − 2f cosh √ − Σh2 . Σcp 3k (31) 0 where Σ 0 and E 0 are fixed values. For a In the kinematic method, the power rate is written: Ptot /V = Σh0 Eh + Σcp Ecp cp h 14,000-element mesh and the porosity rate f = 6.4%, the curves in Fig. 22 are obtained using MOSEK as the solver. √ 3 2 = Σ 2 + Σ 2 with: Σ In plane strain (i.e. Ez = 0 here), Σeqv is given by: Σeqv gps = Σh − Σz , Σps = 2 (Σx − Σy ). gps ps The corresponding Gurson criterion is not analytical, the plane strain condition being given by: Fig. 22. Criteria in macroscopic plane stress (left) and strain (right). 14,000-element mesh. f = 6.4%. M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 817 
2 f (Σh + Σz /2) = 0. sinh (Σ − Σ ) + √ √ z h 3k 2 3k k 3 (32) Then, for each given Σh , the value of Σps is obtained after solving the resulting non-linear system in (Σx , Σy , Σz ) with Matlab, obtaining the plane strain Gurson criterion of Fig. 22. On other hand, using a 14,000-element mesh and a porosity rate f = 6.4% with Ez = 0 (recall here Exz = Eyz = 0 are implicit), static and kinematic values of Fig. 22 are obtained using MOSEK as the solver. As expected, the plane strain and plane stress Gurson criteria can be taken as an approximation of the real criterion, although it is a little not as good as in the general 3D case. 5. Conclusion After recalling the results obtained for porous materials with cylindrical cavities, we have proposed two complementary approaches, axisymmetric and 3D, to determine the yield criterion of a porous material with spherical cavities according to the Gurson model. In each of these two approaches, the static and kinematic methods of limit analysis have been used, via a discretization of the model in finite elements. In the case studied herein of a Mises material, it can be pointed out that the high degree of efficiency reached by our optimization programs is mainly attributable to recent optimization algorithms of the interior point type in both linear and conic optimization such as XA and MOSEK codes. For porous materials with spherical cavities, the Gurson criterion is checked and seems to be a good analytic expression contrary to the case of materials with cylindrical cavities where Gurson’s expression is too restrictive. More precisely for spherical cavities, the 3D static approach confirms the excellent quality of the results obtained with the axisymmetric kinematic method; moreover the latter is confirmed by the 3D kinematic approach. Therefore, the 3D static and the axisymmetric kinematic approaches (or the 3D kinematic approach) are efficient tools for forecasting the yield criterion of these porous materials, at least for a Mises material. Further investigation is needed by using other criteria, for example Tresca and Coulomb instead of Mises materials, resulting in more complex, non-conic programming problems. The study of RVE with several cavities is also in progress. Acknowledgements This work was partly supported by a financial support of Electricité de France (EDF, G. Rousselier). The support of the Région Rhône-Alpes for this work as an Emergence Project is also gratefully acknowledged here, as well as Professor J.-B. Leblond for fruitful comments about this subject. Appendix A. Numerical values of Σeqv under 3D loading with homogeneous boundary strain and stress rates for f = 6.4%, 1%, 10%. Table 1 Numerical values of Σeqv under 3D loading with homogeneous boundary strain and stress rates for f = 6.4% Σm 0 0.8 1.6 2 2.4 2.8 3 3.1 max Σm 6.4% Strain loading Gurson static kinematic Stress loading static kinematic 1.6212 1.5913 1.4821 1.3759 1.2057 0.9086 0.6465 0.4309 3.1741 1.5848 1.5589 1.4417 1.3241 1.1330 0.7942 0.5057 0.2618 3.1086 1.5911 1.5651 1.4532 1.3420 1.1621 0.8528 0.5978 0.4077 3.1915 1.4421 1.4276 1.2836 1.1661 0.9642 0.6791 0.3968 – 3.1085 1.4963 1.4749 1.3235 1.1922 1.0106 0.7396 0.5120 0.3613 3.1911 818 M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 Table 2 Numerical values of Σeqv under 3D loading with homogeneous boundary strain and stress rates for f = 1% Σm 0 1.6 3.2 4 4.4 4.8 5 5.1 5.2 5.3 max Σm 1% Strain loading Stress loading Gurson static kinematic static kinematic 1.7147 1.6949 1.5872 1.4286 1.2824 1.0410 0.8493 0.7177 0.5389 0.2129 5.3176 1.7 1.6784 1.5435 1.3452 1.1568 0.8450 0.5467 0.3856 – – 5.1746 1.7086 1.6894 1.5805 1.4271 1.2888 1.0670 0.9013 0.7966 0.6698 0.5084 5.4388 1.6937 1.6732 1.5016 1.2745 1.0705 0.7698 0.5472 0.2001 – – 5.1791 1.7048 1.6835 1.5511 1.3713 1.2214 1.0009 0.8467 0.7509 0.6355 0.4872 5.4385 Table 3 Numerical values of Σeqv under 3D loading with homogeneous boundary strain and stress rates for f = 10% Σm 0 0.8 1.6 2 2.2 2.4 2.5 max Σm 10% Strain loading Gurson static kinematic Stress loading static kinematic 1.5588 1.5100 1.3250 1.1319 0.9844 0.7712 0.6173 2.6588 1.4805 1.4624 1.2669 1.0373 0.8598 0.6315 0.4621 2.5650 1.5206 1.4774 1.2868 1.0820 0.9269 0.7106 0.5601 2.6316 1.3774 1.3430 1.1550 0.9436 0.7918 0.5973 0.4318 2.6229 1.4150 1.3849 1.1872 0.9912 0.8511 0.6608 0.5305 2.6704 References Ben-Tal, A., Nemirovski, A., 2001. On polyhedral approximations of the second-order cone. Math. Oper. Res. 26, 193–205. Bigoni, D., Piccolroaz, A., 2003. Yield criteria for quasibrittle and frictional materials. Int. J. Solids Struct. 41, 2855–2878. Chakrabarty, J., 1987. Theory of Plasticity. MacGraw-Hill. Francescato, P., Pastor, J., Riveill-Reydet, B., 2004. Ductile failure of cylindrically porous materials. Part 1: Plane stress problem and experimental results. Eur. J. Mech. A Solids 23, 181–190. Garajeu, M., 1995. Contributions à l’étude du comportement non linéaire de milieux poreux avec ou sans renforts. PhD Thesis, Université de la Méditerranée, France. Gologanu, M., Leblond, J.-B., Perrin, G., Devaux, J., 1997. Recent extentions of Gurson’s model for porous ductile metals. In: Suquet, P. (Ed.), Continuum Micromechanics. In: CISM Courses and Lectures. Springer, New York, pp. 61–130. Gurson, A.L., 1977. Continuum theory of ductile rupture by void nucleation and growth – Part I: Yield criteria and flow rules for porous ductile media. ASME J. Engrg. Mater. Technol. 99, 2–15. Kim, K.T., Kwon, Y.S., 1992. Strain hardening response of sintered porous iron tubes with various initial porosities under combined tension and torsion. J. Engrg. Mater. Technol. 114, 213–217. Leblond, J.-B., 2003. Mécanique de la rupture fragile et ductile. Etudes en mécanique des matériaux et des structures. Hermes Science. MOSEK ApS, Denmark, 2002. The MOSEK optimization tools version 2.5 (revision 51) – User’s manual and reference. Download on http: //www.mosek.com/. Pastor, J., 1983. Application de la théorie de l’analyse limite aux milieux isotropes et orthotropes de révolution. Thèse d’Etat, ISMG, Grenoble. Pastor, J., Francescato, P., Trillat, M., Loute, E., Rousselier, G., 2004. Ductile failure of cylindrically porous materials. Part 2: Other cases of symmetry. Eur. J. Mech. A Solids 23, 191–201. Pastor, J., Loute, E., Thai, T.H., 2002. Limit analysis and new methods of optimization. In: Mestat, P. (Ed.), Vth NUMGE (Vth European Conference Numerical Methods in Geotechnical Engineering). Presses des Ponts et Chaussées. Pastor, J., Turgeman, S., 1982. Limit analysis in axisymmetrical problems: Numerical determination of complete statical solutions. Int. J. Mech. Sci. 24, 95–117. Pastor, J., Turgeman, S., Boehler, J.P., 1990. Solution of anisotropic plasticity problems by using associated isotropic problems. Int. J. Plasticity 6, 143–168. Pellegrini, Y.-P., 2002. Plasticity criterion for porous medium with cylindrical void. C. R. Mecanique 330, 763–768. Richmond, O., Smelser, R.E., 1985. Alcoal Technical Center Memorandum, March 7. M. Trillat, J. Pastor / European Journal of Mechanics A/Solids 24 (2005) 800–819 Rousselier, G., 2001. Approche locale de la mécanique de la rupture: état de l’art, Rapport du projet GALA. Sun, Y., Wang, D., 1989. A lower bound approach to the yield loci of porous materials. Acta Mech. Sinica 5 (3), 237–243. Tvergaard, V., 1981. Influence of voids on shear band instabilities under plane strain conditions. Int. J. Fracture 17 (4), 389–407. 819