Numerical limit analysis bounds for ductile porous media with oblate voids

Mechanics Research Communications 38 (2011) 350–354

Contents lists available at ScienceDirect

Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom

Numerical limit analysis bounds for ductile porous media with oblate voids F. Pastor a , D. Kondo c,∗ , J. Pastor b a

Laboratoire de mécanique de Lille (LML), UMR 8107 CNRS, 59 655 Villeneuve d’Ascq, France Université de Savoie, POLYTECH’Annecy-Chambéry, Laboratoire LOCIE, 73 376 Le Bourget du Lac, France c Institut D’Alembert, UMR 7190 CNRS, UMR 7190 CNRS, UPMC, 4, Place Jussieu, 75252 Paris Cedex, France b

a r t i c l e

i n f o

Article history: Received 9 March 2010 Received in revised form 10 November 2010 Available online 9 February 2011 Keywords: Gurson-type models Spheroidal voids Micromechanics Limit analysis Upper and lower bounds Conic programming

a b s t r a c t The paper is devoted to a numerical limit analysis of a hollow spheroidal model with a von Mises solid matrix. To this purpose, existing kinematic and static 3D-FEM codes for the case of spherical cavities have been modified and improved to account for the model of a spheroidal cavity confocal with the external spheroidal boundary. The optimized conic programming formulations and the resulting codes appear to be very efficient. This framework is then applied to the derivation of numerical upper and lower anisotropic bounds in the case of an oblate void. The numerical results obtained from a series of tests are presented and allow to assess the accuracy of closed-form expressions of the macroscopic criteria proposed by Gologanu et al. (1994, 1997) for porous media with oblate voids. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction

and Pastor et al. (2004), the two LA approaches made it possible to numerically determine the yield criteria of a cylindrically porous material, proving also that the Gurson criterion is approximate in the plane strain case. On the contrary, in the subsequent work (Trillat and Pastor, 2005) the Gurson criterion is shown to be relevant for materials with spherical voids. The main advantage of these LA numerical approaches is that they provide rigorous lower and upper bounds to the macroscopic criterion together with their controllability a posteriori from the final optimal solution. This capability to control numerical or analytical results is central and was used in Pastor and Castaneda (2002), and in Thoré et al. (2009) for example. In the present paper, we briefly describe the extension of the 3D static and kinematic codes of Trillat and Pastor (2005) for a von Mises matrix containing a spheroidal cavity (confocal with the external boundary). Then, we compare the predictions of the two closed-form approximate expressions proposed by Gologanu et al. (1994, 1997), to the numerical bounds obtained for the macroscopic yield surface under uniform macroscopic strain rates imposed at the external boundary of the model.

The celebrated Gurson criterion (Gurson, 1977) for ductile porous media is based on a micro-macro approach and on the limit analysis (LA) theory. The Gurson approach is based on the consideration of a hollow von Mises sphere or cylinder subjected at its external boundary to a uniform macroscopic strain rate. Gurson’s analysis consists in the use of the LA kinematic approach of the hollow sphere in order to obtain an upper bound (as pointed out in Leblond, 2003) to the macroscopic criterion of the isotropic porous material, at least in the sense of the Hashin’s Composite Sphere Assemblage. A heuristic parametric refinement of Gurson’s model has been proposed later in Tvergaard (1981) and Tvergaard and Needleman (1984) to define the widely used Gurson–Tvergaard–Needleman (GTN) model. More recently, several extensions1 of the Gurson model have been proposed, the most probably important developments being those accounting for voidshape effects (Gologanu et al., 1997; Garajeu and Suquet, 1997; Monchiet et al., 2007). On the other hand, using finite element discretizations, both static and kinematic methods of LA have been elaborated for Gurson problems since Thai-The et al. (1998). In Francescato et al. (2004)

2. The hollow spheroid model

