Quadrilateral elements for the solution of elasto-plastic finite strain problems

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2001; 51:883–917 Quadrilateral elements for the solution of elasto-plastic .nite strain problems Jos0eM.A.C0esar de S0a∗;†;‡;Pedro M. A. Areias§and Renato M. Natal Jorge¶ Instituto de Engenharia Mecˆ anica (IDMEC);Faculdade de Engenharia da Universidade do Porto; Rua dos Bragas;4050-123 Porto;Portugal SUMMARY In this paper two plane strain quadrilateral elements with two and four variables, are proposed. These elements are applied to the analysis of .nite strain elasto-plastic problems. The elements are based on the enhanced strain and B-bar methodologies and possess a stabilizing term. The pressure and dilatation .elds are assumed to be constant in each element’s domain and the deformation gradient is enriched with additional variables, as in the enhanced strain methodology. The formulation is deduced from a four-.eld functional, based on the imposition of two constraints: annulment of the enhanced part of the deformation gradient and the relation between the assumed dilatation and the deformation gradient determinant. The discretized form of equilibrium is presented, and the analytical linearization is deduced, to ensure the asymptotically quadratic rate of convergence in the Newton–Raphson method. The proposed formulation for the enhanced terms is carried out in the isoparametric domain and does not need the usually adopted procedure of evaluating the Jacobian matrix in the centre of the element. The elements are very e@ective for the particular class of problems analysed and do not present any locking or instability tendencies, as illustrated by various representative examples. Copyright ?2001 John Wiley & Sons, Ltd. KEY WORDS: .nite strains; elasto-plasticity; enhanced strains; .nite elements 1. INTRODUCTION The de.ciencies of the standard low-order plane strain elements, related to locking behaviour, are visible in many problems, such as in the analysis of near-incompressible situations (including elasto-plastic problems as investigated by Nagtegaal et al. [1]) and in bending-dominated situations. ∗Correspondence to: Jos0eC0esar de S0a, Department of Mechanical Engineering, University of Porto, Rua dos Bragas, 4099 Porto Cedex, Portugal †E-mail: [email protected] ‡Associate Professor §Ph.D. Student ¶Assistant Professor Contract=grant sponsor: PRAXIS; contract=grant number: XXI=BD=18533=98 Received 11 April 2000 Copyright ?2001 John Wiley & Sons, Ltd. Revised 26 October 2000 884 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE However, due to their simplicity, this class of elements is attractive for non-linear analyses, and therefore, some improvements have been carried out over the years to eliminate that undesirable behaviour. The use of higher order elements (quadratic and cubic) usually avoids these diKculties, particularly when based on a mixed formulation [2]. Nevertheless, they are more sensitive to mesh distortion [3] and the contact calculations are more diKcult. Under these conditions, properly formulated low-order elements are suitable to a wide range of applications. For the linear elastic case, there are some reliable formulations as those described in References [4–6], therefore, the non-linear case is addressed in this paper. On the contrary, for the non-linear case, the main diKculty is the trade-o@ between stability and accuracy. In fact, it is the presence of instabilities (the so-called hourglass patterns) in several important cases of compression and tension that has somehow compromised promising formulations (for example the enhanced strain element of Simo and Armero [7]). The enhanced strain method, theoretically justi.ed for non-linear problems by Simo and Armero [7], consists in the enrichment of the deformation gradient with additional variables. It is noticeable that, to satisfy the patch test, it is necessary to perform an ad hoc modi.cation of the enhanced term (classi.ed as a trick by Lautersztajn and Samuelsson in Reference [8]), based on the central evaluation of the Jacobian co-ordinate. As an extension of the standard element, the enhanced strain method is very convenient for the computational implementation of elasto-plastic material models, because the algorithmic treatment is unchanged. On the contrary, the hybrid formulations (such as presented by Pian and Sumihara [5]) and two-.eld u–pformulations (described, for example, by Brink and Stein [9]) usually demand speci.c treatment. Despite the generally good results obtained with the enhanced strain elements (for elastic problems see Reference [10]), they present some imperfections such as increased mesh distortion sensitivity and marked instabilities that appear in numerous practical problems [11;3;12]; even in the original work, Simo and Armero [7] recognize the presence of spurious modes in situations of high dilatational deformations. Owing to this drawback, several authors proposed modi.cations to the original work. Simo et al. [3] introduced a .ve-point integration rule, as opposed to the original four-point rule, with the argument that the original rule under-integrates the enhanced strain element. Notwithstanding, the instabilities are only mildly attenuated (as illustrated by Glaser and Armero [12]). Additionally, the resulting element loses Nexibility and computational eKciency. The work developed by Korelc and Wriggers [13] over a single element’s instability led to the introduction of additional orthogonality conditions and to the modi.cation of the interpolation for the deformation-gradient-enhanced part. These improvements attenuated the compression instabilities. However, the discrete strain–displacement operators lose their original sparsity. Subsequently, Glaser and Armero [12] alleged lack of objectivity in the previous proposal [13] and introduced a correction based on a further evaluation of the deformation gradient in the element’s centre. Additionally, to tackle the tension instability, a stabilization term was added to the original functional and .ve-point or nine-point integration rules were tested. Despite its e@ectiveness in controlling the instabilities, this last contribution increases the analytical complexity of the element, particularly in the evaluation of the tangent sti@ness, and further increases the computational cost due to the additional integration points. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 885 The attenuation of the mesh distortion sensitivity has also drawn attention of some researchers. For this purpose, the truncated Taylor expansion of the shape functions proposed by Wriggers and Hueck [14] and Korelc and Wriggers [15] appears to be an e@ective way of increasing the element’s robustness. It is important to note that the enhanced strain methodology is not the only strategy developed with the aim of improving the results of the standard element. There are some worth noting alternative strategies. The mean-dilatation technique of Nagtegaal et al. [1], and the closely related B-bar technique presented by Hughes [16] and further extended by Simo et al. [17] and Moran et al. [18], are appropriate to the analysis of nearly incompressible problems. These techniques are based on the assumption of an independent interpolation of the dilatation .eld. The derivations are supported by a three-.eld functional [17]. The B-bar technique can be related, under certain conditions, to the selective integration technique (as discussed by Hughes [16]). There is also an ad hoc B-bar procedure to avoid the shear locking (as discussed by Zhu and Cescotto [19]). A particularly important issue is the linearization results of the B-bar elements (addressed by Simo et al. [17] and Moran et al. [18]) that will be detailed in this paper. Although the results obtained with the B-bar formulation are frequently good, certain authors [3] have noticed a di@use response in strain localization problems. Another approach, advocated by various authors, is the reduced integration technique with a stabilizing term (as detailed in References [20–23]) which is very eKcient and avoids locking tendencies. However, this technique presents some drawbacks. In Reference [21], Reese et al. proposes a stabilization technique based on an ‘equivalent parallelogram’ that can induce a dependence between the results and the selected load increments (see Reference [21]), which can be undesirable. Besides this limitation, the single-point integration can be inadequate for elasto-plastic problems (as stated by Zhu and Cescotto [19]). There are also some important contributions regarding the mixed u−pmethods, with pressure and displacement variables (a review is presented by Brink and Stein [9]). These methods are usually based on two-.eld functionals (as adopted by Bathe [24]) or three-.eld functionals (as in the recent paper by Cris.eld [25]). This last formulation is closely related to the B-bar presented by Simo et al. [17]. Recently, some mixed=enhanced strain formulations were introduced, namely with the contributions of Piltner and Taylor [26] with added internal stress variables and Pantuso and Bathe [2] who included dilatation and pressure variables. In another context, Taylor [27] applied a mixed=enhanced formulation to triangular elements. These mixed=enhanced formulations typically imply an increase in the number of variables and therefore, are computationally more expensive than the usual enhanced strain elements. Inspired by these last two works [2;27], the authors of this paper propose two enhanced strain=B-bar elements with assumed dilatation and enhanced distortional part of the deformation gradient. These two elements extend a previous work carried out by C0esar de S0a and Natal Jorge [4]. The fundamental di@erence between the two elements is the number of additional variables: two or four. Therefore, the number of variables does not exceed the number of variables present in the standard enhanced strain element. The spurious mode control is carried out through a stabilization term included in the original functional that imposes the annulment of the enhanced part of the deformation gradient. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 886 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE The variational formulation as well as the discretization are fully detailed. The linearization calculations, necessary to obtain an asymptotically quadratic rate in the Newton–Raphson algorithm are also described. The proposed formulation is tested by several representative numerical examples, with elastic and elasto-plastic materials, covering a variety of situations. All the examples show the e@ectiveness of the two newly proposed elements and the absence of spurious modes. It is worth noting that, during the course of submission of this work, a paper describing an apparently similar enhanced strain=B-bar element formulation has been published by Armero [28]. However, as it will become clear, although there are some coincident points between the two formulations, they present some di@erences. 2. INCOMPRESSIBILITY AND LOCKING A known deformation of a body denoted as Bis said to be incompressible if any region of the body has a constant volume. This condition can be stated by the equation div v(X) = 0, which can be obtained through the conservation of mass equation, imposing 0=where is the mass density in a given con.guration and 0is the mass density in the reference con.guration. According to this de.nition, a given deformation is incompressible when the spatial divergence of the velocity is null in any point Xof B. Denoting J=det[F], Fbeing the deformation gradient, it is possible to restate the incompressibility condition as J(X)=1 (1) An equally important case, due to its practical importance, is that of near-incompressible case. The elasto-plastic materials, for isochoric plastic Nows, are frequently near-incompressible, since the dilatation is purely elastic. When numerical quadrature is used, the incompressibility condition is veri.ed only in each integration point. The numerical integration induces a relaxation in the incompressibility constraint. The locking behaviour, in the incompressible situation, of the standard four-node element is related to the absence of the two hourglass deformation modes [29;30;4]. In fact, as referred by de Souza Neto et al. [31] there is an inherent inability in the interpolation functions in the correct representation of isochoric displacement .elds. For a four node square element, C0esar de S0aet al. [4;29;30] have shown the relation between the incompressibility constraint, imposed at the four integration points, and the absence of three deformation modes: the dilatation mode and the two hourglass modes. The lack of hourglass deformation is responsible for the volumetric locking. That analysis, carried out for the linear case, is further extended here, through the annulment condition of the velocity spatial divergence: div v=0⇒@vi @xi = 0 (2) The discretized form of the velocity .eld, in each element, is stated as vi(1; 2)=Nk(1; 2)vk i, where Nkis the shape function of node k. The nodal velocity variables are denoted as vk i. The co-ordinates 1and 2are denominated local co-ordinates of the element. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 887 According to this discretized form, it is possible to write, for the referred particular case, @vi @xi (1; 2)=@Nk @i (1; 2)vk i(3) Imposing condition (2) in the standard four Gauss integration points (e.g. Reference [24]) whose local co-ordinates are denoted as k 1; k 2,k=1;:::;4, it is possible to write          div v(1 1; 1 2) div v(2 1; 2 2) div v(3 1; 3 2) div v(4 1; 4 2)          =[Q]{v}(4) where [Q] is the matrix described in References [4;30]; making use of the same notation, this matrix is de.ned for a square element in the usual local co-ordinates, i.e. with side size equal to 2, as [Q]=      −a1−a1a1−a2a2a2−a2a1 −a2−a1a2−a2a1a2−a1a1 −a1−a2a1−a1a2a1−a2a2 −a2−a2a2−a1a1a1−a1a2     (5) where a1=0:25(1 + a0);a 2=0:25(1 −a0) and a0=1=√3. The conclusions of C0esar de S0aet al., presented in References [4;29;30] keep their validity in the present .nite strain case. The point-wise incompressibility constraint for a standard square element is stated as [Q]{v}={0}(6) The hourglass modes and the dilatation mode do not belong to the space of solutions of nodal velocity vectors {v}that satisfy (6). If a reduced=selective integration is used, the incompressibility constraint is imposed in the element’s central point, and in the mean dilatation technique, the dilatation is calculated as an average of its values at the four integration points. For the selective integration case, C0esar de S0aet al. [4;29;30] have shown, based on condition (6), that the hourglass modes, previously absent, are present in this case. The [Q] matrix, corresponding to the selective integration case, can be written as [Q]=1 4[−1−11 −111 −1 1] (7) In the mean-dilatation technique, the incompressibility constraint is veri.ed in a global fashion, since the equation ve=Veis satis.ed, which states the coincidence between the material and spatial volumes in each element for any isochoric deformation. As it will become clear in this paper, the mean-dilatation method is strictly equivalent to the selective integration method, if the volume ratio is calculated exactly. Kinematically, it is possible to locally decompose the deformation gradient in a distortional term and a dilatational term (as in Reference [31]). This type of decomposition can be specially useful when there is the intent of de.ning distinct interpolations for these two terms, as it is the case in this paper. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 888 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE The kinematic decomposition can be stated as F=ˆ FFV(8) In agreement with Equation (8), the deformation gradient can be decomposed in a distortional term, such as ˆ F:det[ˆ F] = 1 and a volumetric term FV=IJ1=3(where Idenotes the identity tensor), which can also be denominated as dilatation term because it corresponds to a pure dilatation deformation. In harmony with the previous remarks, the distortional term can be calculated through the following relation: ˆ F=J−1=3F(9) The dilatation term will be considered as ‘assumed’ for the elements proposed in this paper. 3. A FOUR-FIELD FUNCTIONAL The analysis of near-incompressible problems is frequently based on functionals with two or three independent .elds. In the two-.eld case, a penalty parameter is present, related to the bulk modulus [24], and the pressure pand the spatial position xare the assumed independent .elds. In the three-.eld case, the pressure is a Lagrange multiplier and the independent .elds are: the pressure itself p, the dilatation and the spatial position x[17;32]. In this paper, a four-.eld functional is proposed. This functional is an extension of the three-.eld functional presented in Reference [17]. As usual, a particle of the body B, observed in a given reference con.guration, denoted as B, is given by its position in that con.guration, X. It is assumed that the reference con.guration Bis subjected to a conservative body force whose density is denoted by B. The surface forces are also included, identi.ed by the nominal stress vector T Tde.ned in a subset of the reference con.guration boundary, Ut⊂@B. In a subset of the boundary, Ux⊂@B, at least one component of the spatial position is prescribed. It is also assumed that Ut∩Ux=∅. A potential function associated with these conservative forces is written as Vext(x)= −B 0B·xdV−Ut T T·xdV(10) where 0is the mass density in the reference con.guration. The special class of materials considered in this paper is characterized by a stored energy function W, which can be related to the free energy function [33]. This function depends on the right Cauchy–Green deformation tensor C=FTF. The .rst Piola–Kirchho@ stress tensor, denoted as P, can be determined by di@erentiating Wwith respect to F, i.e. P=@FW, and the Kirchho@ stress tensor is related to Pby the equation =PFT. The elastic potential energy is calculated by the following expression: Vint(x)= B W(C)dV(11) Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 889 The total potential is the sum of the internal and external potentials: V(x)=V int(x)+V ext(x) (12) From the imposition of the stationarity condition for potential (12), it is possible to obtain the local form of the momentum equations and the described boundary conditions. The key idea in the following developments is the pre-established independence of the following .elds: •Spatial position: x(X) •Dilatation: (X) •Deformation gradient: F(X) As this independence is assumed, it turns out to be necessary to impose local (i.e. valid for each X∈B) constraints to relate the .elds. These constraints are added to the original form of the functional (12). The considered local constraints can be written as =det[F] (13a) F=∇0x⇔(F−∇ 0x):(F−∇ 0x) = 0 (13b) The scalar components of ∇0xare written as ∇0xij =@xi=@Xj. The second constraint, in its scalar form, allows the use of a single parameter (or a Lagrange multiplier), in contrast with the tensorial form, introduced by Simo and Armero [7] that makes use of a second order tensorial Lagrange multiplier which is the .rst Piola–Kirchho@ stress tensor. In agreement with this particular form for the constraints, the resulting functional possesses two scalar constraints. It is noticeable that, according to the multiplicative decomposition of F, de.ned by Equation (8), it is possible to introduce an assumed deformation gradient, with an assumed dilatation .eld, making use of the following notation: ˜ F=1=3ˆ F(14) where ˆ Fis the distortional part of F: ˆ F=J−1=3F(15) and J= det[F]. According to these developments and noting the de.nition J=det[F], the local constraints can be rewritten as =det[F] (16a) 3 J ˜ F−∇ 0x:3 J ˜ F−∇ 0x= 0 (16b) With this notation, the .rst constraint represents the relation J=and the second one is equivalent to ∇0x=F. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 890 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE The second constraint will be imposed by means of a penalty parameter and the .rst one by a scalar Lagrange multiplier which can be straightforwardly related to the pressure .eld, p.Ifrrepresents the penalty parameter, then, a four-.eld functional is written as V(x;F;p;)= B    W(˜ FT˜ F)+ r 2(F−∇ 0x):(F−∇ 0x)    S +p(J−)   dV+V ext(x) (17) where F=3 (J=)˜ F. A practical advantage of this functional, in comparison with the modi.ed three-.eld principle introduced in Reference [7], is the imposition of the deformation-gradient constraint without the necessity of explicitly assuming the .rst Piola–Kirchho@ stress tensor. In the case of the formulation proposed by Armero [28], the .rst Piola–Kirchho@ stress tensor must be assumed as well. This author proposes a .ve-.eld functional, to subsequently use an orthogonality condition, so that the assumed .rst Piola–Kirchho@ tensor, corresponding to a Lagrange multiplier, is absent from the .nal equilibrium equations. The present form is distinct from the one proposed by Armero [28]. The criterion for the determination of the penalty parameter will be addressed later. The inclusion of the term identi.ed as S, related to the imposition of the constraint (13b), acts as a stabilizing term, where the rparameter allows the adjustment of the relation F=∇0x. This term avoids the necessity of assumptions described by Simo and Armero [7] about the constant Piola–Kirchho@ stress .eld. It is noticeable that in the classical enhanced strain formulation this constraint is generally not veri.ed, as it is assumed that there is an orthogonality between the multiplier tensor and the constraint F=∇0x(see the assumptions established by Simo and Armero in Reference [7]). The stationarity condition of the proposed functional (17) gives the momentum equation, the boundary conditions, and if r→∞, the local imposed constraints. Imposing the annulment of the variation of functional (17), it can be stated that B3  J˜ P+r(F−∇ 0x)−1 3 J−1=3 (˜ P:F)F−T+pJF−T:FdV= 0 (18a) B1 3J J−2=3 ˜ P:F−p dV= 0 (18b) B (J−)p dV= 0 (18c) B−r(F−∇ 0x):∇0xdV+Vext(x) = 0 (18d) where ˜ P=@˜ FWis the .rst Piola–Kirchho@ stress tensor calculated with previously introduced assumed deformation gradient. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 891 These equations were calculated using the de.nition J=det[F] and the following equation: J =JF−T:F(19) From Equation (18c), it is possible to write the following relation: J=det[F]=(20) Isolating the pressure, p, from Equation (18b), the result is p=1 3J˜ P:F=1 3Jtr ˜(21) where ˜ =˜ PFT. Equation (18a) and the previous relations allow the following equation for the .rst Piola–Kirchho@ stress tensor: ˜ P=−r(F−∇ 0x) (22) This equation is somewhat predictable as it relates the introduced penalty parameter with the Lagrange multiplier used in the work of Simo and Armero [7] to impose the relation F−∇ 0x. The stress tensor ˜ Pis the cited Lagrange multiplier, assumed to be constant in the work of Simo and Armero [7]. Clearly, if r→∞, the local constraint F=∇0xshould be satis.ed. The substitution of relation (22) in Equation (18d) gives the virtual work principle in the material form. Making use of the notation ∇0x=∇xF and noting the property ˜ P:(∇xF)= ˜ PFT:∇x, it is possible to write the virtual work principle in the spatial form B :∇xdV=B 0B·xdV+Ut T T·xda(23) which is a weak form of the local momentum equation. The purpose of this exposition was to clarify the introduction functional (17), and to show that it includes the local momentum equations, boundary conditions and constraint equations. It is noticeable that there has not been any discussion about appropriate particular forms for the dilatation , pressure p, and deformation gradient F. This will be the objective of the following sections. 4. ASSUMED FIELDS The general concepts presented in Section 3 will be used to develop two plane strain elements with assumed dilatation and pressure and 2 or 4 internal variables, respectively. These elements have an additional stabilization term, that corresponds to the imposed constraint with the penalty parameter r. It will become apparent that the stabilization term is extremely simple when compared with others (for example References [20;21]), specially in the linearized form. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 898 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE 5. DISCRETIZATION The discretized form for a co-ordinate iof x, denoted as xiis written as xi=NKxk i;N kbeing the shape function of node k:Nk=1 4(1 + k 11)(1 + k 22), with 1and 2being the local co-ordinates and xk ia nodal variable. The spatial derivatives of xican be calculated as xi; j =Nk jxk i(54) and the shape function derivative is written as Nk j=@Nk @xj =@Nk @l @l @xj (55) where @l=@x jis calculated using the chain rule @x j @l =@Nm @l xm j(56) The enhanced term Acan be determined through the internal variables denoted as $k i: Aij =$k iGk j(57) Equation (57) de.nes the term Aas a linear combination of the internal variables. The terms identi.ed as Gk jin (57) are given interpolation functions. In the work of Armero [28], the form adopted for Ais far more involved, as it uses a central evaluation of the deformation gradient, and was .rst proposed by Glaser and Armero [12]. The related spatial tensor, a, can be written as aij =Aik F−1 kj =$l igl j(58) A matrix Gndim×nenh, whose dimension ndim is the spatial dimension of the problem and dimension nenh is the number of additional deformation modes, can be de.ned with the scalar components Gk j. The construction of Ais carried out in the isoparametric domain following the steps described by Simo et al. [3] but using an exact tensorial transformation. The .nal form is simply: G=J−TE(59) where Jis the Jacobian matrix of the transformation between the local co-ordinates and the material co-ordinates: i→Xj. Its scalar components can be written as Jij =@Xi=@j. It is important to note that, as shown by C0esar de S0a and Natal Jorge [4], the proposed interpolation functions satisfy the following equation: B GdV=0(60) so that the elements possess a zero mean Aover the domain, as imposed by Simo and Armero [7]. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 899 It is noticeable that this condition is satis.ed without the need for evaluating the Jacobian in the elements centre which is the case for the standard enhanced strain elements [7;3]. The interpolation matrix Eproposed in this paper has one of the two possible forms: EQi5=−211−2 2 −221−2 1 (61a) or EQi6=−211−2 2 0 0−221−2 1 (61b) written for the Qi5B-bar and Qi6B-bar elements, respectively. These denominations are an extension of those introduced in Reference [4]. As already mentioned, the mode related to the form (61a) can be identi.ed as a bubble mode. In the work of Armero [28], a two variable interpolation matrix is used. Both the interpolation matrices (61) have the additional property of being null at the elements centre. Hence, the following statement is valid: G|1=0; 2=0=0(62) This property allows us to write Equation (53) as =det[∇0x]|1=0; 2=0 (63) The divergence operators and the average divergence operators can be written in the following discretized form: div a=$l igl i(64) div a=$l igl 0i= 0 (65) div x=xl iNl i(66) div x=xl iNl 0i(67) where, according to relation (63), Nl 0k≡Nl k|1=0; 2=0 gl 0k≡gl k|1=0; 2=0 =0 (68) In agreement with these particular forms, the discretized equations derived from the equilibrium equations (52) can be written as B ij Nk jxk i+ij 3Nl 0mxl m−Nl mxl m dV+Vext = 0 (69a) B ij gk j$k i−ij 3gl m$l m+rGn rGp r$n s$p sdV= 0 (69b) An inspection of the discretized forms (69) allows the conclusion that the discrete strain operator is the so-called B-bar matrix discussed in References [17;16]. Therefore, the proposed Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 900 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE elements are of the B-bar type. The internal force vector for node kcan be straightforwardly obtained as a consequence of the discretized equations (69). 6. LINEARIZATION To ensure the asymptotically quadratic convergence rate in the Newton–Raphson method, an exact linearization of the equilibrium equations (52) is carried out. A rate analogy is used, so that the rates can be related to the iterative variations in the x and Qvariables: vXt=Xx ˙ QXt=XQ(70) where vis the spatial velocity nodal vector and ˙ Qis the internal variable time derivative. Equations (52) are linearized according to ˙ VXt∼ =−V (71) The .nal result is obtained using the chain rule. The discretized rates can be determined by some straightforward calculations, as follows. The time derivative of the Kirchho@ stress tensor is calculated using the Truesdell rate, which is denoted as ! T. The Truesdell rate of the Kirchho@ stress can be related to the time derivative of Kirchho@ stress as in Reference [36]: ˙ ˜ =! T+l˜+˜lT(72) where lis the spatial velocity gradient whose co-ordinates are lij =@vi=@xj. The Truesdell rate can be related to the deformation rate tensor through the spatial elasticity tensor, c, according to the following product: " T=c:˙U(73) where ˙ U=sym(l) is the deformation rate. The spatial velocity gradient is calculated by its de.nition (noting that lis consistent with ˜ , i.e. assumed) as follows: l=˙ ˜ F˜ F−1=˙  3−˙ J 3JI+˙ f=1 3#div v−div v+div ˙ a−div ˙ a$I+˙ f(74) where ˙ f=˙ FF−1is an auxiliary term and ˙ F=∇0v+˙ Ais the time derivative of the enhanced deformation gradient. The deformation rate is simply ˙ U=1 3(div v−div v+div ˙ a−div ˙ a)I+1 2(˙ f+˙ fT) (75) Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 901 The time derivative of ∇xcan be calculated using the chain rule: ˙ ∇x=∇0x˙ F−1=−∇x˙ f(76) and for the time derivative of a, the result is similar: ˙ a=−a˙ f(77) The assumed pressure term Tpcan be related to ˜ using Equation (50): Tp=1 3˜ P:˜ F=1 3PijFij =1 3 ik F−1 jk Fij =1 3 ik ik =1 3˜:I(78) so that its time derivative is simply: ˙ Tp=1 3˙ ˜ :I(79) The time derivatives of the divergence operators are calculated performing the contraction between the identity tensor and the related gradients. After some calculations it is possible to write the following relations: ˙ div a=−aT:˙ f(80a) ˙ div x=−∇xT:˙ f(80b) The time derivatives of the average divergence operators div xand div a, are determined using the previous relations (80) and the chain rule. For ˙ div x, the following equation can be written as ˙ div x=−˙  2VB Jdiv xdV+1 V B ˙ Jdiv xdV+1 V B J˙ div xdV(81) Using the derivative of div xin Equation (80b) and noting the following relations: ˙ =(div v+div ˙ a) (82a) and ˙ J=J(div v+div ˙ a) (82b) it is possible to write the .nal result for ˙ div x: ˙ div x=−(div v+div ˙ a)div x+1 V B (Jdiv v+Jdiv ˙ a)div xdV−1 V B J∇xT:˙ fdV(83) Equation (83) is visibly distinct from the equation derived by Simo et al. [17]. The di@erence is not only the absence of the .eld ain that paper [17], but even the linearized average spatial divergence is clearly di@erent. Introducing the approximation discussed in Section 4.5, the result is much simpler: ˙ div x=−∇xT 0:˙ f0(84a) Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 902 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE ˙ div a=−aT 0:˙ f0(84b) where the subscript 0 indicates that the quantity is evaluated at the element’s central point. Using these results, the time derivative of the Kirchho@ stress tensor can be written as ˙ ˜ =c: ˙  3−˙ J 3JI+˙ f+2 3˙  −˙ J J˜+˙ f˜+˜˙ fT(85) It is now appropriate to introduce the discretized results, necessary for the numerical implementation. The auxiliary term ˙ fcan be written as ˙ f=∇v+˙ a⇔˙ f ij =Nk jvk i+gl j˙$l i(86a) And the term ˙aij can be discretized as ˙aij =−aik ˙ f kj =−aik #Nm jvm k+gl j˙$l k$(86b) The time derivatives of the discretized forms of the spatial divergence operators and divergence averages are written as follows: ˙ div a=−$l jgl iNk jvk i−gm j˙$m i (86c) ˙ div a= 0 (86d) ˙ div x=−xl jNl iNk jvk i−gm j˙$m i (86e) ˙ div x=−xl jNl 0iNk 0jvk i(86f) Therefore, the .nal discretized form of the time derivative of the Kirchho@ stress tensor can be written as ˙ ˜ ij =cijkl %1 3kl &−gl p˙$l p+Nm 0q−Nm q vm q'+1 2(Nm kvm l+Nn lvn k+gs l˙$s k+gr k˙$r l)( +2 3&−gl p˙$l p+Nm 0q−Nm q vm q' ij +#Nm kvm i+gl k˙$l i$ kj +Nm kvm j+gl k˙$jl ik (87) and the linearized form for the assumed pressure is written as ˙ Tp=˙ ˜ ii (88) Finally, the term ˙ Ais calculated as ˙ Aij =˙$k iGk j(89) from which results a particularly simple form for the stabilizing term in the sti@ness matrix. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 903 Now the following auxiliary notations are introduced: ) ∇x=∇x+1 3(div x−div x) (90a) and * a=a+1 3(div a−div a) (90b) From which the linearized form of Equations (52) can be written as follows: B ˙U:c:) ∇x+l˜+˜lT :) ∇x +˜ :−∇x˙ f+1 3I#−∇xT 0:˙ f0+∇xT:˙ f$dV+˙ Vext =−V(x)=Xt(91a) B ˙U:c:* a+(l˜+˜lT):* a+˜:−a˙ f+1 3I#aT:˙ f$+r˙ A:AdV=−V(a)=Xt(91b) The sti@ness matrix is straightforwardly derived from Equations (91) and the discussed time derivatives. 7. STABILIZATION PARAMETER The proposed stabilization term, dependent on the rparameter, allows the control of the instability modes (i.e. hourglass modes) often present in the element formulations of the enhanced strain type. For some element formulations, based on reduced integration or uniform integration (as, for instance, the technique developed by Bonet and Bhargava [20]), there are some alternative ways of stabilizing the element which are also based on some sort of stabilizing term. The experiences carried out by the authors of this paper, with the present elements, allow the conclusion that the technique gives an eKcient control of the element’s instabilities. For the element that possesses 4 internal variables, denoted as Qi6B-bar, the following constant stabilization parameter, ris proposed: r=1 100 .(92) .being the shear modulus of the material. Usually, this value does not increase too much the element’s hourglass sti@ness to the point of damaging the intrinsic good results, but allows the control of instabilities otherwise present in several tests. The technique suggested by Bonet and Bhargava [20], was also tested with successful results, particularly for regular meshes. However, the proposed technique has a natural justi.cation for enhanced strain elements and has a variational support. It seems appropriate to refer that the Qi6B-bar element, in the absence of a stabilizing term possesses a marked stability de.ciency. However, the element with two internal variables, Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 904 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE denoted as Qi5B-bar, was tested here with a null r, with good results. A null ris equivalent to the orthogonality condition set by Simo and Armero [7] between the tensorial Lagrange multiplier and the constraint A=0. This statement can be veri.ed through Equation (22). It is clear that the imposition of a constraint by a multiplier that is assumed to be orthogonal to the constraint is equivalent to not imposing the constraint. 8. NUMERICAL EXAMPLES The following set of numerical examples attempts to validate the new proposed formulations and compare its results with those obtained using well established techniques. A variety of situations is tested, which allow to get a grasp of the elements’ behaviour, namely in compression, tension and bending situations. The presented examples are two-dimensional and obey to the plane strain hypothesis, in agreement with the discussion of Section 4.1. The Newton–Raphson method is employed, with automatic load increasing and line-search. Beginning with the second load-step, the predictor displacement is calculated as a quadratic extrapolation of the displacements, and possesses the peculiarity of including a sti@ness matrix decomposition: n+1Xu=n−1uk1+nuk2+nK−1n+1f−ni k3(93) where n+1Xuis the displacement variation foreseen for the n+ 1 step. The terms n−1uand nu are the calculated displacements at the penultimate and last steps, respectively. The nKmatrix is the tangent sti@ness matrix calculated at the end of the nth step, and the di@erence n+1f−ni=n XZfis the .rst residual corresponding to the step n+1. The parameters k1;k 2and k3are evaluated by the load factors at the steps nand n−1: n Z, n−1 Z, and by the load factor variation nXZ. Hence, k1=z1x1+z2x4+z3x7 k2=−z1+z1x2+z2x5+z3x8 k3=z1x3+z2x6+z3x9 (94) where x1= nZ2 y;x 2= n−1Z2−2n−1ZnZ y x3= n−1ZnZ2−n−1Z2nZ y;x 4=−2 nZ y x5=2 nZ y;x 6= n−1Z2−nZ2 y x7=1 y;x 8=−1 y x9= nZ−n−1Z yy=( n Z−n−1Z)2 Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 905 and z1=1 nXZ z2=z1(nZ+ nXZ) z3=z1(nZ+ nXZ)2(95) This extrapolation often allows a speed-up of the convergence rate in elasto-plastic problems, that can be irregular in the early iterations. However, if divergence occurs, this predictor is temporarily turned o@. 8.1. Material model The elasto-plastic J2model for .nite strains is based on the multiplicative decomposition of the deformation gradient and on the assumption of existence of a hyperelastic potential. The assumed deformation gradient is, in this model, decomposed in agreement with the following hypothesis: ˜ F=˜ Fe˜ Fp(96) where ˜ Feis the elastic term and ˜ Fpthe plastic term of the assumed deformation gradient. The proposed formulation corresponds to that developed by Simo in References [37;32], in its spatial form. The yield condition is given by the equation ’=s−+2 3[Y0+h(1p)]=0 where 1pis the e@ective plastic strain and sis the deviatoric part of the Kirchho@ stress tensor. The constant Y0is the initial yield stress. The strain hardening function is the one presented in Reference [17]: h(1p)=(Y∞−Y0)&1−e−1p'+H1p(97) where Y∞is the saturation stress, is the saturation exponent and His the linear hardening parameter. The hyperelastic potential has an uncoupled form and can be written with the following notation, as referred in Reference [9] (distinct from the hyperelastic potential proposed by Simo et al. [36;7]): W=3U(Je)+W(be) (98) and U(Je)=1 2[ln(Je)]2W(be)=1 2.(tr[be]−3) (99) where be=˜ Fe˜ FeT;be=Je−2=3beand Je=det[˜ Fe]. The elastic material constants 3and .are the standard bulk modulus and shear modulus, respectively. This material model is analytically simple and, due to the use of the radial return algorithm (for implementation details, consult Reference [36]) is computationally eKcient. 8.2. Nomenclature of the tested elements In the examples given here, there are 6 di@erent element formulations. The formulations are based on the calculations presented and on References [4;7]. A 2 ×2 point, Gauss integration rule is used for all the formulations. Although this type of integration Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 906 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE Figure 1. Block compression; geometry, .nite element mesh and boundary conditions. is often said to be reduced for enhanced strain formulations (as discussed by Simo et al. [3]), it is known that the instabilities do not disappear completely (as illustrated by Glaser and Armero [12]) when the number of integration points is increased. Due to this fact, the use of a higher quadrature rule was avoided, as it penalizes the computing times in plasticity problems. The standard formulation of the four-node element is identi.ed by the symbol Q4. The B-bar formulation, in agreement with the exposed formulation but without additional modes, is identi.ed by the symbol Q4B-bar. The fundamental di@erence between this formulation and the well known B-bar element is the discussed (in Section 4.5) approximation of the dilatation .eld. The well established enhanced strain formulation, originally proposed for geometrically non linear problems by Simo and Armero [7], is identi.ed by its original denomination Q1E4. The non-linear extension of the recently proposed element by C0esar de S0a and Natal Jorge [4] is denoted Qi6. The B-bar element with two internal variables, corresponding to the bubble mode, is identi.ed as Qi5B-bar. The B-bar extension of the Qi6 element is named Qi6B-bar. Finally, we note that the truncated Taylor expansion of the element-shape functions, as discussed in Reference [14] was tested. Some slight di@erences were noted for coarse meshes. However, the results presented were obtained with the standard-shape functions. 8.3. Block upsetting This problem consists of a compression of an elasto-plastic block. The corresponding model is presented in Figure 1. It is noticeable that a similar problem was inspected in a paper by Cris.eld et al. [38], for an upsetting displacement, corresponding to 30 per cent of the block’s height. In that paper the presence of hourglass instabilities was detected, although a .ve-point integration rule was used. A 5 ×10 element mesh was used, taking advantage of the existence of one symmetry plane. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 907 Figure 2. Block compression; deformed meshes of the 6 analysed element formulations. Schonauer et al. [39], using an elasto-plastic formulation based on the Hencky strain tensor, managed to reduce the instabilities presented in Reference [38], maintaining the .ve-point quadrature rule. The mesh proposed by Cris.eld et al. appears to be (as discussed by Schonauer et al. [39]), excessively coarse to model adequately the shear band zone. With the mesh used in this work, it is possible to obtain those bands with just 10 per cent upsetting, as it will be veri.ed. The material properties are the same as in Cris.eld et al., with shear modulus .=92:53 and bulk modulus 3= 200:47. The yield stress is Y0=4:81 and the material is considered to be elastic=perfectly plastic. The deformed meshes related to the 6 formulations are exhibited in Figure 2. The presence of slight instabilities, on some mesh regions, exhibited by the Qi6 and Q1E4 elements is clear. The distortions present in these regions tend to increase with the imposed vertical displacement. It is also visible that all the other elements present less distorted meshes. The Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 914 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE Figure 13. Plastic localization problem; deformed meshes and e@ective plastic strain contour plots. The test allows the examination of the plastic necking zones. The necking zones present, with certain element formulations, several shear bands inclined 45◦in relation to the loading axis. A 200 element mesh is used, with 12.826 width and 53.334 height, according to References [7;31]. An imperfection is introduced through the reduction of 1.8 per cent in the width at the specimen’s centre (as described by de Souza Neto et al. [31]). The specimen bar is discretized in 1=4 of the geometry, owing to the presence of two symmetry planes. The analysis is carried out with displacement control with a total imposed stretch of Tu=4:7. The material properties are the same as in Reference [7], with the elastic properties 3= 164:21, .=80:1983. The initial yield stress is Y0=0:45 and the hardening properties are Y∞=0:715, =16:93, H=−0:012924. Using distinct material properties, Glaser and Armero [12] detected the existence of instabilities with the enhanced strain elements, even with .veand nine-point integration rule. The deformed meshes related to the 6 examined elements are presented in Figure 13. The e@ective plastic strain contour plots are also represented. It is visible that the Qi6 and Q1E4 elements present some instabilities and the e@ective plastic strain contour plot shows some irregularity. The Q4B-bar, Qi5B-bar and Qi6B-bar elements do not show any instabilities. It is also observable that the plastic areas in elements Qi5B-bar and Qi6B-bar present a higher e@ective plastic strain concentration than the one presented by Q4B-bar element. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 915 Figure 14. Geometry, boundary conditions and material properties for the mesh distortion test. Figure 15. Vertical displacement of the nodes with applied forces, for various values of the distortion parameter. In this case the Q4 element exhibits a severe locking behaviour. 8.8. Mesh distortion test The purpose of this linear elastic test is the inspection of the bending behaviour of the Qi6 element, proposed by C0esar de S0a and Natal Jorge [4], when compared with the standard enhanced strain element, introduced by Simo and Armero [7]. This comparison is particularly relevant as it has been recently stated by Lautersztajn and Samuelsson [8] that the Qi6 element performs rather poorly for plane bending. One other important aspect is the evaluation of the mesh distortion sensitivity of the two formulations in the bending case. Therefore, the bending test presented in Reference [8] is reproduced here. Figure 14 shows the mesh, the boundary conditions and material properties of the test. In Figure 15, the vertical displacement of the nodes with applied forces is plotted and compared with the result obtained using the beam theory. It is clear that, in this test, the Qi6 element performs better than the original Q1E4 (incompatible) element. This conclusion somehow contradicts the statement referred in Reference [8]. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 916 J. M. A. C0 ESAR DE S 0 A, P. M. A. AREIAS AND R. M. N. JORGE 9. CONCLUSIONS In this paper two new elements have been described. These elements, named as Qi5B-bar and Qi6B-bar possess two and four internal variables, respectively. The enhanced deformation gradient, corresponding to an enhanced strain formulation is projected so that the dilatation is constant in the element’s domain. The dilatation .eld results from an average reduced integration, with the purpose of imposing its independence from the internal variables. The variational basis of these elements is a four-.eld functional, which contains a stabilizing term that corresponds to a constraint nullifying the distortional part of the enhanced deformation gradient. Thus, it is more suited to the enhanced strain elements than other previously proposed stabilizers (as the one presented in Reference [20]). The enhanced deformation gradient is calculated with an exact tensorial transformation, as it does not need the usually adopted [7] procedure of evaluating the Jacobian matrix in the centre of the element. With the use of a J2elasto-plastic material model, it is clear that both the Qi6 element (proposed by C0esar de S0a and Natal Jorge [4]), and the Q1E4 element (proposed by Simo and Armero [7]) present pronounced instabilities in practical .nite strain problems. It is, however, noticeable that the Qi6 element has shown superior results than the Q1E4 element in the Cook’s membrane problem. The presented results lead to the conclusion that the Qi5B-bar and Qi6B-bar are suitable to these problems, because they do not present any instabilities or any kind of locking. The Qi5B-bar formulation is particularly attractive, since it has only two internal variables, corresponding to the bubble mode, and the results obtained in the Cook’s membrane problem are nearly as good as the ones presented by the original Q1E4 element. This formulation can be extended to the analysis of tridimensional problems, but in that case the selective integration can be advantageously replaced by an exact evaluation of the assumed dilatation, as discussed by Bonet and Bhargava [20]. REFERENCES 1. Nagtegaal JC, Parks DM, Rice JR. On numerically accurate .nite element solutions in the fully plastic range. Computer Methods in Applied Mechanics and Engineering 1974; 4:153–177. 2. Pantuso D, Bathe K-J. On the stability of mixed .nite elements in large strain analysis of incompressible solids. Finite Elements in Analysis and Design 1997; 28:83–104. 3. Simo JC, Armero F, Taylor RL. Improved versions of assumed enhanced strain tri-linear elements for 3D .nite deformation problems. Computer Methods in Applied Mechanics and Engineering 1993; 110:359–386. 4. C0esar de S0a JMA, Natal Jorge RM. New enhanced strain elements for incompressible problems. International Journal for Numerical Methods in Engineering 1999; 44:229–248. 5. Pian THH, Sumihara K. Rational approach for assumed stress .nite elements. International Journal for Numerical Methods in Engineering 1984; 20:1685–1695. 6. Belytschko T, Bachrach WE. EKcient implementation of quadrilaterals with high coarse-mesh accuracy. Computer Methods in Applied Mechanics and Engineering 1986; 54:279–301. 7. Simo JC, Armero F. Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes. International Journal for Numerical Methods in Engineering 1992; 33:1413–1449. 8. Lautersztajn N, Samuelsson A. Further discussion on four-node isoparametric quadrilateral elements in plane bending. International Journal for Numerical Methods in Engineering 2000; 47:129–140. 9. Brink U, Stein E. On some mixed .nite element methods for incompressible and nearly incompressible .nite elasticity. Computational Mechanics 1996; 19:105–119. 10. Miehe C. Aspects of the formulation and .nite element implementation of large strain isotropic elasticity. International Journal for Numerical Methods in Engineering 1994; 37:1981–2004. 11. de Souza Neto EA, Peric D, Huang GC, Owen DRJ. Remarks on the stability of enhanced strain elements in .nite elasticity and elastoplasticity. In Computational Plasticity:Fundamentals and Applications. Owen DRJ, O˜nate E, Hinton E (eds), Proceedings of the Fourth International Conference,Barcelona,Spain, 1995; 1. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917 QUADRILATERAL ELEMENTS FOR FINITE STRAIN PROBLEMS 917 12. Glaser S, Armero F. On the formulation of enhanced strain .nite elements in .nite deformations. Engineering Computations 1997; 14(7):759–791. 13. Korelc J, Wriggers P. Consistent gradient formulation for a stable enhanced strain method for large deformation. Engineering Computations 1996; 13(1):103–123. 14. Wriggers P, Hueck U. A formulation of the qs6 element for large elastic deformations. International Journal for Numerical Methods in Engineering 1996; 39:1437–1454. 15. Korelc J, Wriggers P. Improved enhanced strain four-node element with Taylor expansion of the shape functions. International Journal for Numerical Methods in Engineering 1997; 40:407–421. 16. Hughes TJR. Generalization of selective integration procedures to anisotropic and nonlinear media. International Journal for Numerical Methods in Engineering 1980; 15(9):1413–1418. 17. Simo JC, Taylor RL, Pister KS. Variational and projection methods for the volume constraint in .nite deformation elasto-plasticity. Computer Methods in Applied Mechanics and Engineering 1985; 51:177–208. 18. Moran B, Ortiz M, Shih CF. Formulation of implicit .nite element methods for multiplicative .nite deformation plasticity. International Journal for Numerical Methods in Engineering 1990; 29:483–514. 19. Zhu YY, Cescotto S. Uni.ed and mixed formulation of the 4-node quadrilateral elements by assumed strain method: application to thermomechanical problems. International Journal for Numerical Methods in Engineering 1995; 38:685–716. 20. Bonet J, Bhargava P. A uniform deformation gradient hexahedron element with arti.cial hourglass control. International Journal for Numerical Methods in Engineering 1995; 38:2809–2828. 21. Reese S, Kussner M, Reddy BD. A new stabilization technique for .nite elements in non-linear elasticity. International Journal for Numerical Methods in Engineering 1999; 44:1617–1652. 22. Belytschko T, Ong JS-J. Hourglass control in linear and nonlinear problems. Computer Methods in Applied Mechanics and Engineering 1984; 43:251–276. 23. Belytschko T, Bindeman LP. Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for nonlinear problems. Computer Methods in Applied Mechanics and Engineering 1991; 88:311–340. 24. Bathe K-J. Finite Element Procedures. Prentice-Hall International Editions, 1996. 25. Cris.eld MA, Norris V. A stabilised large-strain elasto-plastic Q1=P0 method. International Journal for Numerical Methods in Engineering 1999; 46:579–592. 26. Piltner R, Taylor RL. A systematic construction of B-bar functions for linear and non-linear mixed-enhanced .nite elements for plane elasticity problems. International Journal for Numerical Methods in Engineering 1999; 44:615–639. 27. Taylor RL. Triangles, .nite elements and mixed methods. ECCM’99 European Conference on Computational Mechanics, August 31–September 3 Munchen, Germany, 1999. 28. Armero F. On the locking and stability of .nite elements in .nite deformation plane strain problems. Computers and Structures 2000; 75:261–290. 29. C0esar de S0a JMA, Owen DRJ. The imposition of the incompressibility constraint in .nite elements—A review of methods with a new insight to the locking phenomenon. In Numerical Methods for Non-Linear Problems, Taylor C et al. (eds), Proceedings of the Third International Conference, Dubrovnik. Pineridge Press Ltd. Swansea, U.K., 1986. 30. C0esar de S0a JMA. Numerical modeling of incompressible problems in glass forming and rubber technology. Ph.D. Thesis, University College of Swansea, 1986. 31. de Souza Neto EA, Peric D, Dutko M, Owen DRJ. Design of simple low order .nite elements for large strain analysis of nearly incompressible solids. International Journal for Solids and Structures 1996; 33:3277–3296. 32. Simo JC. A framework for .nite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. part ii: computational aspects. Computer Methods in Applied Mechanics and Engineering 1988; 68:1–31. 33. Lubliner J. Plasticity Theory. Macmillan Publishing Company: New York, 1990. 34. C0esar de S0a JMA, Areias PMA, Natal Jorge RM. Analysis of large elasto-plastic deformations with new compatible mode elements. ECCOMAS 2000, Barcelona, 11–14 September, 2000, accepted for publication. 35. Flanagan DP, Belytschko T. A uniform strain hexahedron and quadrilateral with orthogonal hourglass control. International Journal for Numerical Methods in Engineering 1981; 17:679–706. 36. Simo JC, Hughes TJR. Computational Inelasticity,Interdisciplinary Applied Mathematics, vol. 7. Springer: Berlin, 1998. 37. Simo JC. A framework for .nite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. part i: continuum formulation. Computer Methods in Applied Mechanics and Engineering 1988; 66:199–219. 38. Cris.eld MA, Moita GF, Jelenic G, Lyons LPR. Enhanced lower-order element formulation for large strains. In Computational Plasticity:Fundamentals and Applications, Owen DRJ, O˜nate E, Hinton E (eds), Proceedings of the Fourth International Conference, Barcelona, Spain, 1995; 1. 39. Shonauer M, de Souza Neto EA, Owen DRJ. Hencky tensor based enhanced large strain element for elastoplastic analysis. In Computational Plasticity:Fundamentals and Applications, Owen DRJ, O˜nate E, Hinton E (eds), Proceedings of the Fourth International Conference, Barcelona, Spain, 1995; 1. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2001; 51:883–917