Microsoft Word - 5-509-ed.doc Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 734-738 734 www.etasr.com Siad and Gangloff: A GTN-like Model for Plastic Porous Materials A GTN-like Model for Plastic Porous Materials Larbi Siad EA 4691, University of Reims France larbi.siad@univ-reims.fr Sophie C. Gangloff EA 4691, University of Reims France sophie.gangloff@univ-reims.fr Abstract— An extended version of the well-known Gurson- Tvergaard-Needleman (GTN) isotropic hardening model is presented in this paper. The yield function of the proposed constitutive model possesses the distinctiveness to explicitly depend upon the third stress invariant. The presented constitutive model is used to analyze the necking of a round tensile bar. As long as softening initiation of specimen is not reached, the obtained numerical results highlight similarities and good agreement with those provided by the use of the GTN model. However, discrepancy shows up as soon as specimen failure starts. Keywords-constitutive modelling; ductile failure; necking; porous microstructures; tensile test; numerical simulations I. INTRODUCTION Regarding the mechanical behavior of porous plastic materials, the first micromechanical model which introduces a strong coupling between deformation and damage is the Gurson-Tvergaard-Needleman (GTN) model [1, 2]. The void volume fraction (vvf) f of the material at hand affects the yield function of this model and subsequently the ductile fracture behavior of this material. Indeed, when f increases, the material softens and loses its capability to carry loads An extension of the GTN’s plastic potential ΦGT is  proposed in [3] where the authors, investigating three arrangements of spherical voids over a wide range of f, focused their study on the determination of yield surfaces for porous materials using a huge number of finite element simulations. For each of the considered three void arrangements, namely Simple Cubic (SC), Body-Centered Cubic (BCC) and Face-Centered Cubic (FCC) arrays, many thousand yield points were determined by monotonically increasing the arbitrarily-prescribed macroscopic strain. The obtained yield points were fitted by a new yield function which turned out to be similar to ΦGT for f ranging between a very small value to the percolation threshold fp found to be equal to 0.89. Moreover, all introduced parameters may be expressed as functions of the parameters characterizing the GTN’s plastic potential ΦGT which, in addition, is retrieved from the plastic potential Φ when some of its parameters  are set to be equal to appropriate values (refer to  Section 2.). As a novelty, the yield function  Φ  was found to explicitly depend upon the third stress invariant I3. In this investigation, a constitutive model, referred to the present model in this paper, is used to describe the progressive damage in ductile solids resulting in material softening and, ultimately, loss of its stress carrying capacity. Numerical aspects were addressed concerning the integration of the proposed constitutive rate equations [4]. Based on the backward Euler integration scheme, a numerical algorithm implicit in all variables and a corresponding algorithmic operator was developed. In order to demonstrate the global accuracy and stability of the numerical solution, finite element damage simulations accounting for finite strain and using the proposed model are performed for the traditional ductile solid problems of necking of a round tensile bar [5]. II. THE PROPOSED RATE-INDEPENDENT PLASTICITY MODEL FOR POROUS MATERIALS The most popular model describing constitutive elastic-plastic equations that account for the effect of ductile damage development is that suggested by Gurson and phenomenologically improved by Tvergaard [1, 2]. For a material containing a volume fraction of voids the proposed approximate yield condition reads 2 2 22 GT 1 1 3 q pp (σ ; H)= +2 q csch - -1-q f =0 σ 2σ             (1) where σ is the effective flow stress of the damage-free matrix material which is a function of the effective plastic strain pε , (q1, q2) are the Tvergaard parameters, and H=(H1, H2) is a vector comprising the scalar state variables p 1H =ε and 2H =f . In yield condition (1), the macroscopic Cauchy stress tensor is resolved as 2 3 ' p I q n with n 3 2 q       (2) with p = -1 /(3 tr σ ) representing the hydrostatic pressure, tr refers to the trace of a tensor, I is the second order identity tensor, I is the deviatoric stress tensor, and 1 / 2q = (3 /2σ ':σ ') is the von Mises stress. An extension of the GTN’s plastic potential ΦGT with no extra parameters which fits the numerical data well and is valid for all void volume fractions and triaxial stress states has recently been proposed in [3]. For more detailed explanations of the subject, we refer the reader to this paper. Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 734-738 735 www.etasr.com Siad and Gangloff: A GTN-like Model for Plastic Porous Materials The porous ductile materials studied in [3] contain spherical empty voids arranged in cubic arrays, namely, simple cubic (SC), Body-Centered Cubic (BCC) and Face-Centered Cubic (FCC) arrays. FCA was used to simulate unit cells and the macroscopic yield surfaces of the porous materials were obtained using the probing technique which goal is to obtain a yield function in an analytical expression that can be used in continuum studies. The matrix material is almost rigid, perfectly plastic and unit cells were meshed with cubes. Depending on the unit cell at hand, the void volume fractions considered range from 0.02 to around 0.89 corresponding to the percolation threshold of the matrix material. Let I3 and J3 be the determinant of the stress tensor σ and the third stress invariant of the deviatoric stress tensor σ΄ respectively: 3I d e tσ and 33 1 J trσ ' 3  are related by 2 3 3 3I = J + p q /3 -p . In stress space (p, q, I3), yield points were found by monotonically increasing macroscopic strain with fixed ratios until the macroscopic equivalent stress reaches a maximum. For each of the three cubic unit cells, a least-squares fit with an extension of the GTN’s yield function (1) was found only approximately with deviations becoming more pronounced as pressure or void volume fraction were increased. To account for the observed tri-lobed numerical data pattern, the Gurson-Tvergaard-like yield criterion proposed in [3] can be expressed as 2 22 1 1 3 23 2 1 13 3 3 pp (σ; H) 2 cosh 1 2 I 3 pp p 1 p s(f ) 2 cosh 1 0 3 2                                               (3) Note that for the extended yield function Φ, the q-like parameters α1 and α2 depend on f for the three considered microstructures, that is α1=α1(f) and α2=α2(f). The yield function Φ is linear with respect to I3, with coefficient proportional to the hydrostatic pressure p. The parameter s also depending on f, determines the influence of the new term in the yield condition (3) which reduces to that of the GTN model (1) for s=0, α1=q1 and α2=q2. Whenever the constant s is non-zero, there is an effect of the third stress invariant I3 on the plastic flow. Clearly the yield function contains three functions depending upon f which are slightly different for each of the three cubic microstructures. In [3], the authors stated that their finite element data were well-fit by expressing the parameters α1 and α2 as n m 1 2 p p p p f f f f A , Y Z f f f f                                         (4) where n, m, A, Y and Z are constant parameters. Note that (1) for f=0 and α1=0, the yield condition (3) is nothing else than the classical Von Mises one; and (2) when f=fp and α1=1, the material becomes plastic under any applied load since the matrix material is no longer connected. Regarding the parameter s, if its value is too large, the tri-lobed pattern of the macroscopic yield surface almost degenerates into an equilateral triangle at some values of pressure p. In non-linear solid mechanics, the material behavior is often described by a rate-form constitutive equation. In the particular case of the proposed rate-independent plasticity model, it is assumed that the macroscopic rate of strain tensor ε is written as the sum of an elastic part eε and a plastic part pε , and the stress rate σ depends linearly on the elastic strain rate tensor eε . Subsequently, the corresponding constitutive equations can be written in a rate format as: e p e p p p ε= ε + ε , σ = C : (ε-ε ) ε = λr(σ ;H ), H = λh (ε ;H )         (5) where eC is the elastic moduli tensor which can depend on the current stress state, r is the direction of the plastic flow which depends on the current stress and on a finite set of plastic internal variables H accounting for history effects and h is the direction of the rate of these plastic internal variables. As previously mentioned, H characterizing the current state of the ductile porous material includes the void volume fraction f and the effective plastic strain pε In associated plasticity, r is the gradient of the yield function, that is: ij 3 ij ' ' ' 2 ij ij ik ikj ij 3 3 r ( ; H) (p, q, J ; H) 1 3 2 q 3 p 2 q q J 9 J                             (6) Then the plastic strain rate p , is given by: . . p ' 2 ij 3 1 2 I q I 3 p ' 9 J                    (7) The plastic strain rate pε is trivially decomposed into volumetric p v ε and deviatoric p q ε parts, p p p v q ε =1/3ε I+ε ,   which facilitate development of the integration algorithm. Indeed, these two terms can be expressed in terms of the yield function, using (6), as: . . p v p        and . . p 2 q 3 2 q I ' 9 J           (8) It should be noted at this stage of calculation that, due to the presence of the third stress invariant in the expression of the yield function Φ, and in view of (2), the deviatoric component cannot be put in the form p p q v ε =ε n  where n is the deviatoric Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 734-738 736 www.etasr.com Siad and Gangloff: A GTN-like Model for Plastic Porous Materials strain rate tensor normal to the yield surface Φ=0 and which norm is unity. As a reminder, this form turned out to be very successful for applying the implicit integration scheme based on the Aravas’s method [6]. To determine the plastic multiplier, the loading-unloading conditions should be imposed in a Kuhn-Tucker forma: . . 3 30, (p, q, J ; H) 0, and (p, q, J ; H) 0     (9) implying that during plastic loading Φ=0, . 0  and Φ 0 . This condition, the consistency condition, allows the determination of the plastic multiplier λ which specifies the magnitude of the plastic strain rate [7]. Finally, the void fraction modification f* introduced in [1] was intended to simulate the rapid loss of strength accompanying void coalescence. The calculated porous plastic material response is strongly dependent on the computational procedure for stress calculation, usually called the stress integration. The presence of the third stress invariant in the yield function typically results in a high degree of non-linearity and complex numerical algorithms. For further details dealing with the description of various integration methods, the reader is referred to [7, 8] where more elaborate coverage of the Closest Point Projection Method may be found. The constitutive equations of the proposed model along with the coalescence criterion based on the effective porosity introduced in [5] are integrated using a fully implicit scheme which allows the calculation of the consistent tangent operator [7, 8, 4]. The final stress and hardening parameters are determined solving the non-linear equations iteratively so that the stress increment fulfills the consistency condition. The corresponding algorithm is programmed in appropriate user subroutines which are called by Abaqus, at each element integration point, for each increment, and during each load step. III. ONE NUMERICAL EXAMPLE: SMOOTH BAR TENSILE TEST This section illustrates the use of the present model through the analysis of damage and failure of smooth bar subject to axial tensile test. The traditional necking analysis of a round bar has been investigated. Uniaxial tensile test of a smooth bar is often used to study the ductile void growth in metals for which the formation of a neck triggers the material instability. A detailed finite element analysis of the formation of the crack at the center of the neck and its subsequent growth leading to cup-cone fracture can be found, for instance, in [5]. The formation of the neck results in a triaxial state of stress at the center of the specimen. Due to high hydrostatic tension, the growth of nucleated voids is accelerated and the material instability proceeds with the initiation of fracture at the center of the neck with linkage of adjacent voids. Numerical simulations of the behavior of a rounded smooth bar subject to tensile test have been performed. The rounded bar has an initial length of zo 2L 8L   where the radius L  is assumed to be equal to one unit. Due to symmetry about ρ=0 and z=0, only one quarter of the specimen needs to be analyzed. A slight conical imperfection is introduced near the center of the specimen such that z 0 L 0.998L     . The geometry of the studied specimen and the original used mesh are shown in Figure 1. The geometry is axisymmetric, and the deformation is assumed to be axisymmetric. Fig. 1. Tensile test: geometry of the rounded smooth bar, original used mesh, boundary conditions, and loading. Reduced integration elements (CAX4R within Abaqus) were used for efficient computations. A refined mesh is generated near the center of the rounded bar where softening and intense deformation are expected, whereas a relatively coarse mesh is used in the rest of the specimen where a rather uniform deformation is expected. The kinematic boundary conditions are: symmetry about ρ=0, symmetry about z=0, and all the nodes on the edge z= z L  of the specimen are given a prescribed velocity depending on time z z v v (t ) in the z- direction. The velocity z v increases rapidly from zero at the start of the analysis to its maximum 30 ms-1 in 0.0025 s and drops back linearly to zero in another 0.0025 s. It then remains at zero for the remainder of the analysis. The specimen is made of a steel alloy which work- hardening characteristics are similar to those of a power law hardening material with Young’s modulus E, Poisson’s coefficient v, and initial yield stress σ0. Using the GTN model to simulate ductile failure, a proper choice of its eight micro- mechanistic parameters proves to be crucial. Before loading the material is assumed to be fully dense, that is f0=0.0. The damage parameters N, s and fN related to void nucleation, as well as the Tvergaard constants q1 and q2 were fixed to typical values suggested in [1]. Hereafter are the material parameters used in the numerical simulations: 0 0.00333   , 0.3  , p 0.1 0 (1 0.033 )     , q1=1.5, q2=1.0, 0.3  , 0.1s   , fN=0.04, fc=0.15 and ff =0.25. For the problem under consideration, it is well known that the material deforms homogeneously in the axial direction before the load becomes maximum, after which deformation becomes heterogeneous and necking starts. As a result of which the specimen softens due coalescence of voids and eventually fracture across the neck region. In this context, Figure 2 displays the deformed shape of the tensile bar showing necking occurs at the mid-section of the specimen. Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 734-738 737 www.etasr.com Siad and Gangloff: A GTN-like Model for Plastic Porous Materials Fig. 2. Comparison of vvf contour plots predicted by the present models and the GTN model just before failure of the tensile specimen; . Fig. 3. Simulation of tensile tests on smooth round bar. Predictions of overall nominal stress versus overall nominal strains based on the present models (SC, BCC and FCC microstructures) and GTN model. The global response of the specimen may be specified by the average states of strain and stress within it. The nominal elongation n z z zo L / L   and the relative diameter n z L / L       are chosen to assess the average strains of the specimen. The nominal elongation is used to assess the average strain in sections perpendicular to the tensile direction. The overall nominal stress is given by n 2 z F / L    where F is the resulting overall tensile load carried by the bar and determined from nodal forces. The evolution of n z  as functions of the nominal elongation and the relative diameter are shown in Figure 3, for both the present and the GTN models. First of all, these figures clearly show that up to failure initiation points which are indicated in Figure 3 using the dyadic  symbols, the agreements between predictions of considered models are quite good. In particular, maximum loads for the present model, regardless the considered microstructure, coincide with the maximum load provided by the GTN model which is reached at strains of about 11% for the nominal elongation and about 5.5% for the relative diameter. It is worth noting that the practically coincidence of the numerical predictions also holds in the softening regime ranging from maximum load and the immediate vicinity of failure points. This fact is substantiated by the deformed configuration and the void volume fraction contours depicted in Figure 2 at the final edge displacements of the specimen which are very similar. As expected, the largest amount of softening localizes around the center of the mid-section of the specimen. It is implied, therefore, that failure initiation should be expected in that area. However, the onset of failure is reached at displacements of 1.12Lρο for the GTN model and of 1.21Lρο for the present model with BCC microstructure, (Figure 2). The maximum loads along with the corresponding nominal strains predicted by both the present models and the GTN model are listed in Table I. TABLE I. MAXIMUM LOADS AND CORRESPONDING NOMINAL STRAINS. n z (%) nz 0 max  n (%) n 0 max  GTN model 8.67 1.301 4.87 1.302 SC 10.27 1.302 4.87 1.303 BCC 9.72 1.303 5.02 1.30 FCC 9.45 1.30 4.86 1.301 Deviations between predictions of the present models and the GTN model are observed as soon as failure of the considered specimen starts, that is when the sudden specimen capacity loss occurs (Figure 3). The numerical results seem to suggest that the present model based on SC microstructure is the closer one to the GTN model. The true fracture strain is the most physical experimental indicator for quantifying the resistance to damage and fracture of a material subjected to tensile loading conditions. An average true fracture strain measure f ε using the diameter reduction is usually used to quantify the ductility [9, 10]. It is defined as f ε 2 log(L / L )    It has been well-known that the dependence of f ε on the stress triaxialit τ can be approximated Engineering, Technology & Applied Science Research Vol. 4, No. 6, 2014, 734-738 738 www.etasr.com Siad and Gangloff: A GTN-like Model for Plastic Porous Materials using an exponential function f ε exp( )   where α is a material constant. The stress triaxiality is the hydrostatic stress to equivalent stress ratio. The obtained failure strains are listed in Table ΙΙ. Specimen made up of GTN material has a little lower ductility than the specimen made upof present models material, regardless the considered microstructure (BCC, FCC or SC). TABLE II. COMPARISON OF PREDICTED FAILURE STRAINS. DUCTILITY GTN SC BCC FCC n z f( ) (%) 27.5 28.3 30.2 30.5 n f( ) (%) 36.2 37.2 42.4 42.3 IV. CONCLUSION The main objective of this paper has been to address an extended version of the GTN model based on a pre-existing yield function for porous plastic materials proposed in [3] and its implementation within a finite element code. A fully implicit stress integration scheme is used to integrate the present model. As a numerical example, the present model is used to perform axisymmetric simulations of the ductile necking failure of a smooth round bar. Similar values for the damage parameters (nucleation and coalescence) have been used for both the present model and the GTN model in order to compare their ability to predict void growth to coalescence and the corresponding failure mechanism. The obtained numerical results are summarized as follows: • The present model has shown the potential to accurately predict all the major features of ductile failure of solids (stress triaxiality dependence …). • Up to the failure initiation of the specimen, the predictions based on the finite element analysis incorporating the present model are in a close agreement with those provided by the GTN model, confirming the potential of the former constitutive model to fulfill to the requirement of transferability between different loading conditions. • Noticeable disagreement has been observed between predictions of both models at and beyond failure points of specimen. The necking round bar problem has shown that during the sudden specimen capacity loss, the present model based on SC microstructure is the closer one to the GTN model. REFERENCES [1] V. Tvergaard, “Materials failure by void growth to coalescence”, in: Advances in Applied Mechanics, Vol. 27, pp. 83–151, Academic Press, New York, 1994 [2] A. A. Benzerga, J. -B. Leblond, “Ductile fracture by void growth to coalescence”, in: Advances in Applied Mechanics, Vol. 44, pp. 169– 305, Academic Press, New York, 2010 [3] D. L. S. McElwain, A. P. Roberts, A. H. Wilkins, “Yield criterion of porous materials subjected to complex stress states”, Acta Materialia, Vol. 54, No. 8, pp. 1995-2002, 2006 [4] L. Siad, “Ductile damage analysis based on a J3-GTN-like model”, Journal of Multiscale Modelling, Vol. 3, pp. 1-27, 2011 [5] V. Tvergaard, A. Needleman, “Analysis of the cup-cone fracture in a round tensile bar”, Acta Metallurgica, Vol. 32, No. 1, pp.157–169, 1984 [6] N. Aravas, “On the numerical integration of a class of pressure- dependent plasticity models”, International Journal for Numerical Methods in Engineering, Vol. 24, No. 7, pp.1395–1416, 1987 [7] J. Simo, .J.R. Hughes, Computational inelasticity, Springer, 1998 [8] A. Anandarajah, Computational Methods in Elasticity and Plasticity. Solids and Porous Media, Springer, 2010 [9] J. Besson, ‘Eprouvettes axisymétriques entaillées’’, in D. François, éditeur, Essais mécaniques et lois de comportement, pp. 319–351. Hermès Science, Paris, 2001 [10] N. Bonora, “Identification and measurements of ductile damage parameters”, The Journal of Strain Analysis for Engineering Design, Vol. 34, No. 6, pp. 463–478, 1999