∗ Corresponding author. Tel.: +33 144275485; fax: +33 144275259. E-mail addresses: [email protected] (F. Pastor), [email protected] (D. Kondo), [email protected] (J. Pastor). 1 Mention can also be made of models taking into account plastic anisotropy (Benzerga and Besson, 2001; Monchiet et al., 2008).

The considered hollow spheroid model is made up of a single spheroidal cavity embedded in a confocal spheroidal cell. The solid matrix is an isotropic, homogeneous, and rigid-plastic von Mises material. Fig. 1 presents the geometrical model, where the given aspect ratio a1 /b1 and porosity f allow to determine the characteristics a2 and b2 of the confocal spheroidal boundary. Let us consider

0093-6413/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2011.02.001

F. Pastor et al. / Mechanics Research Communications 38 (2011) 350–354

351

Fig. 1. The hollow spheroid model (a1 /b1 = 0.5, f = 0.1).

the three-dimensional point of view, and note ˙ and E the macroscopic stress and strain rate tensors. These quantities are related to the microscopic fields by the average over the model of volume V: ˙ij =

1 V



ij dV

Eij =

V

1 V



dij dV

(1)

V

where d denotes the strain rate tensor whose components are dij = (ui,j + uj,i )/2, and u the velocity vector. Under the Hill–Mandel boundary conditions, here ui = Eij xj (in which x represents the position vector) on the external boundary, the overall virtual dissipated power Ptot = ˙ ij Eij can be written as follows: Ptot = V (˙m Em + ˙ps Eps + ˙gps Egps + ˙yz 2Eyz + ˙zx 2Ezx + ˙xy 2Exy ),

(2)

where the macroscopic stresses (the loading parameters in terms of limit analysis) and the associated strain rates here are defined as: ˙x + ˙y 1 (˙x + ˙y + ˙z ), ˙gps = − ˙z , 3 2 √ 3 = (˙x − ˙y ), ˙yz , ˙zx , ˙xy 2

˙m = ˙ps

Em = (Ex + Ey + Ez ), 1 Eps = √ (Ex − Ey ), 3

Egps

2 = 3



(Ex + Ey ) − Ez 2

2Eyz , 2Ezx , 2Exy .

material. From the matrix isotropy and the spheroidal geometry of the model, the resulting material is transversally isotropic around the axis z. Here, is investigated the macroscopic criterion g(˙) in the (Oxyz) anisotropy frame. To compare with Gologanu’s axisymmetric results, which are also 3D results, we search for the projection of g(˙) in the (˙ gps , ˙ m ) plane by optimizing ˙ gps for fixed average stresses ˙ m , other stress components defined in (3) being free. Then ∂g/∂˙ ij = 0 = 2Eij for / j, and ∂g/∂˙ ps = 0 = Eps since the macroscopic material complies i = to the normality rule. As a final result, loadings can be restricted to principal macroscopic strain rates E (as well as ˙ since (Oxyz) is a transverse isotropy frame) with Eps = 0. Moreover, all the axes in the horizontal plane of Fig. 2 are equivalent; therefore, in the previous projection problem we also impose ˙ x = ˙ y , although this is not mandatory as well as Ex = Ey in fact. Indeed, when non imposed a priori, these equalities are always verified in the optimal solutions, giving by the way a good control of the mesh quality. Finally, the overall external power Ptot here reduces to: Ptot = V (˙m Em + ˙gps Egps ).

(3)

 ,

(4)

In these definitions the subscripts (gps for generalized plane strain, and ps for plane strain) were defined in Pastor et al. (2004), as ˙ gps = 0 is the usual relation in plane strain for the von Mises

(5)

Therefore, adapting the technique defined in Trillat and Pastor (2005) for spherical cavities, the eighth of the hollow spheroid is meshed into tetrahedral elements as shown in Figs. 2 and 3. This mesh respects the symmetries of the problem as the vertical coordinate planes are equivalent regarding the distribution of elements, giving rise to a well conditioned numerical problem. Note that the macroscopic equivalent stress ˙ eqv is, in the present case, linked to ˙ gps by: 2 ˙eq v=

