VIEWS: 8 PAGES: 20 POSTED ON: 7/29/2011
ECCM-2001 European Conference on Computational Mechanics June 26-29, 2001 Cracow, Poland Density–dependent ﬁnite–strain plasticity. Computational issues. A. P´rez–Foguet∗ , A. Rodr´ e ıguez–Ferran and A. Huerta a Departament de Matem`tica Aplicada III, E.T.S. de Ingenieros de Caminos, Canales y Puertos e Universitat Polit`cnica de Catalunya, Jordi Girona 1, E-08034 Barcelona, Spain. e–mail: antonio.huerta@upc.es Key words: Finite strain multiplicative plasticity; Cauchy stresses; consistent tan- gent operators; geomaterial plastic models; powder compaction; Arbitrary Lagrangian– Eulerian (ALE) formulation Abstract. The consistent tangent matrix for density–dependent plastic models within the theory of isotropic multiplicative hyperelastoplasticity is presented here. Plastic equations expressed as general functions of the Kirchhoﬀ stresses and density are considered. They include the Cauchy–based plastic models as a particular case. The standard exponential return–mapping algorithm is applied, with the density playing the role of a ﬁxed parameter during the nonlinear plastic corrector problem. The consistent tangent matrix has the same structure as in the usual density–independent plastic models. A simple additional term takes into account the inﬂuence of the density on the plastic corrector problem. Quadratic convergence results are shown for several representative examples involving geomaterial and powder constitutive models. e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ 1 INTRODUCTION Consistent tangent matrices [1] are an essential ingredient for the eﬃcient solution via implicit methods of complex problems in nonlinear computational mechanics. Consis- tent tangent matrices are needed to solve the elastoplastic boundary value problem with quadratic convergence, via a full Newton-Raphson linearization. They are computed from the consistent tangent moduli at Gauss–point level. The expression of the consistent tan- gent moduli for a wide variety of material models can be found in the literature. In the context of ﬁnite–strain modelling of isotropic materials, both multiplicative hy- perelastoplasticity theory with logarithmic strain measures and time–integration based on the exponential mapping are a standard approach [2, 3]. From a computational point of view, its attractiveness stems from the fact that it leads to the same return–mapping algorithm of the inﬁnitesimal theory, and, therefore, the same linearization of the plastic corrector nonlinear equations. A key issue in the application of the multiplicative hyperelastoplasticity theory to geo- materials is the choice of the stress measure for the constitutive equations. This aﬀects directly how to model the inﬂuence of the density, or, equivalently, the volumetric de- formation, on the plastic behaviour. Recently, Meschke and Liu [4] have presented a re-formulation of the return–mapping algorithm in terms of the Cauchy stresses. They present a discussion on the role of the stress measure that needs to be chosen as the argument of the plastic model and comparative examples. The corresponding consistent tangent moduli are also included. A more general way to model the inﬂuence of the volumetric deformation on the plastic be- haviour is to consider density–dependent plastic equations, that is, yield function and ﬂow rules expressed in terms of the density and the Kirchhoﬀ stresses. Here, it is shown that the standard exponential return–mapping algorithm applies to this type of models. More- over, the consistent tangent moduli are presented. The Cauchy–based models presented in [4] are, implicitly, density–dependent plastic models due to the relationship between Cauchy and Kirchhoﬀ stresses. Therefore, the approach presented here also applies to this type of models. On the other hand, in powder compaction simulations, plasticity models expressed as a function of the stresses and the density are usual. In this context, some of these models are formulated within the multiplicative hyperelastoplasticity theory, with the Kirchhoﬀ stresses as a reference measure [5, 6]. In these cases the elastic deformation is assumed to be small with respect to the total deformation, and neglected. Speciﬁc return–mapping algorithms are devised, and, in some cases the corresponding consistent tangent moduli are presented. The approach presented here, without simpliﬁcations of the general kinematic framework, has been applied successfully to powder compaction modelling in [7]. In this reference it is shown that considering the inﬂuence of large elastic deformations in this type of problems does not represent any drawback from a modelling point of view. An outline of the paper follows. The problem is stated in section 2. The constitutive model and the standard numerical time–integration algorithm are brieﬂy reviewed. The 2 ECCM-2001, Cracow, Poland consistent linearization of the algorithm is presented in section 3. First, the expression for density–independent models is shown. After that, the extension to the density–dependent case is devised. In section 4, representative examples are discussed in detail and the convergence results are highlighted. Three elastoplastic models with diﬀerent degree of computational diﬃculty are used: a Drucker–Prager model, an elliptic model and a cone– cap model. The main conclusions are summarized in section 5. 2 PROBLEM STATEMENT This section presents a brief description of the problem. First, the constitutive equations are summarized. After that, key issues of the integration scheme are described. See Simo [2], Meschke and Liu [4] and references therein for further details on the general framework. Here, the dependence of the plastic model on the density is included in the yield function and the ﬂow rule in a general way. Moreover, it is shown that the density has no inﬂuence on the standard exponential return–mapping algorithm. 2.1 Constitutive model Let Ω ⊂ Rndim (ndim = 2, 3) be the material conﬁguration of a continuum body with particles labelled by X ∈ Ω. The motion of Ω is described by the one–parameter family of mappings ϕt : Ω → Rndim with t ∈ [0, T ] the time–like parameter. Let Ωt = ϕt (Ω) be the placement of the body at time t and x = ϕt (X) = ϕ(X, t) ∈ Ωt the position of the material particle X. In that context, the deformation gradient, ∂ϕ F (X, t) = (X, t) , (1) ∂X is assumed to be locally decomposed into elastic and plastic parts as F = F e F p . Uncoupled isotropic hyperelastic behaviour is assumed. Therefore, the local thermody- namic state is deﬁned by means of the elastic left Cauchy–Green tensor be = F e F eT and a set of strain–like scalar internal variables α ∈ Rnα (the superscript T means transpose). The Kirchhoﬀ stress tensor, τ , and the stress–like internal variables, q, are given by dW e e dW p τ =2 b and q=− , (2) dbe dα where W e and W p are the elastic and plastic parts of the free–energy function per unit of undeformed volume. The Cauchy stress tensor is given by τ σ= . (3) det (F ) The mass conservation equation reads, in material formulation, ρ0 (X) ρ(X, t) = , (4) det (F ) 3 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ where ρ(X, t) and ρ0 (X) are the densities of the particle X at times t and t = 0 respec- tively. In powder compaction modelling, equation (4) is rewritten in dimensionless format, in terms of the relative density, η(X, t) [5, 8]. Both sides of equation (4) are divided by the solid density of the compacted material, which is taken as a reference value. This leads to the equivalent expression η0 (X) η(X, t) = . (5) det (F ) The plastic response of the material is assumed isotropic and density–dependent. The dependence on the density is incorporated through η. A generic yield function f (τ , q, η) = 0 (6) and ﬂow rules Lv be = −2 γ mτ (τ , q, η) be ˙ and α = γ mq (τ , q, η) , ˙ ˙ (7) are considered, with Lv the Lie derivative with respect to the spatial velocity, v, mτ and mq the corresponding ﬂow directions, and γ the plastic multiplier. The plastic multiplier ˙ is determined with the classical Kuhn-Tucker conditions: γ≥0 , ˙ f (τ , q, η) ≤ 0 and γf (τ , q, η) = 0 . ˙ (8) Remark 1 The general constitutive model just presented includes as particular cases the plastic models expressed in terms of Cauchy stresses, Cauchy–based plastic models [4]. This can be shown by considering a generic isotropic yield function expressed on terms of the Cauchy stress tensor, f (σ, q) = 0 , (9) and using equations (3) and (5), which relate Cauchy and Kirchhoﬀ stresses through η σ = τ . (10) η0 Then, a density–dependent yield function expressed in terms of the Kirchhoﬀ stress tensor, ˜ f (τ , q, η/η0 ), can be deﬁned simply by substituting equation (10) into (9): ˜ f (σ, q) = f (τ , q, η/η0 ) = 0 . (11) ˜ The function f is not unique. Thus, the most convenient expression from a computational ˜ point of view should be chosen. On the other hand, f depends on η only through η/η0 . So an arbitrary scale can be chosen, for instance η0 = 1 for all particles X. The other components of the Cauchy–based models, i.e. hyperelastic relationships and ﬂow rules, are equal to those presented for density–dependent Kirchhoﬀ–based models, see equa- tion (7). This can be shown following the developments of Meschke and Liu [4] with a yield ˜ function f given by equation (11). Therefore, the numerical time–integration algorithm and consistent tangent moduli presented in the following apply to both types of models. 4 ECCM-2001, Cracow, Poland 2.2 Numerical time–integration The evolution from time n t to time n+1 t = n t + ∆t of the diﬀerent magnitudes associated with a prescribed material particle X is computed by means of the time–integration of the state variables be , α and η. Let us assume that n x = x(X, n t), n be , n α and n η are known values referred to time n t and that the incremental displacement n+1 ∆u = n+1 x − n x with n+1 x = x(X, n+1 t), is given. Then, the new values n+1 be , n+1 α and n+1 η for time n+1 t are computed using the corresponding values for time n t and the incremental deformation gradient, n+1 ∂ n+1 x ∂ n+1 ∆u f = = Indim + , (12) ∂ nx ∂ nx which relates the deformation gradients at times n t and n+1 t, n F and n+1 F , through the relationship n+1 F = n+1 f n F . Indim denotes a identity matrix of order ndim . The relative density, η, is integrated exactly (in the sense that no numerical time– integration scheme is used) because of the Lagrangian expression of the mass conservation principle, equation (5), which leads to n n+1 η0 η η = n+1 F ) = n+1 f ) . (13) det ( det ( The values of n+1 be and n+1 α are obtained by means of the standard elastic predictor– plastic corrector split strategy applied to equations (6–8). Remarkably, the dependence of the constitutive equations on the density does not modify the algorithm. The value of n+1 η is given by equation (13), and therefore it plays the role of a ﬁxed parameter. The result of the elastic predictor step is the so–called trial state. It is deﬁned by be = tr n+1 f n be n+1 f T and αtr = n α . (14) If the trial state is admissible, f ( τ tr , qtr , n+1 η) ≤ 0, the state at time n+1 t is set equal to the trial state. If it is not, a plastic corrector step is computed. The plastic correction step requires the approximation of the ﬂow rules, equations (7). A standard approximation consists in the use of the exponential map and the backward Euler integration scheme [2]. Under the previous isotropy assumptions, this approach leads to a nonlinear system of equations with the same structure as that of inﬁnitesimal elastoplasticity. In order to obtain the nonlinear system of equations, three vectors of Rndim are deﬁned: εe , n+1 εe and n+1 τ . The components of εe and n+1 εe are proportional to the eigenvalues tr ¯ tr e of the tensors btr and n+1 be , respectively, and the components of n+1 τ are the eigenvalues ¯ 5 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ n+1 of τ: ndim be = 2 tr exp([εe ]i ) tr ni ⊗ ni tr tr i=1 ndim n+1 e 2 b = exp([n+1 εe ]i ) ni ⊗ ni tr tr (15) i=1 ndim n+1 τ = [n+1 τ ]i ni ⊗ ni , ¯ tr tr i=1 with {ni }i=1,...,ndim the eigenvectors of the three tensors, and {[∗]i }i=1,...,ndim the compo- tr nents of the three vectors. The eigenvectors are the same for the three tensors because of the isotropy assumptions. After some manipulations, the following expression of the nonlinear system of equations is found: n+1 e ε + ∆γ mτ ( n+1 τ , ¯ ¯ n+1 q, n+1 η) = εe tr − n+1 α + ∆γ mq ( n+1 τ , ¯ n+1 q, n+1 η) = − αtr (16) n+1 n+1 n+1 f( τ, ¯ q, η) = 0 , where ∆γ = ∆t γ is the incremental plastic multiplier, mτ is the ﬂow vector in the ˙ ¯ ndim i i principal direction space (that is, mτ = i=1 [mτ ]i ntr ⊗ ntr ) and the three isotropic ¯ functions mτ , mq and f are expressed as functions of τ . Equations (16) are complemented ¯ ¯ with the restriction ∆γ ≥ 0, equation (2)2 , and e dW τ = ¯ , (17) dεe e with W (εe ) deﬁned so that equation (17) is equivalent to the hyperelastic relationship (2)1 . Once equations (16) are solved, the variables at time n+1 t are fully determined. 3 CONSISTENT TANGENT MODULI The linearization of the previous algorithm with respect to the gradient of the incremen- tal displacement is needed to solve the discrete boundary value problem with quadratic convergence. In the following, the expression of the consistent tangent moduli for density– independent plastic models is brieﬂy reviewed. After that, the general expression including the density inﬂuence is derived. 3.1 Expression for density–independent models The linearized problem is completely speciﬁed once the following consistent tangent mod- uli are given [2]: ndim ndim ndim n+1 c = [ n+1 a]ij ni tr ⊗ ni tr ⊗ nj tr ⊗ nj tr + 2 [ n+1 τ ]i ci , ¯ ˆtr (18) i=1 j=1 i=1 6 ECCM-2001, Cracow, Poland n+1 where a is a matrix of order ndim deﬁned as n+1 dn+1 τ ¯ a= e . (19) dεtr The tensors {ˆi }i=1,...,ndim can be expressed on the basis {ni }i=1,...,ndim and depend on the ctr tr e values εtr . Therefore they are fully determined from the elastic trial state [2]. Equation (18) corresponds to the linearization of the Kirchhoﬀ stresses at time n+1 t, n+1 τ , with respect to the gradient of the incremental displacements. It can be found by applying the chain rule to equation (15)3 . The ﬁrst term in the right–hand–side of equation (18) corresponds to the linearization of the eigenvalues of n+1 τ , i.e. n+1 τ . The second term ¯ i corresponds to the linearization of its eigenvectors, {ntr }i=1,...,ndim . The matrix n+1 a has an expression identical to the consistent tangent moduli of inﬁnites- e imal elastoplasticity [2]. For elastic steps, ∆γ = 0, it is equal to the Hessian of W , e n+1 d2 W a = . (20) dεe 2 t= n+1 t For plastic steps, ∆γ > 0, it can be expressed as [9] n+1 PT Gm ∇f T GP ˜ ˜ a = PT GP − ˜ , (21) ˜ ∇f T Gm with all quantities evaluated at time n+1 t and where the following notation has been intro- duced: P is a projection matrix of ndim +nα rows and ndim columns, PT = Indim , 0ndim ,nα , ˜ 0ndim ,nα is a null matrix of ndim rows and nα columns, G is a matrix of order ndim + nα −1 equal to G−1 + ∆γ∇m , with G a symmetric matrix deﬁned as e d2 W dεe 2 0 G = d2 W p , (22) 0 dα2 m is a ndim + nα vector deﬁned as mT = mT , mT , and ∇ refers to the derivatives with τ ¯ q respect to τ ¯ and q, ∂f ∂mτ¯ ∂mτ¯ ∂τ ∂τ ¯ ∂q ∇f = ¯ ∂f and ∇m = ∂mq ∂mq . (23) ∂q ∂τ ¯ ∂q 3.2 Extension to the density–dependent case In the following, the consistent linearization of the numerical–time integration scheme applied to density–dependent plastic models is presented. Although the dependence on the density of the plastic equations does not modify the standard integration algorithm, it aﬀects the consistent tangent moduli. 7 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ In a plastic step, the relative density, n+1 η, aﬀects the plastic response and thus the value of n+1 τ . It depends on the gradient of the incremental displacement through the relation n n+1 η η = ∂ n+1 ∆u , (24) det Indim + ∂nx obtained substituting equation (12) in (13). The relative density is ﬁxed during the return– mapping algorithm because it is integrated exactly (i.e. there is no need of a plastic cor- rection for n+1 η). For this reason, the standard integration algorithm remains unchanged. However, the value of n+1 τ depends on the incremental displacement through n+1 η, see equations (16) and (24). This inﬂuence has to be taken into account in the expression of the consistent tangent moduli to solve the boundary value problem with quadratic convergence. The general expression of n+1 c, equation (18), is valid for density–dependent plastic mod- els. The dependence on the density does not aﬀect the application of the chain rule nor the eigenvectors of be , {ni }i=1,...,ndim , and, therefore, neither the expression of {ˆi }i=1,...,ndim . tr tr ctr n+1 Consequently, the inﬂuence of the density is restricted to the value of τ and the matrix ¯ n+1 a. The value of n+1 τ is obtained directly from the numerical time–integration scheme, ¯ equation (16). On the other hand, n+1 a only depends on the trial state for elastic steps. In summary, only the new expression of n+1 a for plastic steps must be obtained to extend equation (18) to density–dependent plastic models. In the following, the consistent linearization of the plastic corrector step is presented. It is shown that the new, more general, expression of the moduli n+1 a is composed of two terms, one equal to the expression for density–independent models, equation (21), and another one which adds the inﬂuence of the density. In order to do that it is useful to rephrase the dependence of n+1 η on n+1 f , equation (13), as n+1 n η = η exp − tr(εe ) , ˆ tr (25) with n η = n η det ( n be ) a known value from the previous time step and where tr(∗) ˆ means trace of ∗ (that is, ndim [∗]i ). Equation (25) is found by applying the determinant i=1 function to both sides of equation (14)1 and substituting equation (15)1 into it. From the deﬁnition of n+1 a, equation (19), and the application of the chain rule a new more convenient expression of n+1 a is found: e n+1 dn+1 τ dn+1 εe ¯ d2 W dn+1 εe a = n+1 e = . (26) d ε dεe tr dεe 2 t= n+1 t dεe tr Equation (26) shows that n+1 a is determined once the total inﬂuence of εe on n+1 εe tr is found. This inﬂuence is given by the the nonlinear system of equations (16) and the relationship between n+1 η and εe , equation (25). Thus, it can be computed by linearizing tr 8 ECCM-2001, Cracow, Poland the nonlinear system of equations n+1 e ε + ∆γ mτ ( n+1 τ , ¯ ¯ n+1 q, n+1 η) = εe tr − n+1 α + ∆γ mq ( n+1 τ , ¯ n+1 q, n+1 η) = − αtr (27) f ( n+1 τ , ¯ n+1 q, n+1 η) = 0 n+1 η = n η exp − tr(εe ) ˆ tr , n+1 n+1 n+1 e n+1 and the relationships (17) and (2)2 , this links τ and ¯ q with ε and α. The expression of n+1 a is found substituting the linearization of equations (27) into equa- tion (26). It has two parts: n+1 n+1 ep n+1 ep a = aη−indep + aη . (28) The ﬁrst part, n+1 aep η−indep , is equal to the standard one for density–independent plastic models, equation (21). The second part, n+1 aep , takes into account the inﬂuence of the η density. After some manipulations it can be expressed as ˜ ˜ Gm∇f G ∂m ∂f T Gm˜ n+1 ep aη = η ∆γ PT G − ˜ 11,ndim + η P 11,ndim . (29) ∇f T Gm ∂η ˜ ∂η ˜ ∇f T Gm where 11,ndim is a matrix of one row and ndim columns with coeﬃcients equal to 1. In the density–independent case, symmetric tangent moduli are obtained for associa- tive material models, m = ∇f . On the contrary, unsymmetric moduli are found with all density–dependent material models because n+1 aep is unsymmetric except for very η particular problems (see section 4). For this reason, in density–dependent plasticity, un- symmetric linear solvers have to be used in order to keep the characteristic quadratic convergence of the Newton-Raphson method. On the other hand, it is important to remark that the expression of n+1 a for density– dependent plastic models can be computed with just a few more matrix–vector products than the standard one for density–independent plastic models (compare equations (29) and (21)), and, as expected, the additional information ∂mτ ¯ ∂mq ∂f , and . (30) ∂η ∂η ∂η 4 EXAMPLES In this section, examples with three elastoplastic models are presented: a Drucker–Prager model, an elliptic model and a cone–cap model. The three models have no internal vari- ables, nα = 0, a usual feature in the density–dependent models found in the literature. This fact simpliﬁes some equations presented in the previous section. However, the structure and the dependence on the density of the consistent tangent moduli remain unchanged, see equations (21) and (29), and, therefore, the examples are fully representative. The 9 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ 12,0 50,0 Cauchy DP model o = 46.10 Kirchhoff DP model 10,0 Cauchy DP model 40,0 Load / deformed area [kPa] Kirchhoff DP model Load / undeformed area [kPa] 8,0 30,0 o 6,0 = 27.46 20,0 o 4,0 = 46.10 o o = 27.46 =0 10,0 o 2,0 =0 0,0 0,0 0,0 0,2 0,4 0,6 0,8 0,0 0,2 0,4 0,6 0,8 Vertical displacement [cm] Vertical displacement [cm] (a) (b) Figure 1: Load–displacement results of the uniaxial compression test. Cauchy–based and Kirchhoﬀ–based Drucker–Prager models with diﬀerent friction angles, φ. (a) Load divided by deformed area and (b) load divided by undeformed area. three plastic models are complemented with Hencky’s hyperelastic law, which leads to a linear relationship between τ and εe equal to that of linear elasticity between stresses and ¯ small strains. The use of other hyperelastic laws is straightforward. The examples with the Drucker–Prager model show that the proposed approach (numeri- cal time–integration and consistent tangent moduli) is valid for Cauchy–based elastoplas- tic models [4]. The examples with the elliptic and the cone–cap models show that the consistent tangent moduli also give quadratic convergence for general density–dependent constitutive laws. The yield function and the ﬂow vector of the three models are expressed as functions of the Kirchhoﬀ stresses and the relative density, as presented in section 2. Elliptic models are widely used in porous material modelling [10], speciﬁcally in powder compaction simulations [5]. The formulation used in this work is based on that presented in [5]. There, the elastic strains are assumed small and the kinematic description of the model is simpliﬁed. A speciﬁc return mapping algorithm and the corresponding consistent tangent moduli are presented. In this work, the general ﬁnite hyperelastoplasticity frame- work is used. The results of both approaches are compared in [11]. Here, the compaction of a ﬂanged component is solved using an Arbitrary Lagrangian–Eulerian formulation for multiplicative elastoplasticity [12]. Quadratic convergence is obtained in all cases. The cone–cap model is an extension of the elliptic model. It is deﬁned by a density– dependent Drucker–Prager yield surface and a non-centered ellipse [11]. The yield function and the ﬂow rule are similar to other cone–cap models used recently in powder compaction simulations [6, 13, 8, 14]. However, the cone–cap model used here is deﬁned, like the elliptic model, in the general ﬁnite hyperelastoplasticity framework presented in section 2, without kinematic simplifying assumptions. The results obtained in a particular example are compared with those of the elliptic model. Quadratic convergence results are also shown. 10 ECCM-2001, Cracow, Poland Kirchhoff - Drucker Prager model Cauchy - Drucker Prager model (a) Iteration B = 27.46o (b) Iteration 1 2 3 4 5 6 1 2 3 4 5 6 1,E+01 1,E+01 Vertical displ. 0,01 Vertical displ. 0,01 Vertical displ. 0,05 Vertical displ. 0,05 1,E-02 Vertical displ. 0,25 1,E-02 Vertical displ. 0,25 Vertical displ. 0,50 Vertical displ. 0,50 Vertical displ. 0,75 Vertical displ. 0,75 1,E-05 1,E-05 Relative error Relative error 1,E-08 1,E-08 1,E-11 1,E-11 1,E-14 1,E-14 1,E-17 1,E-17 (c) Iteration B = 46.10o (d) Iteration 1 2 3 4 5 6 7 8 1 2 3 4 5 6 1,E+01 1,E+01 Vertical displ. 0,01 Vertical displ. 0,01 Vertical displ. 0,05 Vertical displ. 0,05 1,E-02 Vertical displ. 0,25 1,E-02 Vertical displ. 0,25 Vertical displ. 0,50 Vertical displ. 0,50 Vertical displ. 0,75 Vertical displ. 0,75 1,E-05 1,E-05 Relative error Relative error 1,E-08 1,E-08 1,E-11 1,E-11 1,E-14 1,E-14 1,E-17 1,E-17 Figure 2: Convergence results for diﬀerent vertical displacements. Kirchhoﬀ–based and Cauchy–based Drucker–Prager models (left and right respectively) with diﬀerent friction angles, φ. 4.1 Drucker–Prager model simulations The following examples involve a Cauchy–based Drucker–Prager model (CDP). They il- lustrate that the proposed approach for the numerical time–integration, standard return– mapping algorithm, and the new consistent tangent moduli, equation (28), are valid for Cauchy–based elastoplastic models. Simulations with the standard density–independent Kirchhoﬀ–based Drucker–Prager model (KDP) are also included for comparative pur- poses. Both models, CDP and KDP, are associative; the ﬂow vectors are equal to the par- tial derivatives of the corresponding yield function with respect to the Kirchhoﬀ stresses. The two models are completely diﬀerent from a modelling point of view, because diﬀerent stress measures are used, so qualitatively diﬀerent results are expected. With the KDP model, the plastic material behaviour is assumed to be density–independent. With the CDP model, on the other hand, density dependence is accounted for in an indirect fashion, because σ depends both on η and τ , see equation (10). See [4] for a discussion on the role of the stress measure in the plastic equations. A square of 1 cm of length is subjected to a plane strain uniaxial compression test [4]. The domain is modelled by a single bilinear element. A height reduction of 75% is imposed with ﬁve increments of 1% and 14 increments of 5%. The material parameters are a Young’s modulus, E, equal to 10 GPa, a Poisson’s ratio, ν, equal to 0.3, and a cohesion 11 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ Cauchy - Drucker Prager model B = 46.10o Iteration Iteration 1 2 3 4 5 6 7 8 9 10 11 1 6 11 16 21 26 31 1,E+00 1,E+00 Vertical displ. 0,01 Vertical displ. 0,45 Vertical displ. 0,05 Vertical displ. 0,50 Vertical displ. 0,10 Vertical displ. 0,55 Vertical displ. 0,25 Vertical displ. 0,60 1,E-03 1,E-03 Vertical displ. 0,65 Vertical displ. 0,40 Relative error Relative error 1,E-06 1,E-06 1,E-09 1,E-09 1,E-12 1,E-12 Figure 3: Convergence results for diﬀerent vertical displacements obtained with the stan- dard (density–independent) consistent tangent moduli, n+1 a = n+1 aepη−indep . q 100 = 1.0 = 0.9 60 = 0.8 20 = 0.55 = 0.7 = 0.6 p –120 –80 –40 0 40 80 120 Figure 4: Traces of the elliptic yield function on the meridian plane qτ –pτ for diﬀerent ¯ ¯ relative densities, η. of 2.7 3/2 GPa. The results for three values of the friction angle, φ = 0◦ , 27.46◦ , and 46.1◦ , are computed [4]. The curves of load divided by deformed and undeformed area versus vertical displacement are depicted in ﬁgure 1. The KDP model with friction angles greater than 0◦ exhibits a post–peak softening in the load divided by deformed area— vertical displacement relationship. The softening increases with the friction angle. As Meschke and Liu [4] indicate, this behaviour is related with the inelastic volume change produced by the ﬂow rule. On the other hand, the results of both models, KDP and CDP, for φ = 0◦ diﬀer only on the elastic volume change, which is not signiﬁcant for the chosen material parameters. The convergence results of both models for diﬀerent vertical displacements and two fric- tion angles are depicted in ﬁgure 2. All results are quadratic. The inﬂuence of neglecting the density contribution on the consistent tangent moduli is shown in ﬁgure 3. The con- vergence results for several load increments and the CDP model with a friction angle of 46.1◦ are depicted. In this case, non-convergence is detected after a vertical displacement of 0.65 cm. In all convergent load increments, after the initial two or three iterations, the convergence is clearly linear. The inﬂuence of the density on the consistent tangent moduli is evident comparing ﬁgures 3 and 2(d). 12 ECCM-2001, Cracow, Poland E 2 000. [MPa] φmin 45◦ ν 0.37 φmax 60◦ σy 90. [MPa] Cref 15. [MPa] η0 0.489 nc 2. n1 1. n2 2.7 Table 1: Sets of material parameters: (E, ν, σy , η0 , n1 , n2 ) for the elliptic model and (E, ν, σy , η0 , n1 , n2 , φmin , φmax , Cref , nc ) for the cone–cap model. F E 11.7 mm G C D 13.7 mm A B 6.3 mm 4.6 mm 15.6 mm Figure 5: Flanged component. Problem deﬁnition [13] and computational mesh. 4.2 Elliptic model simulations In the following, an example of density–dependent multiplicative elastoplasticity with an elliptic yield function is analyzed. The traces of the yield function on the meridian plane for diﬀerent relative densities, η, are depicted in ﬁgure 4. The ﬂow vector is equal to the partial derivatives of the yield function with respect to the Kirchhoﬀ stresses. The material parameters are summarized in table 1. They are calibrated in [5] by comparing the results of numerical simulations with experiments. The example is the frictionless compaction of a rotational ﬂanged component. The component is modelled by an axisymmetric rep- resentation as illustrated in ﬁgure 5 [13, 14]. In the simulation presented in this work the friction eﬀects are neglected. Therefore, the results are a qualitative approximation to the real compaction process, where the friction eﬀects have to be taken into account. Complete compaction simulations with this material model and friction can be found in [7]. A structured mesh of 170 eight–noded elements with reduced integration (four Gauss– points per element) is used. The radial displacement of segments BC, DE and FA and the vertical movement of segments AB and CD are set equal to zero (see ﬁgure 5). A vertical 13 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ 0.78 0.79 0.80 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 Figure 6: Elliptic model. Final relative density, η, distribution after a top displacement of 6.06 mm. Iteration Iteration 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1,E+00 1,E+00 Load level 0.6 Load level 0.7 1,E-03 1,E-03 Load level 0.8 Load level 0.9 Load level 1.0 Relative error Relative error 1,E-06 1,E-06 1,E-09 1,E-09 Load level 2.E-4 Load level 4.E-4 1,E-12 1,E-12 Load level 6.E-4 Load level 8.E-4 Load level 1.E-3 1,E-15 1,E-15 Figure 7: Elliptic model. Convergence results for diﬀerent load levels. movement of −6.06 mm is imposed to the segment EF in 80 non-uniform increments. Small increments are needed at the beginning of the test (10 of 0.01%, 19 of 0.1%, 2 of 1%) in order to ensure convergence of the boundary value problem. Larger increments are used for the rest of the test (48 of 2%). This is directly related with the high curvature of the yield surface for η close to η0 , see ﬁgure 4. The Arbitrary Lagrangian–Eulerian formulation presented in [12] for multiplicative elastoplasticity is used to reduce the mesh distortion. The mesh region ABCG is Eulerian, and equal height elements are prescribed in the mesh region CDEFG. The ﬁnal distribution of the relative density is shown in ﬁgure 6. The results are qualita- tively similar to those presented in [13], although the lack of friction produces a relatively high compacted region near point F which is not reported in similar test simulations. The convergence results for diﬀerent load levels are depicted in ﬁgure 7. Quadratic convergence is found in all cases. The ﬁnal distribution of Cauchy stresses at Gauss–point level on the meridian plane is 14 ECCM-2001, Cracow, Poland q 200 p 800 Figure 8: Elliptic model. Final distribution of Cauchy stresses at Gauss–point level on the meridian plane qσ –pσ . Blue marks indicate plastic steps and green marks indicate elastic ¯ ¯ steps. q 100 = 1.0 60 = 0.9 20 = 0.8 = 0.6 = 0.7 p 0 40 80 120 160 200 = 0.55 Figure 9: Traces of the cone–cap yield function on the meridian plane qτ –pτ for diﬀerent ¯ ¯ relative densities, η. depicted in ﬁgure 8. This ﬁgure will be useful for comparing the results obtained with this model with those of the cone–cap simulation. During the major part of the load process all Gauss points are in plastic load, blue marks; however at the end of the test a small number of Gauss points, 6, are under elastic load, green marks. On the other hand, the Gauss points subjected to higher and lower Cauchy mean pressure are located respectively above and below point C. 4.3 Cone–cap model simulations In this subsection, the frictionless compaction of a ﬂanged component (ﬁgure 5) is mod- elled with a cone–cap model. The traces of the yield function on the meridian plane for diﬀerent relative densities are depicted in ﬁgure 9. The material parameters are summa- rized in table 1. The dependence on density of the cap yield function is similar to that of the elliptic model, compare ﬁgures 9 and 4. The cone is deﬁned by a cohesion which increases up to Cref for η = 1, and a friction angle which varies parabolically from φmin at η = η0 to φmax at η = 1. The friction angle range and the dependence on the density has been set following [15], with the aim of illustrating the convergence properties of the presented consistent tangent moduli. Lower constant (i.e. density–independent) values are used in [13], and higher constant values are reported in [16]. The ﬂow vector is deﬁned so there are no grey regions on the stress space [17] and therefore no corner return–mapping algorithms [6] are needed. Complete expressions of plastic equations are presented in [11]. Note that this model represents a more demanding test 15 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ qI Load level 3.E-5 Load level 1.E-4 qI 0.5 0.5 pI pI 0.6 0.6 Load level 5.E-4 Load level 7.E-4 qI qI 0.5 0.5 pI pI 0.6 0.6 qI Load level 1.0 200 pI 800 Figure 10: Cone–cap model. Distribution of Cauchy stresses at Gauss–point level on the meridian plane qσ –pσ for diﬀerent load levels. Blue marks indicate plastic steps on the ¯ ¯ cap region, red marks plastic steps on the cone region and green marks elastic steps. for the consistent tangent moduli than the elliptic model because the dependence on the density of the yield function and the ﬂow rule is more complicated. The numerical parameters of the simulation are the same that for the elliptic model, except for the number of load increments that has been increased to 103 (smaller initial load increments have been needed to obtain convergence of the boundary value problem). The distribution of Cauchy stresses at Gauss–point level on the meridian plane for diﬀerent load levels is depicted in ﬁgure 10. Red and blue marks correspond to points under plastic regime on the cone and cap regions respectively. Green marks denote elastic regime. The inﬂuence of the cone yield surface on the stress distribution is signiﬁcant in the ﬁrst steps of the test, but it reduces drastically to 1—3 Gauss points for the rest of test. The Gauss points with plastic loading in the cone region are those of the elements located on the left side of point C. The distribution of Cauchy stresses at the end of the test (ﬁgure 10, load level 1.0) is quite similar to that obtained with the elliptic model, ﬁgure 8. The distribution of the relative density at the end of the test is shown in ﬁgure 11. The results 16 ECCM-2001, Cracow, Poland 0.78 0.79 0.80 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 Figure 11: Cone–cap model. Final relative density, η, distribution after a top displacement of 6.06 mm. 450 16 Elliptic model Rel. difference [%] 400 14 Cone-cap model Abs. difference [MPa] 350 - elliptic model 12 Mean pressure [MPa] 300 10 250 8 200 Cone-cap model 6 150 100 4 50 2 0 0 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 Top displacement [mm] Top displacement [mm] (a) (b) Figure 12: Elliptic and the cone–cap models: (a) Mean pressure of the top punch versus top displacement; (b) relative and absolute diﬀerence of the two curves. are also similar to those obtained with the elliptic model, ﬁgure 6. The inﬂuence of the cone is restricted to a slight reduction of the relative density in the ABCG region. The evolution of the mean pressure of the top punch during the test with the elliptic model and the cone–cap model are depicted in ﬁgure 12. These results agree with the previous ones: the inﬂuence of the cone is reduced. Moreover, it decreases as the compaction test progresses. No direct comparison of powder compaction simulations with elliptic and cone–cap models has been found in the literature. The cone yield surface incorporates in the constitutive model the particle sliding that occurs at low pressures [13]. According to Cante et al. [6], this eﬀect could be relevant at the transfer stage (the beginning of the industrial process, 17 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ Iteration Iteration 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1,E+00 1,E+00 Load level 3.E-5 Load level 0.5 Load level 1.E-4 Load level 0.7 1,E-03 Load level 5.E-4 1,E-03 Load level 0.9 Load level 7.E-4 Load level 1.0 Relative error Relative error 1,E-06 1,E-06 1,E-09 1,E-09 1,E-12 1,E-12 1,E-15 1,E-15 Figure 13: Cone–cap model. Convergence results for diﬀerent load levels. when an initial movement of the punches leads to the initial conﬁguration of the sample). The present results show that, as expected, the eﬀect of the cone yield surface is signiﬁcant at the beginning of the compaction process, when a low pressure is applied. However, the inﬂuence is highly reduced as the compaction process advances, and the ﬁnal results are quite similar with and without the cone (i.e.with the cone–cap model and with the elliptic model). The inﬂuence of the cone yield function when the friction eﬀects between powder and compaction tools are important has not been established yet. The convergence results for diﬀerent load levels are depicted in ﬁgure 13. These load levels include those used in ﬁgure 10 to show the distribution of Cauchy stresses at Gauss–point level. Quadratic convergence is found in all cases, for load increments with Gauss points under plastic loading in the cone region and with Gauss points under plastic loading in both regions. 5 CONCLUDING REMARKS The consistent tangent moduli for general isotropic multiplicative hyperelastic and density–dependent plastic models have been presented. Moreover, it has been shown that the inclusion of the density in the plastic equations does not modify the standard back- ward Euler return–mapping algorithm based on the exponential map [2]. Both elastic and plastic large deformations are included in the formulation. The yield function and ﬂow rules are assumed to be expressed as general functions of the Kirchhoﬀ stresses and the relative density. This includes the plastic models formulated in terms of Cauchy stresses [4] as a particular case. As in the density–independent case [2], the expression of the consistent tangent moduli is composed of two terms, the geometric part, which only depends on the trial state and the elastic energy function, and the material part, which depends on the plastic model. The inﬂuence of the density is restricted to the material part. At Gauss points under plastic loading, the material part is found to be composed of two terms. The standard one for density–independent plastic models and an additional one which includes the inﬂuence of the density on the plastic equations. The computation of the additional term only 18 ECCM-2001, Cracow, Poland requires some vector–matrix products and the derivatives of the plastic equations with respect to the relative density. Therefore, it does not represent any signiﬁcant increase of the computational cost. On the other hand, the consistent tangent moduli are, in general, unsymmetric. This has to be taken into account when solving the linear system of equations in order to obtain quadratic convergence with the Newton-Raphson method. Some representative examples with a Cauchy–based Drucker–Prager model, a density– dependent elliptic model and a density–dependent cone–cap model have been presented. In all the examples quadratic convergence has been found. Moreover, the inﬂuence of neglecting the density part of the consistent tangent moduli has been analyzed. It turns out that convergence is usually lost. When it is not, only linear convergence is obtained. A frictionless compaction test of an iron powder has been simulated, as an example of a highly demanding boundary value problem from a computational point of view. An Arbitrary Lagrangian–Eulerian formulation for multiplicative elastoplasticity [12] is used to avoid mesh distortion. The results obtained with the elliptic and the cone–cap models have been compared. Quadratic convergence has been found in both cases. Further analysis of compaction tests, including friction eﬀects, are presented in [7]. References [1] J. C. Simo and R. L. Taylor. Consistent tangent operators for rate-independent elastoplasticity. Comp. Meth. Appl. Mech. Engrg., 48, 101–118, (1985). [2] J. C. Simo. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the inﬁnitesimal theory. Comp. Meth. Appl. Mech. Engrg., 99, 61–112, (1992). c [3] D. Peri´, D. R. J. Owen, and M. E. A. Honnor. A model for ﬁnite strain elasto– plasticity based on logarithmic strains: computational issues. Comp. Meth. Appl. Mech. Engrg., 94, 35–61, (1992). [4] G. Meschke and W. N. Liu. A re-formulation of the exponential algorithm for ﬁnite strain plasticity in terms of cauchy stresses. Comp. Meth. Appl. Mech. Engrg., 173, 167–187, (1999). [5] J. Oliver, S. Oller, and J. C. Cante. A plasticity model for simulation of industrial powder compaction processes. Int. J. Sol. Struct., 33(20–22), 3161–3178, (1996). a o e [6] J. C. Cante, J. Oliver, and R. Hern´ndez. Simulaci´n num´rica de las etapas de o transferencia y prensado del proceso de compactaci´n de pulvimateriales. In IV C. e M´t. Num. Ing. SEMNI, Spain, (1999). e ıguez-Ferran, and A. Huerta. Arbitrary Lagrangian Eulerian [7] A. P´rez-Foguet, A. Rodr´ simulation of powder compaction processes with density–dependent plastic models. Submitted to Int. J. Num. Meth. Engrg., (2000). 19 e ıguez–Ferran, A. Huerta A. P´rez–Foguet, A. Rodr´ [8] J. Brandt and L. Nilsson. FE–Simulation of compaction and solid–state sintering of cemented carbides. Mech. Cohes. Fric. Mat., 3(2), 181–205, (1998). [9] M. Ortiz and J. B. Martin. Symmetry-preserving return mapping algorithms and incrementally extremal paths: a uniﬁcation of concepts. Int. J. Num. Meth. Engrg., 28, 1839–1853, (1989). [10] A. R. Ragab and Ch. A. R. Saleh. Evaluation of constitutive models for voided solids. Int. J. Plast., 15, 1041–1065, (1999). e [11] A. P´rez-Foguet, A. Rodr´ıguez-Ferran, and A. Huerta. Consistent tangent matrices for density–dependent ﬁnite plasticity models. Accepted in Mech. Cohes. Fric. Mat., (2000). [12] A. Rodr´ e ıguez-Ferran, A. P´rez-Foguet, and A. Huerta. Arbitrary Lagrangian– Eulerian (ALE) formulation for hyperelastoplasticity. Accepted in Int. J. Num. Meth. Engrg., (2000). [13] R. W. Lewis and A. R. Khoei. Numerical modelling of large deformation in metal powder forming. Comp. Meth. Appl. Mech. Engrg., 159, 291–328, (1998). [14] A. R. Khoei and R. W. Lewis. Adaptive ﬁnite element remeshing in a large defor- mation analysis of metal powder forming. Int. J. Num. Meth. Engrg., 45, 801–820, (1999). [15] K. Shinohara, M. Oida, and B. Golman. Eﬀect of particle shape on angle of internal friction by triaxial compression test. Powder Technology, 107, 131–136, (2000). e [16] P. Dor´mus and F. Toussaint. Simple tests standard procedure for the character- isation of green compacted powder. Institut National Polytechnique de Grenoble, Laboratoire Sols Solides Structures, (2000). e [17] A. P´rez-Foguet and A. Huerta. Plastic ﬂow potential for the cone region of the MRS–Lade model. J. Engrg. Mech., 125(3), 364–367, (1999). 20