Limit analysis of unidirectional porous media

M~hanics Research Communications, Vol. 25, No. 5, pp. 535-542, 1998 Copyright © 1998 Elsevier Scioncc Ltd Printed in the USA. All rights reserved 0093.6413/98 $19.00 + .00

Pergamon

Pll s0093-6413(98)00070-6

Limit Analysis of Unidirectional Porous Media The-Hung THAI, Pascal FRANCESCATO and Joseph PASTOR Laboratoire Matrriaux Composites (LaMaCo), ESIGEC, Universit6 de Savoie 73376 Le Bourget du Lac Cedex, France (Received5 September 1997; acceptedfor print 5 February 1998) Introduction. Metal rupture often occurs during great plastic strain caused by nucleation and the growth of cavities and microfissures. It is therefore necessary to model the behavior of porous materials in order to describe the ductile rupture of metals, as Gurson [ 1], the pioneer in the field, has done by a kinematic approach within the theoretical framework of Limit Analysis. In this paper we have undertaken the study of the behavior of cylindrical porous material using numerical formulations based on both the kinematic and static approaches of the above-mentioned theory. We demonstrate that Gurson's results are incorrect in the case of plane strain and correct in the case of generalized plane strain, the latter case and the axisymmetric constituting a first 3-dimensional approach. Our results are also compared to those from criteria proposed by various authors, allowing for an unequivocal positioning of the latter values from the Limit Analysis point of view. Presentation of the oroblem and aoplication of Limit Analysis. Among the models of porous material behavior proposed by several authors where the actual porous material is replaced by an "equivalent" homogeneous material with an additional internal parameter (porosity rate3'), Gurson's [1] is the most widely accepted. It is expressed in its original form and for cylindrical cavity materials as follows: --

Z

eqv

2

f

(1) Ceq v = (1 +

where

Y'm

3f + 24¢)2 forplane strain, Ceq v = 1 foraxisymmetriccase.

is the average macroscopic

s t r e s s , ~'eqv

is the Mises equivalent macroscopic stress as