3 2 2 = (˙x − ˙z ) , dev(˙) : dev(˙) = ˙gps 2

where dev(˙) is the deviatoric part of ˙.

Fig. 2. General view and Oxz plane of a 896-tetrahedon mesh (a1 /b1 = 0.5, f = 0.1).

(6)

352

F. Pastor et al. / Mechanics Research Communications 38 (2011) 350–354

a conic constraint for the conic optimizer mosek:

 5 

x2 ≤ x6 = 2k. j

(9)

j=1

3.2. Numerical implementation

Fig. 3. Decomposition of the prism into two tetrahedrons and three pyramids of four tetrahedrons each (after Trillat and Pastor, 2005).

Hereafter, we first briefly present both lower and upper limit analysis approaches. Since the basis of the present numerical tools is its versions presented in Trillat and Pastor (2005) for spherical cavities, we only detail here the modifications and improvements implemented for the present case of spheroidal voids. In a second step, the obtained results are analyzed and compared to those provided in Gologanu et al. (1994, 1997) for oblate models under uniform strain rate boundary conditions.

3. Limit analysis: static method 3.1. The von Mises criterion expression As classically, the criterion is written as: f () =



J2

1 with J2 = tr(s2 ) and 2

1 s =  − tr()ı, 3

(7)

where ı is the second order unit tensor. Then in a (x, y, z) reference frame, the full 3D criterion reads:



2 √ 3



x + y − z 2

2

2

2

2

2

+ (x − y ) + (2yz ) + (2zx ) + (2xy ) ≤ 2k.

(8)

√ The constant k is the limit in pure shear, also given by 0 / 3 where  0 is the tensile strength of the von Mises material. It should be noted that (8) can be written, after obvious changes of variables, as

For each aspect ratio a1 /b1 , the inner and outer matrix boundaries of the spheroid mesh are adapted from the spherical case to obtain their confocal forms in the final mesh. Figs. 2 and 3 illustrate the case of a mesh with 4 layers of triangle based prisms, each layer constituted by 4 × 4 prisms of 3 × 4 tetrahedrons and 2 extremal tetrahedrons each, i.e., 896 tetrahedrons for the whole mesh. As the problem boundaries are not homothetic and from the practically finite number of triangles of the resulting polyhedral mesh boundaries, the porosity of the resulting mesh results not exactly the input one. Then, in a first step for each case of porosity and aspect ratio, the distribution of the angle ˛ (see Fig. 2 right) is optimized to precisely retrieve the desired porosity by progressively concentrating this distribution towards the more curved zone. In a final step, the radial distribution of the spheroid layers is also optimized to obtain the best value for the isotropic loading (˙ gps = 0) with the static code. The local stress field is chosen as linearly varying in x, y, z in the tetrahedral element, then with one 6-component tensor  located at each vertex. Consequently this stress field can be discontinuous trough any element boundary, a feature proven to be indispensable in the finite element static approach (Pastor, 1978; Bottero et al., 1980). Finally, to reduce the size of the constraint matrix of the numerical problem, a change of variables  → (x0 , . . ., x5 ) is now performed, where x0 = tr and x1 –x5 defined in (9), so that only the definition of x0 is needed as a new constraint. To get a microscopic stress field statically and plastically admissible, the following conditions are implemented: • Definition of macroscopic stress components ˙ ij as averages of ␴ij , plus the definitions of ˙ m , ˙ gps , ˙ ps , i.e. six conditions. • In each element, the equilibrium equations, ␴ij,j = 0, generate three conditions. • Continuity conditions: the stress vector is continuous across every discontinuity surface; each discontinuity triangle generates nine conditions. • Boundary conditions: the components of stress vector Ti =  ij nj are null at each apex of the triangles of the cavity boundary, here also generating nine conditions per triangle. • Symmetry conditions: the microscopic tangential stresses are null at the triangle apices of the coordinate planes, i.e. six conditions for each triangle. All the previous conditions form a matrix of linear equality constraints. To enforce the stress field to be plastically admissible, the criterion (8) is imposed at each apex of the tetrahedron; hence, due to its convexity, this criterion is fulfilled anywhere in the element. For each tetrahedron the four conic inequalities are directly enhanced by mosek by simply indicating the names of the variables xi involved in (9). The final numerical problem is a constrained conic programming one which can be handled by mosek. To define here the objective functional, ˙ m is initialized to the desired value and ˙ gps is optimized; when ˙ m approaches its maximum value, then ˙ gps is fixed and ˙ m is maximized for better convergence of the optimization process.

