INTERNATIONALJOURNALFOR NUMERICALMETHODS IN ENGINEERING, VOL. 39,619-633 (1996) FINITE ELEMENT METHOD FOR GRADIENT PLASTICITY AT LARGE STRAINS XIKUI LI AND S. CESCOTTO' National Laborarory for Structural Analysis of Industrial Equipment, Dalian University of Technology, Dalian 116023, P.R. China SUMMARY A finite element method for gradient elasto-plastic continuum in which the yield strength of strain hardening/softeningmaterials not only depends on the effective plastic strain but also on its Laplacian is presented. The consistent integration algorithm to update the stress and the internal state variable at integration points and the consistent compliance matrix for the gradient plasticity are formulated in the non-local sense. The methodology to derive the finite element formulation for the gradient plasticity at large strains presented in this paper is applicable to general finite element analysis; the formulation in the context of the two-dimensional four-noded mixed finite element with one integration point and mean von Mises yield criterion is particularly derived. Numerical examples are tested to demonstrate the capability and performance of the present finite element method at large strain in solving for the strain localization problem. KEY WORDS gradient plasticity; strain localization; mixed strain element method; large strain; consistent algorithm 1. INTRODUCTION Localization of deformation into narrow bands of intense straining caused by strain softening is a characteristic feature of plastic deformation. It has been experimentally observed in many engineering materials such as metals, polymers, concrete and geological materials. Localization phenomena are often associated with a significant reduction of the load-carrying capacity of the structures; hence, the onset of localization is naturally considered as the point of threshold to predict rupture of engineering structures. Over the last decade considerable effort has been devoted to obtain a comprehensive understanding of the problem and to describe the behaviour quantitatively. Numerous attempts to simulate the behaviour using classical plastic continuum theories have been unsatisfactory.The principal defect of the classical plastic continuum models is that they do not incorporate an internal length scale or higher-order continuum structure and, therefore, the governing equations under quasi-static loading lose ellipticity. Consequently, as strain softening behaviour is incorporated into a computational model the ill-posedness of the problem will lead to pathologically mesh-dependent results. Furthermore, the energy dissipated at strain softening is incorrectly predicted to be zero, and the finite element solutions converge to incorrect, physically meaningless ones as the element mesh is refined.' *Permanent address: M.S.M.Department, University of Liege, Quai Banning 6,4000,Liege, Belgium CCC 0029-5981 /96/040619- 15 0 1996 by John Wiley & Sons, Ltd. Received 22 August 1994 Revised 29 January 1995 620 X. LI AND S. CESCOTTO To remedy the situation generalized continuum theories, in which internal length scales or high-order continuum structures are introduced as regularization mechanisms for loss of ellipticity, should be adopted. Among those are gradient continuum theories,'-' Cosserat continuum theory" and rate-dependent continua.'' In the present work we will utilize the gradient continuum approach. Many efforts have been devoted to developing the models using the gradient continuum theories. Bazant presented the non-local gradient formulation with a Laplacian of strain2 and applied the formulationin solving a one-dimensionallocalization p r ~ b l e mAifantis4* .~ show that the introduction of the second gradient of deformation into the expression of the flow stress for plastic materials can preserve ellipticity of the governing equations. Lasry and Beiytschko6 introduced additional higher-order terms in the strain displacement relation for the strain softeningportions of the domains for the proposed non-local formulations. Muhlhaus and Aifantis' proposed their non-local gradient model on the basis of a variational principle for gradient plasticity, in which higher-order spatial gradients of the effective plastic strain are included in the yield condition. Recently De Borst and his c o - w o r k e r ~presented ~*~ the formulations and algorithms for the gradient plastic continuum in a finite element context and in these the yield strength depends not only on the effective plastic strain but also on its Laplacian. It makes the presented formulations and algorithms able to solve two- or three-dimensional boundary value problems. As the Laplacian of the effective plastic strain is involved in the yield function and the increment of the plastic multiplier cannot be obtained at a local (integration point) level, in contrast to the conventional approach in computational plasticity, the plastic multiplier is taken as an independent global variable. The consistency condition of the yield function is enforced only in a weak form; in addition, it is not satisfied at each iteration but only at the end of a load step. With the use of the gradient plastic model including Laplacian of the effective plastic strain in the yield condition, a finite element method for gradient elasto-plastic continuum is presented in this work. The consistent integration algorithm to update the stress and the internal state variable at integration points and the consistent compliance matrix for the gradient plasticity are formulated in the non-local sense. Even though the yield condition at each integration point is no longer local, it is satisfied in a pointwise fashion and at each iteration of a load step. The plastic multipliers are internal state variables locally defined at each integration point. The Laplacian of the effective plastic strain at an integration point is evaluated from the derivatives of the polynomial, which is interpolated by using the values of the plastic multipliers at neighbouring integration points. As the plastic multipliers are not taken as independent global unknowns there is no requirement on C'-continuous shape functions for the interpolation of the plastic multiplier. The methodology to derive the finite element formulations for the gradient plasticity at large strains presented in this paper is applicable to general finite element analysis. Even the formulations in this paper are particularly derived in the context of the two-dimensional four-noded mixed finite element with one integration point and von Mises yield criterion," the methodology can be applied to general cases with the use of displacement-based finite elements, more general yield conditions and non-local softening law such as that for quasi-brittle materials exhibiting microcracking and internal friction. 2. CONSTITUTIVE EQUATIONS OF THE GRADIENT DEPENDENT VON MISES PLASTICITY MODEL We start with the rate constitutive equation for elasto-plastic materials: GRADIENT PLASTICITY AT LARGE STRAINS 621 where 8, %, 8, are total, elastic and plastic strain rates, respectively. D, is the matrix of elastic moduli. The rate of plastic strain vector can be expressed as tp where the rate of plastic multiplier associated plasticity is defined by = in (2) is a non-negative scalar and the flow vector n for the n = af/aa (3) Here the yield functionf for the gradient plasticity is non-local and can be assumed in the form f =f(a, E P , V%P) (4) The assumed effectiveplastic strain Ep can be calculated by the integration of incremental effective plastic strain AEP: where AEp is defined by As the von Mises yield criterion is considered, the flow vector n takes the form 3 n=-s (7) 28 with s being the deviatoric stress vector and 6 the effective stress defined by 5 = (3STS)l’Z (8) The non-local von Mises yield surface can simply be written as f (a,E P , VZEP) = @(a) - ay(EP,VZEP) = 0 (9) Here ay is the current yield stress and is assumed as piecewise linear function of EP as well as a linear function of the Laplacian of the effective plastic strain: ay = uyo + h,(EP)P - cV%P (10) where ayois the initial yield strength, h, is the hardening/softening modulus and is a function of Ep, c is a non-local material parameter defined by c = af/avZEp According to equations (2) and (6)-(8) it is given that . . 1=g p The variation of the yield function with time can then be written as f = 8 - h,X +cVZi (11) 622 X. 11 AND S. CESCOTTO Consider a load step from time t to time t + At. Integration of equation (1) with the time gives - AADent+Ar Gf+At =a;+~t where the trial elastic stress at time t (14) + At is defined by GF+At = De(&r+Af- The flow vector at the point eE in the stress space is where sEis the deviatoric component of 8 and 5 E = (3SETSE)1/2 It can be verified for plane strain, axisymmetric and three-dimensional solid that E nf+At = nt+Ar (18) Pre-multiplying by nT both sides of equation (14) and using equation (18) gives - 4 + a t = 8,"+~t KEAJ (19) where the effective elastic stiffness (Gis the shear modulus) is given by KE= nTD,n = 3G (20) Substitution of the variation 8 evaluated from equation (19) into equation (13) results in f= - (KE + hp)X + cV2A (21) It is observed that, unlike in classical (local) computational plasticity, incremental plastic multiplier cannot be determined at a local (integration point) level. 3. INTEGRATION ALGORITHMS FOR STRESS UPDATES AND CONSISTENT TANGENT MATRIX Consider a typical time step [t, t + At] for an arbitrary material point i in plastic yielding. The yield condition at the current iteration k can be written in an iterative manner as 5 . k =fi,k- 1 f SA = A , k - 1 - ( K E+ hp)id(AAi) + CV26(AAi) = O (22) where As shown in equation (22) the yield condition at an integration point is related to not only the local plastic variable but also to that in the vicinity of the integration point. To evaluate V2 A&, we consider the neighbouring points around the integration point i. Let N i be the number of the neighbouring integration points, including the point i itself. As the plastic multipliers at these N iintegration points are used to compute V2AJi, we are able to express GRADIENT PLASTICITY AT LARGE STRAINS 623 V2 A l i in terms of Adj ( jE N i ) in the form To determine the coefficients gij, a complete second-order polynomial in two co-ordinates is assumed to represent the function of plastic multiplier around point i, i.e. A 1 = aTv (25) where a = [al a2 a3 a4 a5 aJ (26) v = [ 1 x y x2 x y y2-J‘ (27) The coefficient vector a is obtained by minimization of N, Error = (AA - Ahj)’ (28) j= 1 with AA = [A& . . . A&JT denoting the values of Al, at the neighbouring integration points of the integration point i, and we get the result of vector a in the form Ga = M A A where G= C in equation (30) stands for Cy:,. Matrix M is defined by M = [vI ~2 . . . vj . . . vNi] with vj being the vector v at jth neighbouring integration point. Hence, the coefficientsa can be expressed in terms of the values of the AA of the neighbouring integration points Ni a = G - l M A A = G-l j= 1 Aljvj Though V2 AA can be evaluated, with equations (25), (26) and (32), as V2 AI = 2(a4 + aa) (33) 624 X. Ll AND S. CESCOITO to express V2A& explicitly in terms of AAj (j E Ni) in the form (24)we should determine coefficient gij given by gij = g'vj (34) gT = 2(4th row of G-' + 6th row of G-') (35) where It is remarked that the coefficients gij only depend on the co-ordinates of the neighbouring integration points. In case of small displacements, they are computed once for all and used at each time step to evaluate the Laplacian by (24). However, in case of large strains, they should be updated when geometry changes are significant. The choice of the neighbouring integration points is simply based on the mesh topology. Generally, the integration points of the neighbouring elements are used. However, if these points are not sufficient to allow the computation of G(i.e. if G is singular or ill-conditioned),indirect neighbours (i.e. the neighbours of the neighbours) must be considered. It is worthwhile to point out that the polynomial approximation assumed for the plastic multiplier function in the neighbourhood of the integration point i is used only as a means to determine approximately the Laplacian of the plastic multiplier at the point, and does not mean that the plastic multiplier itself is approximated in a least-squares sense. Substitution of expression (24) into equation (22) gives ' for each integration point i (i = from 1 o N,, N , is the total number of the integration points in the finite element mesh). To ensure quadratic convergence of the global Newton iterative procedure, a consistent tangent matrix for the gradient-dependent von Mises elasto-plasticity has to be formulated. For classical continuum models we can write for each integration point i (i = 1 to Nt),where iri and iiare the stress and the strain rates at the local point i, Dep,is the local consistent tangent matrix. It is noted that the consistent tangent matrix at an integration point for the gradient plasticity is no longer local, the stress rate 6i at point i depends not only on the internal state variable at the local point i, but also on those at the neighbouring integration points. Consequently, in contrast to the form (37), the consistent tangent constitutive relation of the point i will be given in the form N. j= 1 where N, stands for the number of integration points within the elasto-plastic zone, in which point i lies. However, it can be numerically observed that only those Depvrelating to integration pointj, which are among the neighbouring ones of the point i, give the essential contributions to the constitutive relation (38). Hence, it is accurate enough to assume that N. GRADIENT PLASTICITY AT LARGE STRAINS 625 Based on this fundamental consideration, let us consider the non-local consistent tangent matrices Dep, (j = 1,N i ) at integration point z, with the use of the mixed four-node element.” Let L(i) and N ( i ) be the set of the neighbouring integration points and the number of these points for integration point i. Define Li(k) = L ( i ) n L ( k )and let Ni(k) be the number of points in the intersection Li(k). The variation of equation (19),together with the consistency condition f k = 0 and equation (13), for each point k E L(i) at time t -k At gives j = l,jeL.i(k) Solution of the equation set (40) gwes the rates of the effective plastic strains for each point k E L(i) N (i) g:.f+At = Pkj$Er+Af (41) j= 1 Substitution of equation (41) into the rate form of the equation (19) for point i gives where In the strain driven finite element analysis the rate of the trial elastic stress +Af is entirely determined by the strain rate. A constitutive relation between hEand the strain rate can be built up. To do this, we first review the mixed finite element12 used in the present work briefly. The strain of the mixed element is chosen in the form E =f + h,$ (44) and correspondingly, the stress takes the form u = 6 + h*,UX (45) where a overbar denotes a constant field, the second term on the right-hand side in each of equations (44)and (45) varies with co-ordinates. As the strain vector for the plane strain analysis with the form E = [ex. eyy E,, cxyITis assumed, we have with h=fA(q where A is the area of the finite element. 626 X LI AND S CESCOTT'O The constant fields (6,E) are uncoupled with the anti-hourglass mode (a",sX).The stress m and the strain E can be further decomposed into deviatoric and volumetric parts by where s =S + h,,a' (49) and s=e+d" where + h,.Er = E - = f-T'f e =8 &" @ (53) withm= [l 1 1 0IT. Since the anti-hourglass mode is introduced into the element, one integration point scheme is used for two-dimensional four-noded element while neither zero-energy modes nor locking phenomena for shear or incompressible materials arises. The effective stress defined in equation (8) is replaced by the mean value of it over an element as 1 A cF2 = - / A ($sTs)dA (W Correspondingly,CE redefined as the mean value of it over the element and its derivative can be given by 1 (cEIZ = - JA A (tsE'sE)dA where Define the ratio The rate of the deviatoric stress can be divided into a local and a non-local part as si = sf:+ 5;' (58) GRADIENT PLASTICITY AT LARGE STRAINS 621 where After some classical deductions the integration of the virtual work done by the stress rate 6 over the element i gives The first and the second terms on the right-hand side of the equation represent the local and non-local tangent stiffness, with the local consistent compliance matrices C k i = A&P - B3SiST) CieSi= - 2jI,~,(af)~H~ C:evi= - 2P3HiaT$ and the non-local consistent compliance matrices The coefficients in equations (62) can be given by and the matrix P for the plane strain and the axisymmetric solid problems is given by 3 -3 1 -4 0 2 (65) 0 2 To analyse geometrically non-linear problems with finite strains and/or finite rotations by utilizing the mixed element, a co-rotational formulation is adopted. For each element a local 628 X. LI AND S.CESCO’M’O co-ordinate system is embedded at the centre (the integration point) of the element and rotates with it. The element displacements,strains and stresses are evaluated with reference to the local co-ordinates at the current time.” Consider a load step from time t to t At. The algorithm described in the proceeding equations can be summarized as follows: + (1) Compute the element stiffness matrices according to the local and the non-local consistent tangent matrices given by (62) and (63) for each element. (2) Solve for the global displacement increments and then obtain the total displacements at time t + At referred to the global co-ordinates. (3) Determine the orientation of local axes at time t + At for each element according to the total displacements at time t At, and the element displacementsat the last converged load step referred to the local co-ordinate system of the element at time t. The incremental strains at time t + At referred to the local axes at time t + At is then obtained.12 (4) Check the yield function for each integration point according to equation (23), i.e. check if the point is in elasticity, or on the current yield surface or outside the current yield surface. Then 6(AAi)for each integration point i is calculated non-locally according to equation (36). Since this equation is linear in the unknowns and is diagonal dominated, Gauss-Seidel iterations are used at this level. ( 5 ) Accumulate the effective plastic strain, its Laplacian and update the stress vector for each integration point. (6) Check global convergence criterion. If convergence does not occur go to (1). + 4. NUMERICAL EXAMPLES We consider a rectangular bar with initial length Lo = 5 and initial width Bo, meshed by N x M four-noded rectangular mixed elements, shown in Figure l(a). The bar fixed at the bottom and loaded by an increasing prescribed horizontal displacement at the top, is analysed as a plane strain problem. In order to have a pure shear problem, the vertical displacements of all the nodes are prevented and the horizontal displacements of nodes pertaining to the same horizontal line are imposed to be equal. It is noted that at least N = 3 elements over the width Bo of the bar are required. Indeed, with N = 1 or 2 and only one integration point per mixed element, it would be impossible to approximate the plastic multiplier by (25). To demonstrate that pathological mesh dependence has been overcome by using the present finite element method for gradient plasticity, we consider first an example. Four test cases with different dimensions in x co-ordinate and element meshes are specified Case 1: Bo = 3-0, M = 5, N = 3; Case 2 Bo = 1-75, M = 15, N = 5; Case 3: Bo = 1-0,M = 25, N = 5; Case 4 Bo = 0.7, M = 35, N = 5. The values of material elastic properties used for the entire bar are E = 210000 and p = 0.0, The initial yield strength ayo= 240 except for the elements which are between y = - 05 and y = 0.5. The initial yield strength for these elements has been weakened to the value of cy0= 230. As the effectivestress at a material point reaches oy0,a softening modulus h, = - loo00 is used for all elements, including the weakened elements. A value of the non-local material parameter c = 2000,which is related to the internal length I by8 is used for the four test cases. + GRADIENT PLASTICITY AT LARGE STRAINS I I Lo 629 - - X 1 - %- Figure l(a). A two-dimensional bar with length L and width B subjected to a shear stress Numerical results show that the thickness of the localization zone and the effective plastic strain distribution along y axis over the zone rapidly converge to a unique solution, which is independent of the finite element mesh and is only related to the softening modulus h, and the non-local material parameter c. This is illustrated in Figures 1(b) and 1(c). Figure 1(b) illustrates the convergence of the shear stress-displacement (at the nodes .on the top face) curves. The solution for 3 x 5 element mesh is softer, however, the solutions for the other three cases with 5 x 15, 5 x 25, 5 x 35 elements converge to a unique solution. Particularly, the curves for 5 x 25 and 5 x 35 element meshes almost overlap. Figure 1(c) illustrates the effective plastic strain at the top distributions along y axis for a prescribed horizontal displacement Au = 8-56 x face of the bar. Because of symmetry only the distribution on the half of the bar is illustrated. It is observed that the effective plastic strain profiles for 5 x 15, 5 x 25 and 5 x 35 element meshes almost coincide,though the profile for 3 x 5 element mesh deviates somewhat from the converged solution. It is remarked that the localization zones gradually grow up with the increasing prescribed displacement in the test cases. The thicknesses of the converged localization zones for the latter three test cases agree very well with the analytical solution w determined by8 w = 2 z J X (66) The softer behaviour of the coarse mesh is easily understood considering the strain distribution (Figure l(c)). Indeed, the plastic strain at y = 0 (the integration point of the central element) is much higher than that with the finer meshes: 0.037 instead of 0.026 approximately. From (10)and 630 X 160.00 LI AND S. CJ3CO'lTO 3 Displacement (E-2) Figure 1 (b). Shear stress curves with increasing displacement Au of the top face: (0)3 x 5 elements,(A) 5 x 15 elements, ( x ) 5 x 25 elements, (0)5 x 35 elements I -0.30 Coordinate in Y aixs Figure l(c). Effectiveplastic strain distributionalong Y axis of the bar subjected to a pure shear loading with a weakened part between Y = 0.5 and Y = 0.5 for different discretizations: (0)3 x 5 elements, (A) 5 x 15 elements, ( x ) 5 x 25 elements, (0) 5 x 35 elements - GRADIENT PLASTICITY AT LARGE STRAINS 631 I 0 Q) v, 50.00 Displacement Figure 2(a). Shear stress curves with increasing displacement Au of the top face for different v d u a of the non-local material plastic parameter c: ( X ) c = 450, (A) c = 275, (0) c = 100 with the chosen values of h, and c, the yield stress is much lower for the coarse mesh (despiteof the attenuating effect of the non-local term), which explains the lower shear stress obtained. To show the dependence of the thickness of localization zone on the non-local material parameter c relating to the internal length, we consider a second example. Three test cases, with the same dimensions, element mesh and the initial yield strength as those used for the test case 2 in the first example, are executed. To display the performance of the present method at extremely large strain, we use the material elastic properties E = 100o0,p = 0.0,and the softening modulus h, = - 2500 for the entire bar. The non-local material parameters used for the three cases are: c = 450, c = 275, c = 100, respectively. The shear stress-displacement curves are given in Figure 2(a). Figure 2(b)illustrates the thicknesses of the localization zones and the effectiveplastic strain distributions within the zones for the three test cases. It is observed that using larger value of the internal length (i.e. the non-local material parameter c) results in the wider localization zone and the lower peak value of the effective plastic strain. It can be explained by the fact that more material points in the wider localization zone contribute their residual load-carrying capacity to resist the further plastic deformation due to the gradient (non-local) effect. 5. CLOSURE A finite element method for the gradient plasticity, in which the yield function depends not only on the local internal state variable, but also its spatial gradient, is proposed. The major novel character of the present approach is that the plastic multiplier in plasticity at each integration point is still treated as the internal state variable, determined by the yield criterion at the local point in the non-local sense. Hence the non-local yield condition and the constitutive equation 632 X. LI AND S. CESCOlTO Figure 2(b). Effective plastic strain distribution along Y axis of the bar subjected to a pure shear loadingwith a weakened part between Y = - 0 5 and Y = 0.5 for different values of the non-local material plastic parameter c: ( x ) c = 450, (A) c = 275, (0)c = 100 can be fulfilled in a pointwise fashion at each iteration of a load step, even though as a non-local quantity, the spatial gradient of the local internal state variable is introduced into the yield condition and it is approximately determined by the local plastic multipliers at neighbouring integration points surrounding each individual integration point in a least-squares sense. ACKNOWLEDGEMENTS The authors are pleased to acknowledge the support of this work by the National Science Foundation of China and the Education Committee of China. The second author would also like to thank the FNRS (Belgium National Science Foundation) for its support during his sabbatical leave at Dalian University of Technology. REFERENCES 1. Z. P.Bazant, T. Belytschko and T. P.Chang, ‘Continuum theory for strain softening’, J . Eng. Mech. ASCE, 110, 1666-1692 (1984). 2. Z. P.Bazant, ‘Imbricate continuum and its variational derivation’, J. Eng. Mech. ASCE, 110,1693-1712 (1984). 3. Z. P. Bazant and L. Cedolin, Stability of Struczures: Elastic, Inelastic, Fracture and Damage Theories, Oxford University Press, New York, 1991. 4. E. C. Aifantis, ‘On the microstructural origin of certain inelastic models’, Trans. ASME. J . Eng. Mat. Tech. 106, 326-330 (1986). 5. E. C. Aifantis, ‘The physics of plastic deformation’, Int. J. Plasticity, 3,211-247 (1988). 6. D.Lasry and T. Belytschko, ‘Localization limiters in transient problems’, Int. J. Solids Struct., 24, 581-597 (1988). 7. H.B. Muhlhaus and E. C. Aifantis, ‘A variational principle for gradient plasticity’, I#. J. Solids S m r . , 29,845-857 (1991). GRADIENT PLASTICITY AT LARGE STRAINS 8. 633 R. de Borst and H. B. Muhlhaus, ‘Gradient-dependent plasticity: formulation and algorithmic aspects’, Int. j . numer. methods eng., 35, 521-539 (1992). 9. L. J. Sluys, R.de Borst and H. B. Muhlhaus, ‘Wave propagation, localization and dispersion in a gradient-dependent medium’, Int. J . Solids Srruct., 30, 1153-1171 (1993). 10. H. B. Muhlhaus, ‘Application of Cosserat theory in numerical solutions of limit load problems’, Ing.-Archiv, 59, 124-137 (1989). 11. A. Needleman, ‘Material rate dependence and mesh sensitivity in localization problems’, Comp. Methods Appl. Mech. Eng., 67, 69-86 (1988). 12. Ph. Jetteur and S. Cescotto, ‘A mixed finite element for the analysis of large inelastic strains’, 1nt.j. numer. methods eng., 31,229-239 (1991). 13. J. C. S i o and R. L. Taylor, ‘Consistent tangent operations for rate independent plasticity’, Comp. Methods Appl. Mech. Eng., 48, 101-118 (1985). 14. 0.C . Zienkiewicz and R. L. Taylor, The Finite Element Method, Vol. 1, 4th edn, Mcgraw-Hill, New York, 1989. 15. 0. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, Vol. 2,4th edn, Mcgraw-Hill, New York, 1991.