defined precisely in (12) and (13) and k the simple shearing limit of Mises material. This model was elaborated by idealizing the porous material as a single cavity within a homothetic cell composed of a rigid-plastic Mises material. As Gurson has done, let us now turn to the heterogeneity problem of volume given as ~ (Fig. l) by a homogenization technique where the boundary conditions are the average strain rate E set at the boundary, i.e., u = E.x. The set K~ of loads F admissible by structure ~e, where • defines the heterogeneity scale and G ~ is the plasticity criterion for the material, is defined by: K e = ( F [ 3 ~ S . A . underF, ~(y)E Ge(y_), V y E D.e}

(2)

Homogenization here is the replacement of the determination of the set K e by the determination of the homologous set of a homogeneous structure ~ postulated as equivalent. For this homogeneous structure the field of admissible stress is defined by G h°m 535

as

follows:

536

T.-H. THAI, P. FRANCESCATO and J. PASTOR G h°rn = ( __Z=<=o'>Ev I q ~ S.A. under E at the boundary, ~ • G e }

(3)

where a is the stress field in the elementary volume EV, Figure 1. G h°m is therefore obtained by the resolution of a Limit Analysis problem defined on the elementary volume. As in [1] we define the stresses and macroscopic strains by:

[

(4)

~Eij=(vij)=l fvVijdV=2-~ ~r(Uinj+ujni)dS where V is the volume of the EV, F its outer boundary and n the outer normal at F (Fig. 1 ). If we consider a field of kinematically admissible (K.A.) strain rates v , a statistically admissible (S.A.) stress field _ff and the boundary conditions foreseen here, then the equality of Hilrs macrohomogeneity asserts that Zij are loading parameters and Eij are kinematic parameters associated in the Limit Analysis sense. G~(y_)

Ghorn

r-J

f~e

VE

Figure.1.Homogenizationof the heterogeneousstructure. Due to lack of space we will detail here only the method developed in plane strain and refer the reader to [6,8] for the generalized plane strain case, for which the procedure is analagous, and for the three-dimensional case. The plane strain problem. As Gurson did, we examined a circular EV with a surface S and a contour F in plane strain. A Cartesian frame Oxyz is attached to the EV (Fig. 1), which is subjected to F at u = u d = E . x. Displacement velocity u in the EV, the stress fields and the macroscopic and microscopic strain rate fields are defined by: u x = u x (x, y) Uy = Uy (X, y) Uz = 0

::::}

f {E) = (Exx, Eyy, 2Exy) T

{X} = <~xx' Xyy, Xxy, Xzz)T

t

(G} = (GXX, (~yy, Gxy, Gzz)T

(5) {V} = (Vxx , Vyy, 2Vxy)T

Mises' isotropic plasticity criterion is expressed as: f ( g ) = tr_82- 2k 2 < 0

~

f ( g ) = (~x-Gy) 2+ (2'txy) 2 < (2k) 2

where s is the deviatoring part of the stress tensor t~ and k the shearing limit. I

(6)

LIMIT ANALYSIS OF UNIDIRECTIONAL POROUS MEDIA

537

In space {trx - cry ,2Z~y} it is represented by a circle of a radius of 2k, The criterion is then linearized externally to retain the external character of the kinematic method by an outer polygon with n sides: a xcos2-~- C~yCOS2-~ + 2XxySin-~- 2k < 0

, i = l to n

(7)

or internally to retain the internal character of the static method by an inner polygon: (I x cos2-~ - ffyCOS2-~ + 2X×ysin-~ - 2kcos~ < 0

, i = l to n

(8)

With regards to the kinematic approach, we used meshes with Lagrange Pl-type finite elements, with continuous velocity fields on the mesh (Fig. 2). The admissible nature of the strain rate field v is expressed, for the linearized criterion (7), by equations linking for each element i the components Vxx, Vyy and Vxy and the plastic multipliers ~,~i)of the element. Furthermore, given the symmetries of the states of the average strain considered, the study can be reduced to one quarter of EV with specific boundary conditions. Discretizing EV in Nt Lagrange P1 triangular finite elements, each of a surface Si (where i varies from 1 to Nt), the functional to optimize is given by the following expression [8]:

P = Z i= 1

2k

~.~i) Si r=1

(9)

The kinematic method means solving, for different E ~, the following minimization problem: S*~: E d = Ph°m(Ed) <_ Min { P / v admissiblel, = E } (10) The expression of ph°m(E) / S*, S* being the total surface of EV, it therefore takes on the following form: p h o m ( ~ ) / S* = 3 ~ m E m + ~dpEdp + ~ t E t

(11) with: -

2

' Em =

3

'

'

- ~

Exy

(12)

The ~eqv and Eeqv expressions, necessary for comparison with Gurson, are being verified in the present case:

4 1 (Exx + Eyy)2+ 3 (Exx - Eyy)2 +3E2y] = E 2 + E 2 p+E2 Ec2qv= 9~" 2 = Z~p + Z 2 E~q~ 3 (Zxx - Eyy)2 + 3Zx,

(13)

538

T.-H. THAI, P. FRANCESCATO and J. PASTOR

After calculations [6] we obtained the following result: (P"°m(l~)/ S*- 3Y-mEre)2

In order to apply the kinematic method, assume that we have a set of fixed macroscopic strain rate directions E~. For each direction, the kinematic programme combined with linear optimization software gives the value of linearized dissipated power (1 l). It is then a question of finding, for each fixed •o and in the space limited externally by the kinematic hyperplanes within an appropriate frame, the point such that Eeqv is the smallest, a point situated on one of the hyperplanes [6]. The goal is indeed to trace, point by point, the approach obtained in the plane (Zm/k,Xeq,/k). For each value E °, the kinematic value associated with r.q~ is determined solving: ~'ph°m/ S* - 3X'0mEm } - - , ~ =~ •

Zeqv( xO ) = Min (

Ei~

~Ed2p+ E2

(15)

This allows a direct comparison with Gurson's results. With regard to the static method, at each apex p of a finite element is assigned a vector (&} (in numerical notation) of which the three components are those of the tensor cr :

((~p) = ((~xp,•yp, ~xPy)T

(16)

The components of the macroscopic tensor Z, which are here our three loading parameters E x, Zy, X~., are calculated directly by the average of the cr0. The conditions which make the field tr admissible in each finite element are composed of statistically admissible (S.A.) conditions (equilibrium equations, continuities of an and 1;ntbetween adjacent elements and boundary conditions) and the plastically admissible (P.A.) conditions (see [9] for more details). The Mises criterion is internally linearized according to the expression (8). On the boundary of the hole, and therefore at the apexes of each concerned element, the stress vector is set to zero.

I

Mesh M1

Mesh M2

Figure 2. Meshes for the elementary cell of porous materials

LIMIT ANALYSIS OF UNIDIRECTIONAL POROUS MEDIA

539

In the following computations the number of planes linearizing the criterion n is equal to 48. Numerical tests show that a greater value of n does not significantly improve the results given the increase in CPU time. Kinematic and static results. Plane strain

The kinematic meshes used for this problem are those of a quarter circle presented in Figure 2 with a porosity rate f equal to 0.16 for the sake of convenience. The number of planes linearizing the criterion is set at n= 48 and the number of fixed directions E~ is equal to 15. The resulting problems of linear optimization are solved with IBM's OSL optimization code on SUN SPARC stations. For the static method the construction of the macroscopic plasticity domain has been done point by point and the convex envelope of these different points is internal to G horn. Since Zxy is taken as zero on the total EV, ~'~eqvi s reduced to Edp = 4~" (Z x - Zx)/2, the functional to be maximized for each fixed ]~m " Problems are solved as above on SUN stations. The M1 and M2 meshes (Fig. 2) which were made discontinuous are used in the static calculations, as the number of planes linearizing the criterion is also set at n = 48. Table 1 displays the comparison between the numerical predictions obtained from the two meshes in Fig. 2 and those given by Gurson in the case of plane strain.

Zmlk

0.000 0.250 0.500 0.750 1.000 1.250 1.500 1.741 1.750 1.794 1.833 1.867 1.898

GURSON 0.9827 0.9757 0.9539 0.9147 0.8532 0.7594 0.6111 0.3247 ................. .0000

Eeqv [ k SmrC MI 1.3111 1.2933 1.2378 1.1426 0.9999 0.7722 0.4376 0.0000 .........

Zeqv / k static M2 1.3186 1.3028 1.2533 1.1669 1.0352 0.8263 0.5179

0.0912 0.0000

X,qv/ k kinematic M1 1.3953 1.3953 1.3397 1.2609 1.1265 0.9441 0.7005 0.3655 .................

~'eqv Ik kinematic M2 1.3733 1.3733 1.3156 1.2368 1.1039 0.9124 0.6423

0.2784

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

................................. .........................

0.0000

0.0000 .........

Table I. Comparison in plain strain case with meshes M1 and M2 ( f = O.16, n = 48).

The graphs in Figure 3 show the results obtained for Zeqv/k as a function of ]~m]k and those of the other authors cited above. It can be noted that there is good agreement between the kinematic and static results computed by our method. On the other hand, significant differences arise between our results and those given by Gurson, and also with Guennouni's criteria which were obtained, let us not forget, for a periodic square EV. Note that when Y.m/k is weak, Gurson's criterion, although

540

T.-H. THAI, P. FRANCESCATO and J. PASTOR

obtained for a kinematic method, underestimates the real criterion. For Ym/k = 0 for example, the differences are on the order of 25% (1.318 as opposed to 0.983) between our static value and the value calculated by Gurson, lower despite its being kinematic. This substantiates a result announced in [5] concerning the fact that Gurson's criterion transgresses a Hashin-Strikman-type bound. An admissible explicit stress field obtained for Zm/k = 0 by the static numerical method presented in [8] for a mesh with only 20 Lagrange P1 triangles confirms the gap between the two methods (1.155 as opposed to 0.983).

l 1.5

1

J i................................... '................................... i ii' : . ......... ...... -

~ ~

-

.......

_ "

i

i

i ...................................

i ...................................

. -

. ...................

i

1 1

~

.................................

~. . . . . . . . . . . . . . . . . . . . . .

..................................

i ........

~- . . . . . . . . .

.

Needleman Tvergaard& Needleman

..... n •

. . . . . . . . . . . :.,,.;.......~ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

:

Ouennouni

static values kinematicvalues

.4

~. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.:

i

0.5

L

,

'

0

0.5

i

i

i

1.5

~

N._~k'

~, rL~ I

l~llt

~.

i-

i "

2

I

t

~,m

lh

2.5

Figure 3. Comparison of results in plane strain (f = 0.16).

Generalized plane strain As above, the detail of the results and the tables of values are presented in [6] ; figure 4 shows that Gurson's axisymmetrical result is in fact now situated beyond the static approach and is often improved by our kinematic approach. Note that the present and the axisymmetrical approaches can be compared if both are considered from a three-dimensional point of view. The general problem is a three-dimensional mechanical problem (six components for the stress and strain rate tensors) with a plane mesh [ 11 ]. The generalized plain strain case (with only four components for each tensor) leads to a solution of the real problem : the obtained static and kinematic solutions are easily extended to three dimensional case.

LIMIT ANALYSIS OF UNIDIRECTIONAL POROUS M E D I A

]~eqv]k "t .................................... i................. r.....

2

......... i ....... i ..... ............................................................... ~................... •................... ' '" ~ , • .

1.5

-

Gurson (axisymetric) Tvergaard Koplick & Needleman Tvergaard&Needleman Guennouni kinematicvalues

[]

................~.................i...............",...~.~....... ~ ~ ..........

0.5

•, ~

................

g.

~

i

:

:

i

: ................

" ................

" ................

~ ..............

0

0.5

F i g u r e 4. C o m p a r i s o n

541

s,auc

,uos

....i.................!.............. i.................: . ,,

"r..

k ............

:

"-

1.5

' ................

.: .................

2

i

2.25

o f t h e r e s u l t s in g e n e r a l i z e d p l a n e s t r a i n ( f = 0 . 1 6 ) .

Discussion and conclusion The main goal of this paper concerns the verification, in particular by the static method, of the known results for cylindrical voids in a given case of porosity. Extending this work to other porosity rates would be desirable (this work is in progress), using fairly large ratios however because of numerical difficulties stemming from too much of a difference in the size of mesh elements First of all, it is important to note that the kinematic and static values obtained by our calculations are always coherent, as well as being systematically verified: an example of an explicit stress field confirming our conclusion for plane strain can be found in [8]. It can be seen that being placed in generalized plane strain is beneficial with regard to the axisymmetrical hypothesis which, it is true, allows analytical computations more easily, such as those carried out in [10]. A possible extension of this work is a finer analysis of the continuous displacement velocity fields (less costly numerically) to more precisely define the Gurson expression and examine the appearance of strain concentration zones. Then or in parallel, with a fine d i s c o n t i n u o u s

kinematic

mesh (but continuous at the outer boundary), it would be possible to observe where the origins of discontinuities occur, which would announce an onset of intercavity fissuration, hence the beginning of coalescence. Using the static method the analysis of stress fields would give complementary information on the distribution of plastified zones, providing that, as in plane strain, the static and kinematic values are sufficiently close. In conjunction with the kinematic approach it would be possible to deduce information on the appearance of shearing bands promoting coalescence of cavities. Finally, a study at different porosity rates should be carried out in order to thoroughly verify the validity of the proposed criteria, as the influence of this rate seems important from certain critical values that our r i g o r o u s approaches would be better able to define.

542

T.-H. THAI, P. F R A N C E S C A T O and J. PASTOR

References

[1 ] Gurson A.L. - Continuum Theory of Ductile Rupture by Void Nucleation and Growth : Part I J. Eng. Mat. Tech., 1977, Vol. 99, pp. 2-15. [2] Tvergaard V. and Needleman A. - Analysis of cup-cone fracture in a round tensile bar - Acta Metall., 1984, Vol. 32, pp. 157 - 169. [3] Tvergaard V. - Influence of voids on shear band instabilities under plane strain conditions - Int. J. Fracture, 1981, Vol. 17, n ° 4, pp. 389 - 407. [4] Koplick J. and Needleman A. - Void growth and coalescence in porous plastic solids - Int. J. Solids Structures, 1988, Vol. 24, n ° 8, pp. 835 - 853. [5] Leblond J.B., Perrin G. and Suquet P. - Exact results and approximate models for porous visco-

plastic solids - Int. J. Plasticity, Vol. I0, n ° 3, 1994, pp. 213-235. [6] Thai T.H. - Analyse limite : application aux structures et aux mat~riaux poreux - Th~se de Doctorat de l'Universit6 de Savoie, f6vrier 1997. [7] Guennouni T. - Fronti&e d'~coulement des mat~riaux h~t~rog~nes ~ constituants rigides

parfaitementplastiques. Cas des mat~riauxporeux oufissur~s - J. de M6c. Th. et Appl., 1987, Vol. 6, n ° 4, pp. 571-615. [8] Thai T. H., Francescato P., Pastor J. - Analyse Limite : application aux mat&iaux poreux - Proc. Int. Conf. Engineering Mechanics Today, 1997, Hanoi, Vietnam. [9] Pastor J. - Analyse limite : d~termination num~rique de solutions statiques completes. Application

au talus vertical - J. M6c. Appl,, 1978, n ° 2, pp. 167 - 196. [10] C I S M courses and lectures- Continuum Micromechanics- edited by P. Suquet, 1996, Udine. [11] Francescato P. and Pastor J. - Lower and upper numerical bounds to the off-axis strength of

unidirectional fiber-reinforced composites by limit analysis methods - Eur. Jl. of Mech. - A / Solids, 1997, Vol.16, n ° 2, pp. 213-234.