F. Pastor et al. / Mechanics Research Communications 38 (2011) 350–354

353

4. Limit analysis: kinematic method 4.1. Dissipated powers Let us recall that, from the virtual power principle, the total dissipated power Ptot here reads: ˙m Em + ˙gps Egps =

Ptot (Pvol + Pdisc ) = V V

(10)

where the volumic dissipated power Pvol = now defined as:



(d)=2k

 Vm

(d)dV and (d) is

2  √ 2 3 1 2 + d2 + d2 . + (dxx + dyy ) (dxx − dyy ) + dyz xy zx 2 2 (11)

The power dissipated by the velocity jump [u] on the discontinuities is given by:



Pdisc =



Sd



k |ut | dS =

(u)dS = Sd

k



ut1 2 + ut2 2 dS

4.2. Numerical implementation

• Incompressibility condition: the condition tr(d) = 0 in each element gives rise to one condition, as the velocity field is linear within the tetrahedron. • Continuity of the normal velocity: one condition un = 0 at each apex of a triangular discontinuity surface. • Symmetry conditions: one condition un = 0 on the apices of the coordinate planes of normal n. • Boundary conditions: three equalities ui = Eij xj for each apex of the external boundary triangles. All previous conditions finally give rise to a set of linear equality constraints. Concerning the definition of the functional to be optimized, by taking into account (11), we can upper bound the volumetric dissipated power in the tetrahedron by writing: (13)

for each tetrahedron of volume Vel . The first inequality in (13) gives one conic constraint and one non-negative auxiliary variable Y for each element. As the velocity jump [u] is linear onto the discontinuity side, we use the convexity of ([u]) to upper bound Pdisc by writing at each apex i = 1–3 of the triangular side of surface Sside : ([u]i ) ≤ Zi

˙m = min

