Ductile fracture simulations using a multi-surface coupled damage-plasticity model

Ductile fracture simulations using a multi-surface coupled damage-plasticity model XIV International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XIV E. O˜nate, D.R.J. Owen, D. Peric & M. Chiumenti (Eds) DUCTILE FRACTURE SIMULATIONS USING A MULTI-SURFACE COUPLED DAMAGE-PLASTICITY MODEL D. REDDI†AND S. M. KERALAVARMA Department of Aerospace Engineering Indian Institute of Technology Madras, Chennai 600036, India †e-mail: [email protected] Key words: Ductile Fracture, Void Growth, Void Coalescence, Shear Failure Abstract. In this paper, an isotropic porous metal plasticity model accounting for both void growth by diffuse plastic deformation and void ‘coalescence’ by localization of plastic flow in the inter-void ligaments is presented. Predictions for the effective stress-strain response, evolution of damage and the strains to failure are obtained by integrating the model numerically under triaxial proportional loading conditions. The model predictions are compared with results from micromechanical finite element simulations of the average response of voided unit cells under similar loading conditions. It is shown that the model predictions for the failure strains as a function of the loading path are in good qualitative agreement with the results of the cell model simulations. 1 INTRODUCTION Fracture of ductile materials is usually preceded by the localization of plasticity in a failure process zone such as a diffuse neck or a shear band. Material separation occurs by the growth and coalescence of micro-voids inside the process zone that initiate from second phase particles or inclusions [1]. Both the condition for the onset of plastic instabilities and the rate of crack growth within the localization zone are strongly influenced by the local state of stress. Void growth by diffuse plastic flow around the voids depends on the relative magnitude of the hydrostatic stress, while void coalescence occurs due to plastic collapse of the ligament separating neighboring voids, which depends on the ligament thickness and stress components in the transverse plane of the ligament. Denoting the stress state via two commonly used non-dimensional parameters, the stress triaxiality, T, defined as the ratio of the mean and Von Mises effective stresses and the Lode parameter, L, proportional to the determinant of the deviatoric stress, void growth by diffuse plastic flow depends only on T, while void coalescence in general depends on both Tand L. Predictive modeling of ductile fracture therefore requires a physics-based plasticity model that not only accounts for the effect of the loading path on the damage growth rates, but also the condition for the onset of void coalescence and consequent rapid softening due to transitions in the deformation mode at the micro-scale. 1 534 D. Reddi and S. M. Keralavarma Micromechanics-based plasticity models have been developed since the late 60’s to predict the evolution of plastic strains and damage in a porous ductile material in the pre-coalescence phase [2, 3]. The Gurson model [3] in particular has been successfully used to predict several experimental features of ductile fracture [4–6] using a heuristic criterion for the onset of coalescence proposed by Tvergaard and Needleman [4]. In the Gurson-Tvergaard-Needleman (GTN) model, void coalescence is assumed to initiate once the porosity reaches a critical value. When the critical porosity is exceeded, the damage growth rate is accelerated to simulate the rapid material degradation as observed in finite element simulations of void growth using the unit cell model. However, a major limitation of the above phenomenological approach is that both the critical porosity for the onset of coalescence and the form of the post-coalescence damage evolution law are calibrated based on cell model simulations using a limited set of loading paths; mainly axisymmetric loadings. Importantly, predictions for the ductility obtained using the GTN model under proportional loading conditions depend only on the stress triaxiality Tand are independent of the Lode parameter L. While the major influence of the stress triaxiality parameter on the ductility has been understood for a long time, both experimentally [7] and computationally [8], the importance of the Lode parameter has only been appreciated recently, with the publication of experimental data that appears to suggest that the ductility under low triaxiality loading (typically T<1) depends on the Lode parameter with significantly lower ductilities predicted under shear dominated loadings compared to axisymmetric loadings at higher triaxialities [9, 10]. Subsequently, several authors have investigated further the effect of the Lode parameter on the strains to the onset of void coalescence at the micro-scale using three-dimensional cell model simulations of void growth under combined axisymmetric and shear loading [11–13]. A significant reduction in ductility has been consistently observed in these simulations under shear dominated loadings. These results clearly establish that the critical porosity criterion in the GTN model needs to be replaced with a stress-based criterion that includes the observed effect of the Lode parameter Lon the onset of coalescence, since the value of Lallows to distinguish between axisymmetric and shear dominated loadings at the same triaxiality. In the past few years, several authors have attempted to develop such models by extending the Thomason [14] model for void coalescence by internal necking to account for coalescence under combined tension and shear [15–17]. Most recently, Keralavarma [18] proposed a multi-surface plasticity model for void growth and coalescence in an isotropic material by combining the Gurson [3] void growth model with the void coalescence model of Keralavarma and Chockalingam [17], appropriately extended to account for arbitrary orientations of the coalescence band. The objective of this paper is to examine predictions for the material’s intrinsic ductility as a function of the applied loading path predicted by the above multi-surface model under proportional loading conditions, and to compare the predicted trends with those observed in recent cell model simulations under combined tension and shear [11–13]. A brief summary of the model is presented in section 2, followed by comparison of the model predictions with two-dimensional axisymmetric cell model simulations and predictions from the GTN model in section 3.1. Finally, predictions for the variation of the material’s ductility under general proportional loadings as a function 2 535 D. Reddi and S. M. Keralavarma of the loading path parameters Tand Lare presented and discussed in section 3.2. 2 MULTI-SURFACE POROUS PLASTICITY MODEL Consider an elasto-plastic Von Mises material containing a random distribution of equiaxed voids at the micro-scale. Assuming a dilute volume fraction of voids with porosity f1, the yield stress of a macroscopic material element may be predicted using the Gurson [3] yield criterion FG(Σ)≡Σ2 eq ¯σ2+2qf cosh 3 2 Σm ¯σ−1−(qf)2= 0 (1) where Σis the macroscopic stress tensor, Σm=1 3tr(Σ) and Σeq =3 2Σ:Σare the mean and Von Mises equivalent stresses respectively, Σis the deviatoric stress, ¯σis the flow stress of the matrix material and qis a heuristic parameter introduced by Tvergaard [19]. Capital symbol Σis used for the stress to emphasize the fact that Σis the average stress tensor over a porous representative volume element (RVE) Ω; i.e. σ=σΩwhere σis the Cauchy stress. The above yield function can be formally derived using limit analysis by assuming that plastic flow occurs in the entire RVE during yielding. However, at finite values of the porosity, an RVE can also yield by localized plastic flow (coalescence) in a narrow band encompassing the ligaments connecting neighboring voids and band width equal to the void diameter. For a given orientation of the localization (or coalescence) band identified by the unit normal vector n, Keralavarma and Chockalingam [17] derived the following yield function FC(Σ, n)≡3Σ2 sh ¯σ2+2fbcosh Σn Σc−1−f2 b= 0 (2) where Σn=n·Σnand Σsh =n·Σ2n−Σ2 nare respectively the normal and shear stresses on the coalescence plane, fb=f2/3and Σcis a critical stress given by Σc=¯σ6 5√b2+1−b2+f2 b+bln b+b2+f2 b fb(b+√b2+ 1)log 1 fb−1 (3) with the parameter bgiven by b=1 3+5 288 1+ 1 fb−5fb+3f2 b(4) Σcis a positive definite function of fbthat tends to infinity in the limit fb→0+, so that at dilute porosities, the yield surface defined by (2) falls outside the Gurson yield surface (1) for most loading paths. However for finite values of f(typically f>∼0.01), the yield stress predicted by the coalescence criterion can be lower than the Gurson value for several loading paths in stress space. In a statistically isotropic material, the inter-void ligament dimensions are approximately the same along any material direction. Therefore, according to limit analysis 3 536 D. Reddi and S. M. Keralavarma theory, coalescence occurs along the direction nthat yields the smallest value of the yield stress or, equivalently, maximizes the coalescence yield function FCin (2). Thus, the coalescence criterion for an isotropic material is written as FCiso (Σ)≡max nFC(Σ, n) = 0 (5) Performing the above maximization leads to the result that coalescence occurs either on one of the principal stress planes or on a non-principal plane non which the normal stress Σnsatisfies the following equation (see [18]) 2Σn−2 3 ¯σ2 Σcfbsinh Σn Σc=Σ 1+Σ 2(6) where Σ1and Σ2are two unequal principal stresses. Solving the above equation for every unequal pair of principal stresses yields the normal stresses on planes where shear assisted coalescence can occur, and the corresponding shear stress is found from the equation Σ2 sh =Σ1−Σ2 22 −Σn−Σ1+Σ 2 22 (7) The value of the isotropic coalescence function FCiso then corresponds to the maximum value of FCover all (Σm,Σsh) pairs obtained above, and coalescence occurs when FCiso ≥ 0. Finally, combining the Gurson and the isotropic coalescence models using the same multi-surface approach, Keralavarma [18] proposed the following for the effective yield criterion for a porous isotropic material accounting for both void growth and coalescence F(Σ)≡max FG(Σ),FCiso (Σ)= 0 (8) The macroscopic plastic strain rate, Dp, is obtained from the yield function via the normality property, which yields Dp=˙ λN,N=∂F ∂Σ(9) where ˙ λis the plastic multiplier and the direction of plastic flow Ndepends on the active yield surface, F=FGor FC. We have N=         3Σ ¯σ2+qf ¯σsinh 3 2 Σm ¯σI,F=FG 3 ¯σ2[n⊗Σn+Σn⊗n−2Σnn⊗n]+2fb Σcsinh Σn Σcn⊗n, F=FC(n) (10) The evolution of porosity follows from plastic incompressibility of the matrix, which yields ˙ f 1−f= tr(Dp)=       ˙ λ3f ¯σsinh 3 2 Σm ¯σ,F=FG ˙ λ2fb Σcsinh Σn Σc,F=FC(n) (11) 4 537 D. Reddi and S. M. Keralavarma Strain hardening in the matrix is accounted for by making the matrix yield stress a function of the plastic strain as ¯σ=¯σ(p eq), where p eq is an average measure of the equivalent plastic strain in the matrix material, whose evolution is obtained from the equivalence of the plastic power at the macroand micro-scales; i.e. Σ:Dp= (1 −f)¯σ(p eq)˙ p eq (12) In this paper, a power-law relationship between p eq and ¯σis assumed, of the form ¯σ=σ01+p eq 0n (13) where σ0is the initial yield stress, 0is a reference strain and nis the strain hardening exponent. In the following section, rate equations (9)–(12) are integrated along radial loading paths in stress space characterized by constant values of the loading path parameters, the triaxiality Tand the Lode parameter L, to obtain the evolution of the equivalent stress and porosity as a function of the equivalent plastic strain p eq.Tand Lare related to the invariants of the stress tensor as T=Σm Σeq ,L=−27 2 det(Σ) Σ3 eq (14) Tis a measure of the magnitude of the hydrostatic stress relative to the deviatoric stresses, while the value of Lallows to distinguish between different states of stress for the same T.Lis bounded between -1 and +1, with L=−1 for axisymmetric loadings with a major axial stress, L= 0 for pure shear with superposed hydrostatic stress and L= +1 for axisymmetric loading with major radial stresses. 3 RESULTS 3.1 Comparison with axisymmetric cell model simulations The heuristic parameters in the classical GTN model, namely the Tvergaard parameter qand the critical porosity for the onset of coalescence, fc, are usually calibrated by comparison with cell model simulations under axisymmetric loading conditions (e.g. see [8,19]). In the latter, the average response of a transversely isotropic distribution of voids in an elastic-plastic matrix is simulated using finite element analysis of a two-dimensional RVE subjected to proportional axisymmetric loading, as shown in Fig.1. A periodic distribution of voids in the plane of analysis is assumed so that, exploiting the symmetries of the geometry and the applied loading, the average response of the unit cell shown using dashed lines in Fig.1(a) can be obtained from analysis of the quarter cell shown in Fig.1(b), subject to symmetry boundary conditions on the inner boundaries (edges that intersect the voids) and periodicity condition on the outer boundaries. The unit cell is loaded in such a way that the ratio of the principal stresses in the axial and radial directions, or equivalently the triaxiality T, remains constant during the deformation. The major stress is applied in the axial direction, so that the Lode parameter L= − 1 in all the simulations 5 538 D. Reddi and S. M. Keralavarma z r Σzz Σrr (a) (b) Figure 1: (a) Schematic of a transversely isotropic porous material containing a periodic distribution of voids in the plane of analysis, (b) one quarter of a periodic unit cell used in the finite element analysis. presented here. Further details of the axisymmetric cell model simulations can be found in [8]. Fig.2 shows the response of periodic unit cells of the type shown in Fig.1(b) made of an elasto-plastic Von Mises material with power law hardening. The void shape is assumed to be initially spherical with a volume fraction f=0.001. The values of the material properties assumed are Young’s modulus E= 210 GPa, Poisson ratio ν=0.3, initial yield stress σ0= 420 MPa, hardening exponent n=0.1 and reference strain 0=0.002; see Eq.(13). Fig.2(a) shows the equivalent stress-strain response of the unit cell (solid black lines) for axisymmetric loadings with L=−1 and four different values of the triaxiality T=2/3,1,3/2 and 2. The equivalent stress and strains are defined as Σeq =|Σzz −Σrr|,E eq =2 3|Ezz −Err|(15) where Σ=σΩis the average Cauchy stress and Err and Ezz are the logarithmic strains in the radial and axial directions of the cell respectively. The corresponding evolution of the porosity fas a function of equivalent strain is shown in Fig.2(b). The simulations are terminated at the onset of coalescence (indicated by the ×symbol), when the deformation localizes into the transverse ligament between the voids, which also coincides with the unit cell switching to a uniaxial mode of deformation along the zdirection. The figure also shows predictions of the effective stress-strain response and damage growth obtained from the the multi-surface plasticity model summarized in the previous section (blue dashed lines) and the GTN model (red dotted lines). In the pre-coalescence regime, the GTN model is identical to the multi-surface model, since the yield surface coincides with the Gurson yield surface, i.e. F=FG, and the state evolution equations are identical. However, the criterion for the onset of coalescence is different for the two models. In the GTN model, coalescence occurs when the porosity reaches a critical value fc, while coalescence occurs in the multi-surface model when FG≤FCiso = 0; see section 2. In both 6 539 D. Reddi and S. M. Keralavarma 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Σeq/σ0 Eeq T = 2 T = 3/2 T = 1 T = 2/3 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 f Eeq T = 2 T = 3/2 T = 1 T = 2/3 (a) (b) Figure 2: Comparison of the effective response of two dimensional porous unit cells shown in Fig.1 subjected to axisymmetric proportional loading with L=−1 and various values of T(solid lines). Predictions from the GTN (dotted lines) and multi-surface (dashed lines) porous plasticity models are also shown. (a) Equivalent stress vs. strain and (b) porosity vs. equivalent strain. cases, the simulations are terminated when the onset of coalescence is detected. The values of the GTN model parameters are adopted from [8] with q=1.25 and fc=0.03. The same value of the Tvergaard parameter qis also used in the simulations using the multi-surface model. Further, the coalescence yield function of Eq.(2) is heuristically modified in the spirit of Tvergaard’s modification of the original Gurson model to redefine the effective porosity parameter as fb=qbf2/3, where qb=q2/3is assumed in the present calculations. Comparison of the model predictions with the cell model simulations in Fig.2 shows that both the GTN and the multi-surface models provide reasonable predictions for the strains to the onset of coalescence. However, neither model is satisfactory for predicting the correct trends for the porosity at the onset of coalescence, which appears to increase with Tin the cell model simulations, while the multi-surface model predicts the opposite trend. Nevertheless, it is significant that the multi-surface model can capture the correct trends for the ductility in the axisymmetric simulations, without use of the heuristic fc parameter, because the coalescence criterion in the multi-surface model depends on both the triaxiality Tand the Lode parameter L, unlike the GTN model, which depends only on T. Hence, the stress state dependence of the ductility can now be examined under general triaxial loading conditions, as is done in the next section, and compared with the trends reported in the recent literature. 3.2 Loading path dependence of the ductility Fig.3 shows the effective stress strain response predicted by the multi-surface model under triaxial radial loading with T= 1 and several values of the Lode parameter L. Unlike in Fig.2, the post-coalescence response of the material is also included in Fig.3. Notice that the effective stress-strain response depends on the Lode parameter unlike the 7 540 D. Reddi and S. M. Keralavarma 0.6 0.8 1 1.2 1.4 1.6 1.8 0 0.2 0.4 0.6 0.8 1 1.2 Σeq/σ0 εp eq L = -1 L = -0.75 L = -0.5 L = 0 L = 0.5 L = 0.75 L = 1 Figure 3: Equivalent stress vs. strain response predicted by the multi-surface plasticity model under triaxial proportional stressing with T= 1 and various values of L. GTN model, whose predictions essentially collapse into a single curve (not shown). The onset of coalescence is visible as a sharp change in the slope of the stress-strain curve in some of the simulations for large |L|. However, for shear dominated loadings (|L|near zero), coalescence tends to occur early and without a sharp change is the slope of the stress-strain curve. Also, the material undergoes significant stable plastic flow and strain hardening after the onset of coalescence, so that the strain to the onset of coalescence is not an accurate measure of the ‘intrinsic’ ductility of the material. It is more reasonable in such cases to adopt the plastic strain corresponding to the maximum in the stressstrain curve as the material’s ductility. For the present study, we thus adopt the following measure of the ‘failure strain’ p fof the material p f= max {p ult, p coal}(16) where p ult and p coal are the equivalent plastic strains corresponding to the ultimate stress and the onset of coalescence respectively. The results in Fig.3 show that the strain to failure under proportional stressing p f, as defined above, predicted by the multi-surface model shows a non-monotonic trend with respect to Lat fixed T, and the minimum ductility is predicted for shear dominated stress states near L= 0. The Lode parameter dependence of the ductility is further illustrated in Fig.4(a). The figure plots the strain to failure p fas a function of Lfor several representative values of T. Notice that the predicted curves have an approximately convex shape, except for a region near L= 0 where slight concavity is observed. The ductility minimum occurs for shear dominated stress states with small negative values of L. Further, p falso shows a dependence on the sign of L(equivalently the sign of the determinant of the deviatoric stress), with relatively higher ductilities predicted for 8 541 D. Reddi and S. M. Keralavarma 0 0.5 1 1.5 2 2.5 3 -1 -0.5 0 0.5 1 εp f L T = 0.5 T = 0.75 T = 1.0 T = 1.5 T = 2.0 0 0.5 1 1.5 2 2.5 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 εp f T L = -1 L = 0 L = 1 (a) (b) Figure 4: Variation of the strain to failure p funder proportional loading as a function of: (a) the Lode parameter Lfor various values of Tand (b) the triaxiality Tfor various values of L. positive values of L. In contrast, the triaxiality dependence of the ductility in Fig.4(a) exhibits a monotonic trend as expected, with lower failure strains predicted towards higher values of T. The triaxiality dependence of the ductility is also illustrated in Fig.4(b), which plots p fas a function of Tfor several values of L. The shapes of these curves are in accordance with predictions from classical porous plasticity models, except for the significant dependence on L. Recently, several authors have performed cell model simulations in the spirit of section 3.1, using three dimensional unit cells subjected to periodic boundary conditions and combined tensile and shear loads to simulate proportional stressing for arbitrary values of Tand L[11–13]. They report predictions for the strains to the onset of coalescence in remarkable qualitative agreement with the predictions from the multi-surface model in Fig.4. In particular, the shape of the p fvs. Land p fvs. Tcurves and the dependence of the failure strains on the sign of Lare in qualitative agreement with the above cell model simulations, which indicates that the mechanisms of coalescence assumed in the multi-surface model are fundamentally correct. It remains to perform a quantitative comparison and calibration of the model against three dimensional unit cell simulations, with possibly introduction of additional heuristics to correct for some of the discrepancies observed in Fig.2. The results of such a study will be reported in a future publication. 4 CONCLUSION Conclusions from the above study are summarized below. -It is shown that a multi-surface porous plasticity model [18], combining the Gurson model with a void coalescence model accounting for arbitrary orientations of the coalescence band and the effect of the loading path on the coalescence stress, can predict the triaxiality and Lode parameter dependence of the ductility observed in 9 542