ub /V ) − ˙ 0 E (Ptot gps gps

(14)

side ≤ where the Zi are new non-negative auxiliary variables, and Pdisc

S side (Z1 + Z2 + Z3 )/3, resulting in three conic constraints for each discontinuity side. Then, using these definitions and after integraub for tions over the mesh, we substitute the final upper bound Ptot Ptot in the following.

0 Em

,

0 with Em = 1.

(15)

0 ,˙ To obtain other non-zero (˙m gps ) points, we chose the following functional, for better convergence as in the static case:

˙gps = min

The previous mesh type is also used for the kinematic approach. The displacement velocity field is chosen as linearly varying in x, y, z in each tetrahedral element, and any triangular surface common to two tetrahedrons is a potential surface of velocity discontinuity. Then the variables are (ux , uy , uz ) velocity vectors located at the apices of each tetrahedron. The following conditions must be enforced so that a displacement solution field is admissible:

Pvelol ≤ V el Y

In the present case of two loading parameters, ˙ m and ˙ gps , the following functional is used to define the points of the macroscopic criterion for zero or small absolute values of ˙ gps :

(12)

Sd

where Sd is the set of discontinuity surfaces; ut1 and ut2 are the tangential displacement velocity jumps in an orthonormal frame (n, t1 , t2 ) of the discontinuity side of normal n.

(d) ≤ Y ;

Fig. 4. Comparison of the 1994 and 1997 Gologanu criteria with numerical bounds. The aspect ratio is taken as a1 /b1 = 0.5 and the porosity f equal to 0.01.

ub /V ) − ˙ 0 E (Ptot m m 0 Egps

,

0 with Egps = ±1.

(16)

Remark. For both static and kinematic methods the final conic problem is solved using mosek, and the admissible character of the optimal solution field is checked a posteriori. 5. Numerical tests and evaluation of existing macroscopic criteria As the von Mises criterion is even in term√ of stress, we limit the study to nonnegative values of ˙ m , denoting 3k as  0 in in the following. The main objective here is to present the numerical bounds obtained in the case of an oblate cavity and uniform strain rate (E) boundary conditions; two aspect ratios, a1 /b1 = 0.5 and a1 /b1 = 0.2, have been considered. In the static case the mesh involves 13 layers of 12 × 12 prisms (with triangular basis) of 14 tetrahedrons, resulting in a conic programming problem with 733,830 variables and 658,015 constraints. In the kinematic case the mesh is the same, but with 11 × 11 × 11 prisms resulting in 663,570 variables and 445,405 constraints. CPU times are about 3000 s in the static case and 5000 s in the kinematic one, using the release 5 of mosek on a recent Apple Mac Pro (with one core). In both cases the memory limitation does not allow to consider more refined meshes in the performed computations. The CPU times also explains that obtaining the full numerical yield criterion (then with meshes not limited to the height of spheroid, and for various other geometries) with the objective to define some corrective terms to the existing analytical criteria, seems not realistic, at least for the moment. To begin with the comparison of the numerical bounds with the analytical yield criterion of Gologanu et al. (1994) and of the second version in 1997 (Gologanu et al., 1997), we first present the solutions obtained for an aspect ratio a1 /b1 = 0.5, for two different values of porosity. In Figs. 4 and 5 are plotted the results corresponding to f = 0.01 and f = 0.1, respectively. It appears that the two numerical bounds are particularly close each other, giving a precise evaluation of the real, exact criterion. The second approximate yield criterion of Gologanu et al. is closer to the numerical upper bound than the first one, without violating the static approach.

354

F. Pastor et al. / Mechanics Research Communications 38 (2011) 350–354

6. Conclusion

Fig. 5. Comparison of the 1994 and 1997 Gologanu criteria with numerical bounds. The aspect ratio is taken as a1 /b1 = 0.5 and the porosity f equal to 0.1.

In the present study, we have extended previous 3D static and kinematic FEM codes for von Mises matrices to the case of a spheroidal cavity confocal with external boundary. Implementations of both kinematic and static approaches have been described and they deliver close upper and lower bounds to the real criterion surface, respectively. A series of tests concerning with an oblate confocal cavity have been carried out. It concerns two void aspect ratios with two values of porosity in each case. From these tests, it appears that the Gologanu et al. (1997) criterion significantly improves the previous one provided by Gologanu et al. (1994), and it is generally coherent with the present results without significatively violating the static approach. As a final conclusion, it turns out that this second criterion can be obviously selected as the best candidate of analytical approximate one for structural computations involving ductile damage. Acknowledgement This study has benefited of the financial support of the ANR program “MICROFISS (ANR-07-BLAN-0300)”. References

Fig. 6. Comparison of the 1994 and 1997 Gologanu criteria with numerical bounds. The aspect ratio is taken as a1 /b1 = 0.2 and the porosity f equal to 0.01.

Fig. 7. Comparison of the 1994 and 1997 Gologanu criteria with numerical bounds. The aspect ratio is taken as a1 /b1 = 0.2 and the porosity f equal to 0.1.

For an aspect ratio of 0.2, the results are plotted in Figs. 6 and 7. It appears first that the numerical bounds provide yield surfaces satisfactorily close each other, in particular for the lowest porosity. Moreover, for this low porosity, the first approximate yield surface given by Gologanu et al. (1994), based on axisymmetric fields proposed by Lee and Mear (1992), is largely beyond our upper bound and overestimates the exact criterion in a large part of the stress domain for a 0.1 porosity. A contrario the second criterion of Gologanu et al. generally is in agreement with our bounds, without nowhere significatively violating the static approach. Moreover, it can be noted that the Gologanu et al. (1997) criterion provides significant improvements of the first one.

Benzerga, A.A., Besson, J., 2001. Plastic potentials for anisotropic porous solids. Eur. J. Mech. A: Solids 20, 397–434. Bottero, A., Nègre, R., Pastor, J., Turgeman, S., 1980. Finite element method and limit analysis for soil mechanic problems. Comput. Methods Appl. Mech. Eng. 22, 131–149. Francescato, P., Pastor, J., Riveill-Reydet, B., 2004. Ductile failure of cylindrically porous materials. Part I. Plane stress problem and experimental results. Eur. J. Mech. A: Solids 23, 181–190. Garajeu, M., Suquet, P., 1997. Effective properties of porous ideally plastic or viscoplastic materials containing rigid particles. J. Mech. Phys. Solids 45, 873–902. Gologanu, M., Leblond, J., Perrin, G., Devaux, J., 1994. Approximate models for ductile metals containing non-spherical voids—case of axisymmetric oblate ellipsoidal cavities. J. Eng. Mater. Technol. 116, 290–297. Gologanu, M., Leblond, J., Perrin, G., Devaux, J., 1997. Recent extensions of Gurson’s model for porous ductile metals. In: Suquet P. (Ed.), Continuum Micromechanics. Springer Verlag, 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. J. Eng. Mater. Technol. 99, 2–15. Leblond, J.B., 2003. Mécanique de la rupture fragile et ductile Etudes en mécanique des matériaux et des structures. Hermes Sci.. Lee, B.J., Mear, M.E., 1992. Axisymmetric deformation of power-law solids containing a dilute concentration of aligned spheroidal voids. J. Mech. Phys. Solids 40, 1805–1836. Monchiet, V., Charkaluk, E., Kondo, D., 2007. An improvement of Gurson-type models of porous materials by using Eshelby-like trial velocity fields. C. R. Méc. 335, 32–41. Monchiet, V., Cazacu, O., Charkaluk, E., Kondo, D., 2008. Macroscopic yield criteria for plastic anisotropic materials containing spheroidal voids. Int. J. Plast. 24, 1158–1189. MOSEK ApS, 2010. C/O Symbion Science Park, Fruebjergvej 3, Box 16, 2100 Copenhagen Ø, Denmark. http://www.mosek.com. Pastor, J., 1978. Analyse limite: détermination numérique de solutions statiques complètes Application au talus vertical. J. Méc. Appl. (now Eur. J. Mech. A: Solids) 2, 167–196. Pastor, J., Castaneda, P.P., 2002. Yield criteria for porous media in plane strain: second-order estimates versus numerical results. C. R. Méc. 330, 741–747. Pastor, J., Francescato, P., Trillat, M., Loute, E., Rousselier, G., 2004. Ductile failure of cylindrically porous materials. Part II. Other cases of symmetry. Eur. J. Mech. A: Solids 23, 191–201. Thai-The, H., Francescato, P., Pastor, J., 1998. Limit analysis of unidirectional porous media. Mech. Res. Commun. 25, 535–542. Thoré, P., Pastor, F., Pastor, J., Kondo, D., 2009. Closed form solutions for the hollow sphere model with Coulomb and Drucker–Prager materials under isotropic loadings. C. R. Méc. Acad. Sci. Paris 337, 260–267. Trillat, M., Pastor, J., 2005. Limit analysis and Gurson’s model. Eur. J. Mech. A: Solids 24, 800–819. Tvergaard, V., 1981. Influence of voids on shear band instabilities under plane strain conditions. Int. J. Fract. 17, 389–407. Tvergaard, V., Needleman, A., 1984. Analysis of the cup-cone fracture in a round tensile bar. Acta Metall. 32, 157–169.