Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems

Computer Methods in Applied Mechanics and Engineering 110 (1993) 359-386 North-Holland

Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems* J . C . S i m o , F. A r m e r o Division of Applied Mechanics, Department of Mechanical Engineering, Stanford University, Stanford, CA 94305, USA R.L. Taylor SEMM, Department of Civil Engineering, University of California at Berkeley, CA 94720, USA Received 1 March 1993 Improved three-dimensional tri-linear elements for finite deformation problems are developed based on an assumed enhanced strain methodology which, in the linear regime, incorporates the classical method of incompatible modes as a particular case. Three crucial modifications of a recently proposed element, which reduces to Wilson's brick in the linear regime, are introduced to prevent locking response in distorted configurations and to maintain proper rank, while preserving excellent performance in bending dominated and localization problems: (i) a modified quadrature rule; (ii) an additional enhancement of the divergence term; and (iii) a modification of the isoparametric shape function derivatives for the three-dimensional problem. Moreover, these modified shape function derivatives are shown to improve the performance of the standard tri-linear brick in distorted configurations. In addition, a strategy is described to circumvent the memory storage requirements in the static condensation procedure of the enhanced strain parameters. The excellent performance of the improved methodologyis illustrated in representative numerical simulations.

I. Introduction It is well known that the standard bi-linear quad and the tri-linear brick elements exhibit rather poor performance in finite element analysis involving bending dominated problems and locking response in the nearly incompressible limit. It is also well known that these difficulties can be circumvented by resorting to high order elements; e.g., via p-refinement. Locking response, for instance, disappears in high order ( p / > 4) triangles [1]. Nevertheless, low order elements remain quite popular in the solution of practical problems and, in particular, in finite deformation Lagrangian calculations arising in nonlinear solid mechanics. One reason is their inherent simplicity. A more important factor is the lack of robustness exhibited by high order elements in non-linear Lagrangian calculations. In contrast with purely Eulerian calculations, for Lagrangian calculations, the element Jacobian tends to be rather sensitive to the large deformations arising in finite deformation problems. In addition, adaptive mesh-generation and rezoning techniques are simpler and tend to perform better with low order elements. Motivated by these practical considerations, the development of lower order elements that improve upon the performance of the standard Galerkin method has received considerable attention in the finite element literature. Mixed finite element techniques, for instance, offer a viable alternative for circumventing the difficulties arising in the nearly incompressible limit. A fairly comprehensive review

Correspondence to: Prof. J.C. Simo, Division of Applied Mechanics, Durand Building, Stanford, CA 94305, USA. * Supported by CENTRIC Engineering Systems and NCEL under Contract No. 2-DJA-826. 0045-7825/93/$06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved

360

J.C. Simo et al., Improved versions of assumed enhanced strain

of the subject is contained in the monograph of Brezzi and Fortin [2] which includes an in-depth analysis, restricted to the linear regime, of one the most widely mixed formulations: the Q1/P0 element. For plasticity problems, this element arises in the pioneering work of Nagtegaal et al. [3] as the mean dilation approach. Mixed formulations, however, do not improve upon the bending response and, as in the Galerkin method, tend to exhibit an overly diffusive response for softening and localization problems. Moreover, the convenient 'strain-driven' format afforded by the Galerkin finite element method for local integration algorithms of general inelastic problems, is typically lost in mixed formulations. An alternative approach to improve upon low order Galerkin elements is the use of non-conventional formulations based on a local enhancement of the strain field. For strain localization problems in plasticity, for instance, approaches of this type are explored by Ortiz et al. [4], and Belytschko et al. [5] in a variational setting described in [6]. A general strategy for the construction of assumed strain finite element methods within the framework of the three-field mixed method is described in [7]. This approach includes the classical method of incompatible modes, introduced in [8] and described in [9, 10], as a particular case. The methodology allows the systematic design of enhanced elements with improved accuracy on coarse meshes comparable to the classical Wilson's brick and related elements, such as the modified QM6 of Taylor et al. [10, 21] or the one-point quadrature of Kosloff and Frazier [11]. Furthermore, as recently pointed out in [12], conventional mixed methods based on the Helinger-Reissner formulation can be brought into correspondence with assumed-strain enhanced elements. Finally, standard mixed methods such as the Q1/P0 formulation can be recovered within the assumed enhanced strain method. From a practical standpoint, a significant advantage of assumed strain methods lies in the treatment of general constitutive equations, which is identical to that arising in the Galerkin finite element method. The goal of this paper is the formulation in the full finite deformation regime of 3D enhanced tri-linear bricks that exhibit improved performance in coarse meshes, particularly in bending dominated problems, exhibit locking-free response in the incompressible limit, and are well suited for localization problems. Our approach is based on the methodology recently proposed in [13], and extends these results with reference to enhanced elements related to Wiison's brick as follows. (i) For three-dimensional problems, both the original Wilson's brick and the modified QM6 element exhibit (mild) locking in the incompressible limit, a feature pointed out in [12]. We show that locking response can be circumvented with the addition of only three additional internal modes. A general methodology for the design of enhanced modes that automatically preserve frame invariance and lead to satisfaction of a non-linear version of the patch test is also described. (ii) For distorted elements, an additional modification of the gradients of the standard shape functions is described to retain the effectiveness of preceding additional enhancement in distorted elements. The proposed modification exploits a technique for the computation of shape function derivatives used in [14] in the construction of one-point quadrature elements, which does not affect the order of convergence of the element. (iii) The conventional 8-point Gaussian quadrature in three dimensions, or the 4-point quadrature formulas in two dimensions, under-integrates these enhanced elements (and for that matter the original Wilson's brick). For finite deformation problems, the effects of this under-integration manifest themselves in the appearance of severe hourglassing in highly confined situations, particularly in the compressive regimes. A modified quadrature formula is described which circumvents these difficulties. (iv) A drawback common to both incompatible modes and enhanced assumed-strain elements lies in the fairly large local data base necessary to carry out the element static condensation procedure in non-linear problems. A modified algorithm is described which completely eliminates the need for such a large data base. An outline of the remainder of this paper is as follows. Section 2 summarizes the basic results on the standard Galerkin formulation, describes a procedure for the computation of modified shape function derivatives and motivates our subsequent developments by providing a detailed discussion of the structure taken by Galerkin approximations to the deformation gradient. Section 3 outlines the assumed enhanced strain method, describes a general procedure for the construction on enhanced deformation gradients and applies the results to specific examples, including the design of a new enhanced element.

J.C. Simo et al., Improved versions of assumed enhanced strain

361

The motivation for these new enhanced elements comes from the linear case. A detailed step-by step implementation of this enhanced element is given in Section 4, along with the description of a modified static condensation procedure that results in drastic reduction in the element data base. Representative numerical simulations that verify the results of the analysis and illustrate the good performance of the proposed methodology are described in Section 5. Conclusions are drawn in Section 6.

2. Basic notation: shape functions for the tri-linear brick

In this section we collect the basic notation employed throughout this paper and outline the standard (Lagrangian) weak form of the initial boundary value problem (IBVP) for non-linear continuum mechanics. To motivate our subsequent design of enhanced assumed strain tri-linear elements, we also describe in detail the structure of the gradients of the shape functions for the tri-linear brick and the form taken by the standard Galerkin discretization of the deformation gradient. 2. i. The standard Galerkin weak f o r m o f the I B V P

We denote by g2 C R "~m (ndi m = 2 or 3 is the space dimension) the reference configuration of a continuum body with smooth boundary F and particles labeled by X E O, with coordinates X = (X t . . . . . X,,,~m) relative to the standard Cartesian reference system in R"~'. The deformation at time t ~ [0, T] is a smooth map ~o(., t) on g) = O U F which is prescribed on a part F~ as ~o[G = ~b. We write x = ~o(X, t) = X + u(X, t ) ,

where u(., t) = displacement field.

(2.1)

We use the standard notation F = GRADx[~O ] = 13 + G R A D x [ u ] for the deformation gradient. Here 13 is the unit 3 x 3 matrix and GRADx[- ] is the gradient relative to the reference coordinates, with components O[.]/OX r In addition, we denote by P the nominal (non-symmetric) stress tensor and let T = P N be the nominal traction on the boundary F with unit normal N. The nominal traction is prescribed on a part Pr of the boundary as /O, where Fr n £~ = 0 and Fr U F, = F. To state the classical weak form of the momentum equations, we introduce the space 7/" of test functions defined in the standard fashion as (2.2)

• "= {6~(-): 6~o(X) = 0 forX E F~}.

Precise smoothness assumptions will be omitted throughout; see [15] for a detailed account of these issues. Let B and P0 designate the body force per unit of reference volume and the reference density, respectively. Denoting by (-, .) the standard L2-pairing on O of functions, vectors or tensor fields, depending on the context, the weak form of the momentum equations takes the form (po~b,&p) + (P, GRADx[6~] ) = ( B , 6 ~ ) + (T, 6~) r

V & # E °//" .

(2.3a)

Here ( . , ")r denotes the L2-inner product on the boundary F. The first and second terms on the left-hand side of (2.3) represent the virtual work of the intertial forces and the internal stresses, respectively, while the two terms on the right-hand side of (2.3a) give the virtual work of the external loads. The variational equation (2.3a) is supplemented by the initial conditions ~°l,=0 = ~Po and

VI,= o = V 0 ,

(2.3b)

where ~o0(.) and V0(.) are prescribed initial data on O. All that remains to complete the standard formulation of the weak form of the initial boundary value problem is the statement of constitutive equations for the nominal stress. For simplicity, we shall restrict most of our discussion to classical elasticity, governed by the local constitutive equations

J.C. Simo et al.. lmproved versions of assumed enhanced strain

362

P = F20c(V(C ) , where C = FtF

(2.3c)

is the right Cauchy-Green tensor and I~' is the stored energy function (a polyconvex function of the deformation gradient; see [15] for an explanation of this terminology). We remark that all the developments discussed below hold for completely general constitutive models that incorporate inelastic behavior.

2.2. The standard Galerkin discretization The numerical solution by the conventional Galerkin finite element method of the initial boundary value problem defined by (2.3a,b,c) is fairly standard and will be briefly sketched below only to motivate our subsequent developments. Comprehensive expositions can be found in one of the many standard references on the subject (e.g. [16, 17]). Consider the useful discretization of the reference placement I2 into a collection of finite element subdomains 12~, such that 12 = [,_3e_t O,, and let h denote the characteristic size of the finite element discretization. For isoparametric finite elements, the interpolation of the reference geometry and the displacement field over a typical element g2e takes the form e rtnod¢

xh(~) = E

e nnod¢

NA(~)XA

and

uh(~ :, t) = E

A=I

NA(~)UA(t) •

(2.4)

A=I

Here n~od~ is the number of nodes of element 12~ and N A : f-q--->~ are the standard isoparametric shape functions defined on the bi-unit cube: [3 = [ - 1 , 1] x [ - 1 , 1] x [ - 1 , 1]

(2.5)

(for n0im = 3) .

The interpolation over a typical element I2e of the deformation gradient, the test functions and their gradient then take the form e l! n o d e

F~(g,t)= ~

whered~ =X~ + u eA ,

d~a(t)®GRADx[NA],

A=I

(2.6)

e nnode

GRADx[&P',:] = 2

tl n o d e

6d~ ® G R A D x [ N A ] ,

where ~ ~ h = ~

A=I

NA(~)6d~ .

A=I

The formulation of the conventional finite element equations proceeds in the usual fashion and requires an explicit expression for the reference shape function gradients GRADx[NA]. The standard approach employs the chain rule relation

G R A D x [ N A] = J ( ~ ) - ' GRAD~[NA],

whereJ(~:):=

0~

(2.7)

To motivate our subsequent developments, we examine in detail the structure of the derivatives G R A D x [ N A] and the Galerkin approximation (2.6) to the deformation gradient.

2.3. Derivatives of the tri-linear isoparametric shape functions In what follows, attention is restricted to the tri-linear isoparametric brick for which n~ode = 8. Let N(~:) denote the 8 x 1 vector of associated isoparametric shape functions:

N(g)'=[N t N 2 ...

Ns] '

with

1

NA:=-~(I+~:tSC~)(I+~:2scA)(I+~:3~:A),

(2.8)

J.C. Simo et al., Improved versions of assumed enhanced strain

363

where ~:a = ( ~ , seA, ~3) are the vertices of the bi-unit cube UI. Following Belytschko et al. [141, we write N(~) in the form 3

4

N ( ~ ) = b o + ~_~ Xtb t + ~ ~ ( ' ( ~ ) y j ,

(2.9)

J=l

I~I

where ~ J ( ~ ) are hourglass-functions defined as ffl~ 1 ( ~ ) : =

~:2~:3 ,

~2(~:)

~ 3 ( ~ ) : = sc~:2 '

: = ~3~1 ,

~4(~:) : = sc s~zsc3.

(2.10)

For nai m = 2 the only hourglass function is ~ = ~:1~:z. The constant vectors b0, bt and ~,j in R 8 arising in the representation (2.9) of the shape functions can be determined via a straightforward argument (compare with the derivation given in [14]) as follows. Evaluating (2.9) at the nodes of the bi-unit cube and using the well-known property NA(~ B) = 8an, we arrive at the 8 x 8 matrix expression 3

4

18 = b 0 ® l 8 + ~] b t ® Y t + ~ y j ® h ~ , 1=1

(2.11)

J=l

where 18 is the 8 identity matrix, Y,:= [X) X~ "-" X~] t for I = 1,2, 3, and --1--

18:=

1 1 1 1 1 1

-

lm

11

'

-1 11

hi:=

'

h 2 :=

1 1

--1_

--

1--

1

-1 -1 1 -1 1 1

!-17

-1

h3:_- -111_ '

h4_- _-1 _ 111

(2.12)

_ -1_1

Inspection of these expressions reveals that 18 .h J = 0 and h J. h K = 8~.t K. Consequently, using these orthogonality properties, from the matrix expression (2.11) we arrive at the explicit results

b 0=-~

18-

1=1

(18 .YI)b,

and

39=-~

hJ-

1=1

(hJ'Yl)bt

•

(2.13a)

Finally, an easy calculation yields the additional result

bt=Jo-'~ON ~=o ( I = 1 , 2 , 3 ) ,

whereJo=J(~: ) ~=o"

(2.13b)

Expressions (2.13a,b) can be found in [14], where they are taken as the point of departure in the systematic development of one-point quadrature elements. A review of this seminal work is given in [18] and references therein. Developments that closely follow this work have been undertaken by a number of authors; see e.g., [19, 20] among others. To compute from (2.9) and (2.13a,b) the gradients of the shape functions relative to the Cartesian reference coordinates, all that is needed is an explicit expression for the Cartesian derivatives of the functions ~J(~:), J = 1 . . . . . 4. Define the differential operator GRADx[. ] by the expression: GRADx['] :=

Jot GRAD~[.],

where j ( g ) : = det[J(se)] and 1o := j(O),

(2.14)

which is to be compared with (2.7)1. Now recall that for two spatial dimensions (ndim = 2) the only

J.C. Simo et al., Improved versions of assumed enhanced strain

364

hourglass function is ~ that

=

fit 03 = b¢l b¢2 .

A direct computation that uses the preceding relation then shows

GRADx[Yg ] = GRADx[Y( ]

(for nd~m = 2) with G R A D e [ ~ ] =

(2.15)

sc~ .

This result, which appears to be not well known, does not hold for nd~m = 3 except when the Jacobian matrix J ( g ) of the element in the reference configuration is constant•

R E M A R K 1. Observe, however, that for nmm

=

3, no essential loss of accuracy occurs if GRADx[Y£ j]

is replaced by G~x[~

J] :=/.(/~ ) J o t GRADe[YgJ],

(2.16)

since both gradients coincide for elements possessing constant reference Jacobi0n. In fact, the modification embodied by (2.16) will prove essential if the enhanced tri-linear brick elements described below are to remain locking-free for distorted configurations (with non-constant reference Jacobian) in the incompressible limit.

R E M A R K 2. Relation (2.14) first appeared in [21] as an empirical modification in the computation of the incompatible shape functions of [8] that allows the resulting element (known as QM6) to pass the patch test in distorted configurations. Relation (2.14) also arises in [7, 13] in the design of enhanced strain interpolation leading to elements which satisfy the patch test by design.

2.4. The conforming Galerkin approximation to the deformation gradient With the preceding results in hand, the shape function derivatives relative to the Cartesian system in the reference configuration can be written as 4

GRADx[N A] =GRADo[N A] + E TA G R A D x [ ~ J] v constant

J=l"

(A = 1 . . . . . 8),

(2.17)

T • hourglass contribution

where yjA (j = 1. . . . . 4) are the components of the gamma-stabilization vector "yj E It~8 defined by (2.13a)2. Inserting (2.17) into (2.6), we arrive at the following result for the (conforming) Galerkin approximation to the deformation gradient: 4

8

GRADx[~P~] =GRADo[~O~'] + E gj®GRADx[ggJI . constant

J = I ~-

v hourglass contribution

with flj" = J

2

~ lAl dA

,

(2.18)

A=I

where ~o~ denotes the restriction of ~oh to element O e. It follows that the second term in (2.18) can be viewed as an enhancement to the one-point constant deformation gradient defined by F ~ - - G R A D 0 [~oh]:--GRADx [~oh][~=0. This additional term provides the 'correct rank' (stability) to one-point quadrature elements but does not increase the order of convergence of the conventional Galerkin method (see [9, 22] for precise analyses within the context of the linear problem)• The goal of the finite element method described below is precisely to improve upon the representation (2•18) by introducing a further enhancement of the deformation gradient via a mixed method• Although the convergence rate cannot be increased relative to the conventional Galerkin finite element method (or the one-point quadrature for that matter), the resulting elements will be shown to exhibit improved 'bending response' in coarse meshes, locking-free response in the incompressible limit and excellent performance in localization problems•

J.C. Simo et al.. Improved versions of assumed enhanced strain

365

3. Design of improved assumed-strain enhanced elements In this section we present a systematic procedure for the design of assumed enhanced strain interpolations which preserve the key property of frame invariance. First, we recall the key aspects involved in the formulation of assumed enhanced strain elements presented in [13]. Then, we describe a step-by-step procedure for the design of these enhanced strain interpolations and illustrate the construction with two specific examples. The second example provides an extension of Wilson's brick to the full finite deformation regime which removes locking in the nearly incompressible limit but can no longer be interpreted as an incompatible mode method.

3.1. Enhanced assumed-strain finite element method Following a standard convention, we denote by '/'~h C ~ a conforming finite element Galerkin approximation to the space ~ of admissible test function and use the notation G R A D x [ ~ h] to designate the space of associated gradients. We restrict our attention to either tri-linear (for ndi m = 3 ) Galerkin approximations or bi-linear approximations (for ndi m = 2) to the configurations and the test functions, so that the element shape functions NA(~) are defined by (2.8) (for nd~m = 3). The key idea in the formulation of an assumed strain enhanced element is to consider an improved representation of the Galerkin approximation to the deformation gradient of the form F h =GRADx[~O h] + v Galerkin

/~"

,

(3.1)

enhanced

where GRADx[~O h] is the conforming Galerkin approximation defined by expression (2.18) and ffh is an independent (additional) field which is assumed to be discontinuous across elements. In view of (2.18) and (3.1), we observe that both the hourglass term and ffh can be regarded as enhancements of the one-point-deformations gradient. Let ~ h denote the finite element space spanned by the independent enhanced deformation gradients. Taking a three-field Hu-Washizu functional in the variables {~o,F, P} as a point of departure, it can be shown that the two variational equations governing the resulting enhanced finite element method are

(p,)~o , 6~o") + (P", GRAOx[6~oh]) = (B, 6~o") + (T, 8~o")r VS~ohE ~" -. h

(3.2)

(e", 8P") = 0 w P " E Here ph is the stress field computed with the enhanced deformation gradient (3.1) via the constitutive equations. We remark that the third variational equation in this three-field mixed formulation, namely (/~h, 6 p h ) = 0, drops out as a result of the orthogonality condition between stresses and enhanced deformation gradients; see [13] for additional details. For classical elasticity p h = Fh20c(V(C n) where C h = F h t F h is the enhanced right Cauchy-Green tensor. Within the present context, the stability of the method requires the following two conditions: (i) ~h is L,-orthogonal to constant nominal stresses; i.e., (P0, 6/~h) = 0 for any P0 = constant and all 6ffh E ~h. This condition implies that fin must have a zero mean over the element. (ii) The enhanced deformation gradients generated by ~h are not contained in the conforming approximation GRADx[~h]. Equivalently, the condition G R A D x [ ~ j'] f-I ~h -----0 must hold. For the linear static theory (P0 =--0), the two preceding conditions ensure stability and convergence of the resulting finite element method (modulo technical conditions, see [7, 22]). In particular, condition (i) ensures satisfaction of the patch test. For the non-linear static theory, conditions (i) and (ii) imply that the method produces the exact solution for homogeneous deformations (i.e., deformations with constant deformation gradient); a requirement which can be viewed as a non-linear version of the classical patch test. Additionally, a further requirement must be imposed on the enhanced deformation gradients: (iii) The enhanced approximations generated by ~h must be frame-invariant in the sense that any

366

J.C. Simo et al.. Improved versions of assumed enhanced strain

1fh E ~h transforms as ffh ~ Qffh for any rotation matrix Q E SO(3), where x ~ Qx is an arbitrary superposed rigid body motion on the current configuration so that x E g,(12). Since GRADx[~p h] transforms objectively under superposed rigid body motions, this last condition ensures that the enhanced deformation gradient defined by (3.1) also transforms objectively, thus leading to a finite element method that preserves the fundamental requirement of objectivity placed on the constitutive equations. A systematic procedure for the design of interpolations for the enhanced deformation gradient that automatically satisfies the preceding requirements is described below.

3.2. Design of enhanced approximations to the deformation gradient The key idea in a convenient design of interpolations for /~h is to carry out the construction in the isoparametric domain, denoted by I--1, in terms of a transformation 0:[.] linear for fixed ~: ~ I-1, which maps vectors on [] onto vectors on []. The mapping ~:(s¢) is then transformed appropriately to obtain fib, which maps vectors in 12 onto vectors in ~o (12), in such a way that properties (i), (ii) and (iii) above are preserved. Step 1. Assume that one is given a linear map U:[.], whose range and domain are vectors defined on •, subject to the restriction

fa I:(~:) dFl = 0,

(3.3)

and such that U:(~) is not contained in GRAD_~[Wh]. In other words, condition (ii) is satisfied by construction in the isoparametric domain I-7. Step 2. Transform ~:[-] to a mapping [:~[-], whose range and domain consists of vectors defined on the reference configuration 12, via the tensorial transformation J0 J U : ' " J - 10. 0

(3.4)

By construction, since dO = j ( s¢) dV1, this transformation has the property

fa~-e(K)d12=JoJo[fm IF(s¢)dV1]Jo 1 = 0 . The transformation rule in (3.4) is motivated by the result (2.15) for ndi m = 2. ,tO

F3 ) F.

Fig. I. Design of enhanced approximations to the deformation gradient.

(3.5)

J.C. Simo et al., Improved versions o f assumed enhanced strain

367

~ h

Step 3. Transform U:[.] to obtain the two-point linear transformation F , which maps vectors in O onto vectors in ~ph(o), via the formula ph

h~

=F00:,(g),

(3.6)

whereFa:=GRADx['P"lle=0

is the one-point constant deformation gradient. Since F h transforms objectively, (3.6) ensures that the frame-invariance property (iii) is satisfied. Furthermore, combining (3.6) and (3.5), we conclude that /~" dO - F 0 ¢

(3.7)

D:e(~)dO = 0, e

i.e., the enhanced deformation gradient interpolation satisfies condition (i). Condition (ii) is satisfied by design of IF(s~). It follows that the space ~h of enhanced deformation gradients is completely determined once the mapping D:(~) defined on the isoparametric domain [] is constructed. Below we present two representative examples of this construction. 3.3. Example 1: Enhanced tri-linear brick This example generalizes to the non-linear regime the Wilson's-brick, analyzed in [23], following the approach described in [13]. Recall that Wilson's incompatible shape functions for the tri-linear brick are / V ' : = ~ [ ~,: , _, - 1] ,

/~:=_

~, - [ ~ :_=, -

1],

~ ' : = ~_[~:3, 2 1].

(3.8)

Now define the map 0:(-) on [] by the interpolation formula 3

U:(~) := ~ F~ @ GRADe[~ff ] with F~ E R 3 (I = 1, 2, 3).

(3.9)

/=1

The enhanced part of the deformation gradient is then computed according to the formula •

-j

)Foj,)a:(~)jo = I~=31 [Fo~Jor,l® [ j ( ~ Jo' GRADe[/V']]

(3.10)

Since a~ := F 0h JoF~ are constant vectors in R 3, we can adopt at in place of Ft as the independent 'local element parameters' to arrive at the interpolation 3

/~" = ~ a , ® G ~ x [ , ~ '

] , where G ~ x [ / V '

] := j ( ~ J o '

GRADe[N'I.

(3.11)

/=1

The choice (3.9) can be motivated with reference to an undistorted regular element identified with the isoparametric domain (i.e., 12e = [3) as follows. Setting GRAD0[~p"] = [f0 f0 f0], the columns of the conforming tri-linear approximation to the deformation gradient defined by (2.18) then take the form Ii

0

h

0

h

0

q9 e2 = f 2 + ~tfl3 + ~3/~l, + ~t ~:3#4 ,

(3.12)

On the other hand, the enhanced part ffh of the deformation gradient defined by expression (3.9) can be written as

J.C. Simo et al., Improved versions o f assumed enhanced strain

368

ph=[~,0/,

~:20/2 ~:3a31 "

(3.13)

Inspection of (3.12) and (3.13) reveals that the enhancement defined by (3.9) produces a modified deformation gradient whose columns now contain complete linear polynomials. The steps outlined above extend this property to distorted elements in such a way that frame invariance and satisfaction of the patch test are preserved. 3.4. Example 2: Modified enhanced tri-linear brick The further enhancement described below of the element in the preceding example can be motivated by restricting the present non-linear formulation to the infinitesimal theory. Since element distortion is taken into account by the transformation formula (3.4), it suffices to consider an undistorted element whose reference configuration coincides with the bi-unit cube; i.e., D e --F-1. Recall that under these conditions, the linearized version of the enhanced element described above is identical to the classical Wilson's brick [7]. Accordingly, in view of (3.12), the linearized change in volume AVh(~) predicted by the standard Galerkin tri-linear interpolation is of the form

ave(g) V-----~---- tr[Fh] - ndimm° + Ate:' + A2~2 +A3~3, + A 4~1 ~:2 + .A 5 ~2~3 + A 6~:3~:1 ',

(3.14)

terms arising for ndim=3

where V0 is a reference volume and A 0, A~ . . . . . A 6 are constants depending on the nodal displacements ua and the hourglass vectors ~5- Observe that the quadratic terms in (3.14) only arise in the case ndi m = 3. The linearized change in volume in two dimensions (ndi m = 2) is of the form A v h / v h = A o + A I ~ ~ + A2~72; i.e., linear in ( ~ , 62). On the other hand, the linearized change in volume predicted by the enhanced part ffh of the deformation gradient is given by AV"-'h(¢)/Vo = tr[#"l = 7,,¢, + A2~2 + A3~3,

(3.15)

where the constants A ~, A2, ,i, 3 depend on the enhanced modes 0/~, a 2 and 0/3. If one considers the linear incompressible problem and a patch of four elements, a simple argument reveals the following results: (i) For the two-dimensional problem considered in [22], enforcement of the incoml2ressibility constraint determines the constants appearing in (3.14) and (3.15) as A~ = - A ~ and A 2 = - A 2, and the method does not exhibit locking. (ii) For ndi m = 3, an analogous argument yields the additional conditions A4 = A5 = '46 = 0 which can be satisfied only for zero nodal displacements, thus leading to an element that exhibits locking in the incompressible limit. Numerical experiments indicate that this form of locking for the threedimensional Wilson's brick, first identified in [12], is fairly mild and arises only in severely constrained situations. Motivated by the preceding considerations, for ndi m = 3 we modify the enhanced part of the deformation gradient in (3.9) according to the interpolation formula 3

U:(g) := ~ F~ ® GRAD~[Nt] + (F4 • GRAD~[N4])I3,

(3.16)

1=1

where F~ E R 3 (I = 1 . . . . . 4), 13 is the 3 x 3 unit matrix and the shape function ~/4 : I-l---> R which defines the additional enhancement is given by

F~2~3-]

~4 := ~1~2~3 so that G R A D e [ N 4] = | ~ 3 ~ 1 | •

L~, ~2J

(3.17)

369

J.C. Simo et al.. Improved versions of assumed enhanced strain

In contrast with (3.15), the linearized change in volume AI~ h associated with (3.17) now contains the additional term A4 ~t ~2 + A5 ~2~3 + d 6 ~ 3 ~1. For the linear problem, enforcement of the incompressibility constraint then results in the cancellation properties Aj + Aj = 0 for J = 1 . . . . . 6, which produce an element free from locking in the incompressible limit. For distorted elements, this cancellation takes place only if the transformation between G R A D ¢ [ ~ t] and G R A D x [ ~ t] is of. the same form as the transformation rule (3.4) adopted for the enhanced deformation gradient. As pointed out in Section 2.3 this is the case only for ndi m 2. Consequently, to retain the preceding cancellation property for ndi m 3, we adopt the modification indicated in (2.16) for the computation of GRADx[Yft]. Recall that this modification does not alter the order of convergence of the method. All that remains is to derive the expression for i f associated with the map U:(~:) defined by (3.17). Using the identity =

=

10 (F4. ORAD~[~/4]) = a4" G ~ X [ / ~ 4 ] -

j(~)

with a 4 :=JoF4

(3.18)

along with the transformation formula i f " = [Jo/J(g)] GRADo[~°h]j0U:(~:) resulting from the general procedure described above, yields the final result 3

p " = ~ at ® G ~ x [ N ' /=l

(3.19)

l + (a 4 • G ~ x [ N 4 I ) G R A D o [ ~ o h ] .

The actual implementation of the resulting assumed enhanced strain method is described in detail in the following section.

R E M A R K 3. For the infinitesimal theory, classical mixed finite elements can be recovered as assumed enhanced strain elements. For instance, for nd~m = 2 the enhanced interpolations

yield an element closely related to the Q 1 / P 0 element. To see this, it suffices to note that the volume change predicted by the enhanced part of the deformation gradient is again given by (3.15) with constants A~ = a .e 2, ,~,_ = a .e I and A a. = 0. An argument analogous to that outlined above then shows that this element does not exhibit locking in the incompressible limit.

R E M A R K 4. The counterpart of the Q 1 / P 0 (constant volume) element for ndi m 3 can be obtained within the context of the infinitesimal theory by considering the enhanced interpolations =

[F'I

IF(~:)= at" ~:, 12+ot2.i~'.~'/l~, U:d

(3.21)

L I 2JJ "

where a~ and a 2 are constant vectors in R 3.

3.5. Quadrature rules and numerical integration A key point to be considered in the implementation of the enhanced elements described above concerns the numerical quadrature rule to be employed. The two key observations to be made are: (i) The conventional 2 × 2 Gaussian quadrature for nd~m = 2, or the 2 x 2 × 2 Gaussian quadrature for n~im = 3, which are typically used in the implementation of the standard Galerkin method, yield reduced quadrature formulae for the enhanced elements in Examples 1 and 2. (ii) Numerical experiments indicate that no further improvement is obtained, even for inelastic

370

J.C. Simo et al.. Improved versions of assumed enhanced strain

constitutive models, if the quadrature is increased beyond a 3 × 3-point Gaussian formula for ndi m = 2 or a 3 x 3 × 3-point Gaussian formula for ndi m = 3. The preceding observations are confirmed by numerical experiments both in the linear and the non-linear regime. In particular, for finite deformation problems involving highly confined situations in compressive regimes, the effect of under-integration (observation (i)) manifests itself by the presence of severe hourglassing that may lead to nearly useless results. In practice, the following 9-point quadrature formula for ndi m --3, or the corresponding 5-point formula for nd~m = 2, was found to produce nearly identical numerical results to those obtained with the 27-point (ndim = 3) or 9-point (ndim = 2) Gaussian quadrature formulae in observation (ii),

fO

f ( ~ ) d7-1= ~8 WLf(~L)+Wof(~o),

(3.22)

L=I

where the weights are W0 = 32/9, Wz, = 5 / 9 for L = 1 . . . . . 8 and the quadrature points ~L are the vertices of the scaled bi-unit cube [ - a , a] × I - a , ~] × [ - a , a] with a = ~ 3 / 5 . It can be easily shown that for ndi m = 2 , this formula integrates exactly all the monomials up to fourth order in the Pascal triangle with the exception of sc2tsc~ which, incidentally, is integrated exactly by the standard 8-point formula. The reduction in computational effort resulting from the quadrature formula (3.22), relative to the conventional 27-point Gaussian quadrature formula, is clearly substantial.

4. Implementation aspects We describe below a step-by-step implementation of the improved enhanced tri-linear brick presented in Example 2. The key features to be emphasized are sparsity of the symmetric gradient operators (the 'B'-matrices in the finite element jargon), when the formulation is cast in the spatial description, and generality of the constitutive models allowed. In particular, the implementation of inelastic models retains the same structure as in the conventional displacement-based Galerkin finite element method.

4.1. Step 1: Shape function derivatives and enhanced deformation gradient To facilitate the implementation, it is convenient to define the following gradients of the standard shape functions NA(~), A = 1. . . . . 8, and the enhanced shape functions NJ(~), J = 1 . . . . . 4, relative to the reference coordinates GRAD0[N A] = J¢~' GRAD~ [N A]I~:=0 (one-point quadrature derivatives), ~"~x[~J]

= j(~jgt

GRAD¢[/~J],

(4.1) 4

GRADx[N A] = (1 + aq" G ~ x [ N 4 ] )

GRADo[N A] + ~] TA ~ " ~ x [ ; C g J ] . J=l

Recall that yj are the gamma-stabilization vectors defined by (2.13a)_, and ~ J are the hourglass functions defined by (2.10). Combining result (2.17) with expression (3.19) in Example 2 of Section 3, the enhanced deformation gradient is computed as (the superscript h is omitted throughout) 8

3

F = ~ dA ® G R A D x [ N A] + ~] otj ® G"-R'~x[N' ] . A=I

(4.2)

J=l

Observe that this expression differs from that arisin 8 in the implementation of Example 1 described in [13] in that G R A D x [ N A] is replaced here by GRADx[NA]. The modifications involved are:

J.C. Simo et al., Improved versions o f assumed enhanced strain

371

N

(i) The use of G R A D x [ Y ( J] in place of G R A D x [ ~ J ] . Recall that these two expressions coincide for ndi m = 2. However, for ndi m = 3 they coincide only if the isoparametric Jacobian matrix is constant. Recall further that this modification ensures locking-free response for initially distorted elements. (ii) The presence of the addmonal" " enhanced term t~4 • R ~ x [ N " 4 ] (only for ndi m = 3 ) which" prevents locking in the incompressible limit. Finally, we transform the shape function derivatives to the current configuration via the standard formulae qo[Nal = F -' G R A D o [ N a ] ,

A=I .....

V[,V"] = F -' G ~ D x [ N J I ,

J=1,2,3,

? [ N a l = F -' G R A D x [ N a ] ,

A=I .....

8, (4.3) 8.

The spatial shape function derivatives defined by (4.3) are used in the formulation of the method in the spatial description.

4.2. Step 2: Spatial formulation of the weak form Define the Kirchhoff stress tensor by the standard expression 1-:= PF', where P is the nominal stress. The directional derivative of the enhanced deformation gradient, written with the usual abuse in notation as 6F, is easily computed from (4.2) as (4.4) 3 where we have set 6 ~ : = Ej=~ ~/J6aj. Noting that the two variational equations (3.2) that define the enhanced assumed strain method are equivalent to the single equation

+
(4.5)

+
expression (4.4) along with a straightforward manipulation yields the equivalent variational equations
V6~o E ~Fh ,

(T, sym(V[Su-])) n~ = 0

V~£11~l,~tlt2, ~¢1¢3 ~ ~3 ,

(4.6)

( r : sym(Vo[~] ), a 4 • G"-ffA-Dx[N4]) a, = 0 H e r e (. , ")n, = J'n,(')(') dO denotes the L : p a i r i n g on a typical element e = 1, 2 . . . . . E, the symbol : stands for double contraction and sym(-) denotes the symmetric part of the indicated argument. Observe that the inertia and external force contributions do not enter in the local weak forms associated with the internal enhanced element modes. This result is a direct consequence of the three-field Hu-Washizu formulation (see [7]). Translated into vector notation, the variational equations (4.6) yield the non-linear semi-discrete system of equations

M[I + Fi"'(d, ~t~) = F T M , ~

int

F~ ( d , a ~ ) = 0 ,

e=l .....

(4.7) E,

where M is the mass matrix, F ext is the vector of nodal forces, d is the vector of global nodal displacements, a e = [ a t 1 . . . at4]'e is the vector of internal element modes. In addition,

J.C. Simo et al., Improved versions o f assumed enhanced strain

372

(4.8)

"I

"='

are the internal force vectors associated with the nodal displacements and the internal element modes, • • • int A ~int .t . . respectwely. The element contributions F e and Fe to the internal force vector, respectively associated with the element nodes and the internal element modes, are defined by

F~ntA = fa (~[NA])t"r d O ,

A=I .... ,8,

F ,int J = fo ( a [ ~ ] )J , d a ,

J = l . . . . . 3,

~

(4.9)

e ~int4

F e

=

f~

6-ff~x[~"](~o[~]-

,,-1 dr~,

e

In these expressions we have used the standard finite element vector/matrix convention and set T = [TI,

'1"22

'1"33

TI2

'123

"/'31] t .

(4.10)

Similarly, e0[~] is the vector associated with sym(V0[~] ) and written with the same convention as in (4.10), with the exception of a factor of 2 now present in the shear components• Finally, BIN a] and ~[N~] are the conventional (6 x 3) sparse matrices associated with V[NA] and V[~'J], respectively. ~1 According, denoting by N~ and N,.- (i = 1, 2, 3) the components of these two vectors, we have "

OlINA] =

A

-~j

0

0

Ni 0

N,~

0

0

0

N3

Nf

0

N3

NA

0

N~.

0

Ni

and

~[~"] =

0 --j

0

N,-

o

o

0 0 (4.11)

0 0 -

3

0

1-

Here A = 1. . . . . 8 and J = 1, 2, 3. We emphasize that these matrices retain the sparsity structure of the conventional Galerkin finite element method. The effect of the enhancement is reflected only in the procedure to compute their entries, via the modified shape function derivatives (4.3). The only exception occurs in the computation of the fourth enhanced mode for which, in view of (4.9)3, the equivalent discrete gradient operator is given by the rank-one (6 x 3) matrix (4.12) Clearly, this matrix is no longer sparse. Nevertheless, the bulk of the computation is performed with the sparse matrices (4.11).

4.3. Step 3: The material part of the tangent stiffness matrix The element tangent stiffness matrix associated with the static part of the semi-discrete equations (4.9) has the following structure:

J.C. Simo et al.. Improved versions of assumed enhanced strain

K

I_sym

/~m22j

material part Km

(4.13)

I_sym K~,2 -I" geometric part Kg

In each of these two matrices the first column (row) is associated displacements while the second column (row) is associated with element modes. The expressions for the submatrices in the material part of standard. Denoting by c the spatial elasticity tensor (i.e., the tensor), the entries in these submatrices are given by

-A. Kmll

:

~e

373

with the 24 components of the nodal the 12 components of the enhanced the tangent stiffness are completely linearization of the Kirchhoff stress

(~[NA])tc(~[NB]) dO,

A,12 f. (~[NA])tc(B[ ~['t]) dO, Km

(4.14)

e

lJ 22 Km

f,f

(I~[,'V'I)tc(B[NJ]) dO,

where the indices in the preceding expressions take the values A, B E {1. . . . .

8}, I, J E {1 . . . . . 4}.

R E M A R K 5. The formulation for the geometrically linear theory follows immediately by linearization of the preceding results. The simplifications involved follow from the identification of the current and reference configurations. All that is needed is to (a) set F = 13 in expressions (4.3) for the gradients of the shape functions, (b) replace V0[~p] by 13 and set (e0[~p]) = [1 1 1 0 0 0] in all the preceding expressions and (c) delete the geometric part of the stiffness matrix. Finally, the constitutive equations are evaluated with the linearized enhanced strain matrix e : = ½IF + F t] -13, where F is defined by (4.2). 4.4. Step 4: The geometric part of the tangent stiffness matrix To express the final result in a compact form, it is convenient to revert to tensor notation and regard 7 as a 3 x 3 matrix. Let AF denote the directional derivative in the direction A~o E T "h and Aaj E I~3. As in (4.4), we have (4.15)

a e := [ ¢ [ a ~ ] + q [ a ~ + (a,, 4 • G ~ x [ ~ 4 I ) ¢ , I [ ~ , I I F ,

where a f t : = E~= l from (4.4) as

Aott. In addition, let A(6F) denote the second direction derivative o f F computed

m(t~F) "= [(~g4"

~Dx[N4])~0[A¢]

+ (m~g4" G ~ x [ / Q 4 ] ) V o [ t ~ t P ] ]

F.

(4.16)

The geometric stiffness matrix then arises from the term f6del t

K, fad,./ A,,J

= f~,, r-[(~Fr-l)t(a~F-') + a(~e)e-'] d O

(4.17)

Carrying out the computations, we arrive at the following results for the entries in the submatrices of the geometric stiffness. The entries associated with the nodal displacements and the first three enhanced modes have the standard expression:

374

J.C. Simo et al., Improved versions o f assumed enhanced strain

t"

-ABt l Kg

:

Jne V[NA]" (TV[NR]) dS2

KgA~12

:

f~le V[N A] " (~-V[/~/Jl)dg2

Kg -" 22

f,, V[U']"

(4.18)

'

dO,

e

where A, B E {1 . . . . . 8} and I, J E {1,2, 3}. Finally, in view of (4.15), (4.16) and (4.17), the entries involving the fourth enhanced mode take the form of rank-one matrices given by

K Aa = f, 1. [?,,[¢](rV[NA]) + ~-?¢,[NA]]® G-'-R~Dx[/V" ] d a g 12 / ~ z = f¢ [~,,[~](~.~[~tj])]®G-"-~,v[~ql d O ,

(4.19)

The actual computation of the preceding expressions is considerably simplified by noting that the factors ~-~7[Na] and ~'V[/VJ] are precisely the integrands appearing under the integrals that define the internal force vectors (4.9).

4.5• A modified element static condensation procedure Since the residual equations (4.7)2 hold for each element e = 1. . . . . E, the most effective solution procedure is to condense out the 12 internal element modes a " to arrive at an equivalent displacement model involving only the nodal displacement vector at the global level. The difficulty here lies in the non-linearity of equations (4.7). Given initial data {d,,,a~,V,,} at time t,,, the computation of the • • • (k) e(k) (k) solution at time t,,+~ requires an iteration process in which successwe tterates {d,,+ t, a,,+ t, V,,+t} are computed to arrive at the final solution {d,,+t, al]+~,V,,+t }. Although a conventional iterative procedure (say Newton's method) can be easily formulated, see iP, t [13] for the precise details, the actual implementation requires the storage of the element residuals F~ , the 12 element modes a " and their increments in each iteration. Additionally, to avoid repeated computation of the tangent matrices in the static condensation procedure, it is convenient to store the matrix/¢~zJK'z~. It is clear that this process places significant demands on memory allocation that grow quickly with fine meshes, particularly in 3D calculations. The alternative approach described below uses a simple observation which is widely used in the local integration of inelastic models. For a given iterate "~d(k) e(k}~, 17(/,-) t--,,+ ~, a,,+ --,,+~ 1. j, we perform a local sub-iteration to solve exactly the non-linear element equations (4.7)2 associated with the internal element modes. In practice, we have found that only a few local sub-iterations are required to obtain good results. The additional computational cost involved in this technique is more than offset by the drastic reduction in the required element data base. All that is required is storage for the 12-element modes. An outline of the algorithm is given in Table 1. Several modifications of this basic idea are possible. In particular the iteration to convergence within the DO W H I L E loop can be replaced by a fixed iteration LOOP i = 1, N~o¢.t, where N~o¢,~ is a predetermined fixed number. Since the residual ~m,¢U,,,~o~) associated with the enhanced element modes etl+l may not be below the prescribed tolerance after the N~o~,~ iterations, the element residual is modified according to the formula F i"' en+l

~-->F i°' en+[

-

~(Nh'cal)[~(Nh'cal)]-IF *il2n+l

I-*~22 n+l

(4.20)

int(Nl°¢al) I

en+l

"

J.C. Simo et al.. Improved versions of assumed enhanced strain

375

Table I Local static condensation algorithm Let {d,,. V.} be the initial data at time t and {d~.k)t, V~]~} a given iterate within the interval It., t,,, ~]. Fix this iterate and compute a , ~ l for each element .0., by means of the followingsub-iteration (at the element level): Compute GRADx[N"t], G ~ x [ / V ' ] and GRADo[N" ] Compute GRADx[~O~k)~]and F,, := GRAD,,[~.o?,I,]. Initialize a2'?, = a~ and i = 1. DO WHILE IlP~"'(d;'fl,,,~2'11~)11> tol. - C o m p u t e enhanced part ~,1 i. - Compute Ft,:l+~:= GRADx[~t,k+~,I+/~', ,. - Push-forward to compute V[Na], V[/V~]and V,[Na]. - U s e constitutive equations to compute r~'l,, = P~'~IFI,:~,r. - E v a l u a t e P~',"+It (and test for convergence) - Update a,~':~'~ = a,~:; - [/~t.,~,,~~]-~p~.,,~l (if convergence is not attained). - Update sub-iteration counter i = i + 1. END DO

In practice, we have obtained successful computations by defaulting the total number of local sub-iterations to NtocaI = 2. This technique is specially well suited for vectorization. The practical advantages of this strategy are obvious due to the reduction in computational effort combined with the significant savings in the required element data base and vectorization possibilities.

5.

Representative

numerical

simulations

The goal of this section is to investigate numerically the actual performance exhibited by the improved enhanced elements described above. Specific issues of interest in this numerical assessment are a locking-free response in the incompressible limit and a good performance in localization problems, comparable to that exhibited by the 2-D counterpart of the Q 1 / E 9 element in plane strain problems. T o provide an indication of the numerical performance to be expected in actual simulations, we consider four test problems which range from linearized incompressible elasticity to finite strain coupled thermoplasticity and compare the performance of different finite element formulations for distorted and undistorted configurations. In particular, the following finite elements are considered: (1) Q1. This element corresponds to the standard tri-linear brick based on the displacement (or primal) formulation of the problem. The standard 8-point quadrature rule is employed and the modification of the gradients of the shape functions described in Section 2.3 is not introduced. This element is expected to show severe locking in the incompressible limit as well as very poor performance in localization problems. (2) QM1. This element is identical to Q1 except for the computation of the gradient of the isoparametric shape functions, which employs the modification described in Section 2.3. Consequently, these two elements should become identical in undistorted element configurations (i.e., exhibit severe locking). On the other hand, the modified shape function derivatives should improve the response in the incompressible limit in initially distorted configuration. The standard 8-point quadrature rule is employed. (3) Q1/P0. This element arises from a mixed formulation in terms of a tri-linear displacement interpolation and a (discontinuous) constant pressure approximation. For geometrically linear problems this element corresponds to the mean dilation formulation presented in [3] and fails to satisfy the inf-sup condition. The geometrically non-linear version is described in [24]. Again, the conventional 8-point quadrature rule is employed along with a standard computation of the gradients of the isoparametric shape functions. This element is expected to show no locking for the incompressible case, but a very 'diffusive' behavior in localization problems.

J.C. Simo et al., Improved versions o f assumed enhanced strain

376

(4) Q1/E9. This element corresponds to the 3-D enhanced formulation with 9 enhanced modes as described in Section 3.3. The conventional computation of the shape functions derivatives is employed, together with the standard 8-point quadrature rule. The linearized version of this element coincides with the QM6 element of [21]. Thus, as pointed out in [12], a 'mild' locking in the incompressible limit is to be expected. In the geometrically non-linear regime, the 2-D version of this element reduces to the Q 1 / E 4 element (i.e., only 4 enhanced modes) described in [13] and performs remarkably well in localization problems. Box 1 Thermoplastic model: J_,-flow theory 1. Free energy fimction: (F = F~F p, b ~ = F~F ~' and /~e = j-Z/3b,, with J = det F) ~ ( b e, a, O) = W ( b e) + U(J) + M(J, O) + T((9) + K(a, (9). (i) Logarithmic hyperelastic response, W(/I~)=/x

~

[log(XA)] 2 and

U(J)=½KlogZJ,

A=l,3

where '~m = j-I/3)~e --A (Am being the elastic principal stretches). (ii) Thermoelastic coupling (/3 = constant is the coefficient of thermal expansion), M(J, (9) = -3K/3((9 - (grc,) log J . (iii) Thermal contribution (co = constant is the heat capacity),

Vref/.s

(iv) Hardening potential (6 = constant), K(e~, (9) = ½h((9 )ot 2 + (Yo((9) -y=((9))H(e~), where H ( a ) : : { O - ( 1 - e - S a ) / 6 , ,

if6#0,if6=0.

2. Plastic response: (i) Von Mises yield criterion,

~b(¢, q, (9) =

Ildevh'lll - V~

[Yo((9) - q] ~< O.

(ii) Conjugate hardening variable, q := -Oa @ = - [ h ( ( 9 ) a + (Yo((9) - y ~ ( ( 9 ) ) ( 1 - e - ~ ) ] . (iii) Linear thermal softening (%, oJh are constants), Yo((9) = Yo((grCf)[1 - wo((9 - (9~f)], h((9 ) = h((9~,)[1 - ¢0h((9 -- (9~,)] , y~((9) = y=((9~c,)[1 -- W,,((9 -- (9~¢,)1 .

J.C. Simo et al.. Improved versions of assumed enhanced strain

377

Table 2 Thermoplastic model: material properties Bulk modulus Shear modulus Flow stress Linear hardening Saturation hardening Hardening exponent Density Expansion coefficient Conductivity Specificcapacity Flow stress softening Hardening softening

K /.t y. h y~ 8 p. /3 k c, to. to~,

164.206 GPa 80.1938 GPa 0.450 GPa 0.12924 GPa 0.715 GPa 16.93 7800 kg/m -~ 1 x 10-~K 45 N/s K 460 m'-/s'- K 2.0 × 10 3 K2.0 × 10- 3 K- t

(5) Q 1 / E 1 2 . This element corresponds to the enhanced formulation described in Section 3.4, employing a total of 12 enhanced modes. The 9-point quadrature rule described in Section 3.5 is employed, but no modification of the gradient of the isoparametric shape functions is considered. Consequently, this element is expected to exhibit locking-free response in the incompressible limit for undistorted configurations. On the other hand, a mild locking is expected to occur in distorted configurations. Observe that the 2-D plane strain counterparts of this and the Q 1 / E 9 elements are identical except for the different quadrature rules employed: the conventional 4-point formula for the latter and the 5-point formula of Section 3.5 for the former. Whether tha modified quadrature formula alters the good performance observed in localization problems is a key issue addressed below. (6) Q M 1 / E 1 2 . This element is identical to the previous Q 1 / E 1 2 except for the computation of the gradients of the isoparametric shape functions, which employs the modification described in Section 2.3. No locking response is expected to occur for either distorted or undistorted configurations in the incompressible limit. The 2D-plane strain counterpart of this element and the Q 1 / E 1 2 element are identical. The only difference between these two elements and the element labeled Q 1 / E 4 in [13] is the quadrature rule. As pointed out above, the conventional 4-point quadrature formula underintegrates the Q 1 / E 4 element. The main reason for introducing the 'intermediate' elements QM1 and Q 1 / E 1 2 is to assess the relative importance of the different modifications introduced in the preceding sections which led to the final QM1 / E l 2 element. The implementation of all of the preceding elements is carried out within the general framework of finite strain coupled multiplicative thermoplasticity presented in [28], using the algorithmic treatment first described in [27]. The specific constitutive model under consideration is summarized in Box 1. Briefly, it consists of a hyperelastic response characterized by a logarithmic free energy function with an associative flow rule and a Von Mises yield criterion. Table 2 includes the specific values of the material properties. Some of the simulations presented below are performed in the uncoupled (isothermal) case, obtained by considering a fixed temperature which is set to the reference value (i.e., 0 = Orcf). It is worth noting that, due to the strain driven structure of the proposed enhanced formulations, the implementation of coupled models as complicated as the one considered here is rather simple. In fact, the actual implementation retains the same structure arising in the standard Galerkin displacement formulation (Q1), while avoiding the cumbersome difficulties associated with mixed formulations (e.g., O l / P0). An outline of the numerical simulations described below is as follows. In Section 5.1, we consider a single element test and examine the element eigenvalues computed in both undistorted and distorted configurations. The analysis is performed with the initial tangent stiffness and, therefore, is identical to an analysis under the assumption of linearized incompressible elasticity. Section 5.2 describes the results obtained in the upsetting of a cylindrical billet in the context of finite strain (uncoupled) plasticity. Of special interest is an assessment of the performance of the enhanced elements under investigation in problems involving large compressive strains. The 2-D plane strain counterparts of these elements are investigated in Section 5.3 by means of an isothermal, plane strain tension test of an elastoplastic

J.c. Simo et al.. Improved versions of assumed enhanced strain

378

m a t e r i a l w i t h s t r a i n - s o f t e n i n g . T h e p u r p o s e o f t h i s e x a m p l e is t o a s s e s s t h e p e r f o r m a n c e o f t h e p r o p o s e d higher order quadrature

r u l e in t h i s c o n t e x t . F i n a l l y , S e c t i o n 5.4 d e s c r i b e s t h e r e s u l t s o b t a i n e d

full-coupled thermomechanical

in a

necking problem.

5.1. L i n e a r i z e d eigenvalue analysis T h e g o a l o f t h i s s i m u l a t i o n is t o a s s e s s t h e p e r f o r m a n c e

in t h e n e a r l y i n c o m p r e s s i b l e l i m i t o f t h e

d i f f e r e n t e l e m e n t s o u t l i n e d a b o v e via a s p e c t r a l a n a l y s i s o f t h e initial e l a s t i c s t i f f n e s s m a t r i x o f a s i n g l e finite e l e m e n t . T h e u n c o u p l e d m e c h a n i c a l t h e o r y is a s s u m e d w i t h a r a t i o K/gt = 1 x 109 t h a t e f f e c t i v e l y enforces the incompressibility constraint. This analysis provides an assessment of element performance in t h e l i n e a r i n c o m p r e s s i b l e e l a s t i c r e g i m e . T w o d i f f e r e n t c o n f i g u r a t i o n s a r e c o n s i d e r e d . T h e first, s h o w n in Fig. 2 ( a ) , c o r r e s p o n d s t o t h e u n d i s t o r t e d b i - u n i t i s o p a r a m e t r i c b r i c k w h i l e in t h e s e c o n d , a l s o s h o w n in Fig. 2 ( b ) , t h e b r i c k is s e v e r e l y d i s t o r t e d . Tables

3 a n d 4 list t h e c o m p u t e d

eigenvalues with the undistorted

and

distorted

bi-unit cube,

r e s p e c t i v e l y . T h e s e r e s u l t s o n l y i n c l u d e t h e e i g h t e e n n o n - z e r o e i g e n v a l u e s , i . e . , t h e six z e r o e i g e n v a l u e s c o r r e s p o n d i n g t o r i g i d b o d y m o d e s a r e n o t r e p o r t e d . L o c k i n g r e s p o n s e in t h e n e a r l y i n c o m p r e s s i b l e l i m i t c a n b e i n f e r r e d f r o m t h e n u m b e r o f t h e e i g e n v a l u e s w h i c h g o t o i n f i n i t y as K//x ~ co. F o r t h e p r o b l e m at h a n d , t h e s e c r i t i c a l e i g e n v a l u e s a r e h i g h l i g h t e d in T a b l e s 3 a n d 4. A l o c k i n g - f r e e e l e m e n t s h o u l d e x h i b i t o n l y o n e e i g e n v a l u e g o i n g t o i n f i n i t y as K//z---> oo, w h i c h c o r r e s p o n d s t o p u r e d i l a t a t i o n a l mode. For the undistorted element

i n s p e c t i o n o f T a b l e 3 first r e v e a l s t h a t a t o t a l o f s e v e n e i g e n v a l u e s g o t o

a

b

Fig. 2. Single element eigenvalue analysis. (a) Undistorted and (b) distorted bi-unit cube. Table 3 Eigenvalues for the undistorted bi-unit cube (K/# = 1 x 109): elements Q1, QM1 and Q1/E9 are integrated with the standard 8-point quadrature, while the elements O l / E l 2 and Q M 1 / E l 2 employ the special 9-point quadrature rule 01

0.33333 0.33333 0.10000 0.10000 0. I0000 0.13333 0.20000 0,20000 0.20000 0.20000 0.20000 0.11111 0.11111 0.ii111 0.66667 0.66667 0.66667 0.30000

x x x x

10 +° I0 +° 10 +' 10 +'

x I0 +'

10 +I 10 +' x 10 +, x 10 +1 x 10 +' x I0 +' x 10 ÷9 × i0 +') × 10 *9 × i0 +9 x i0 +~ x lO +~ x i0 +,o x x

OM1

01/E9

0.33333 x 10 +° 0.33333 x 10 +° 0.10000 x 10 +' 0. I0000 x 10 +' 0.10000 x I0 +' 0.13333 x I0 +' 0.20000 x 10*' 0.20000 x I0 +' 0.20000 x I0 + ' 0.20000 x 10 ÷' 0.20000 x 10 + '

0.33333 0.33333 0.66667 0.66667 0.66667 0.13333 0.20000 0.20000 0.20000 0.20000 0.20000 0.20000 0.20000 0.20000 0.11111 0.11111 0.11111 0.30000

0. I I I I i x 1 0

0.11111 0.11111 0.66667 0.66667 0.66667 0.30000

x

x × × × x

+9

10 ÷9 10 ÷° 10 ÷9 10 *~ 10 +9 I0 *'°

Ol/E12 x x x x x x x x x x x x x x x x x X

10 *° 10 *° 10 +° 10 *° 10 ÷° 10 "l 10 ÷' 10 +~ I0 +' 10 +' 10 +1 10 ÷' 10+' 10 ÷t 10 ÷9 10 *~ 10 ~9 10 +lO

0.33333 x 0.33333 x 0.66667 x 0.66667 x 0.66667 x 0.66667 x 0.66667 x 0.66667 x O. 13333 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.30000 X

OM1/E12 10 +° 10 +° 10 ÷° 10 +° 10 *° 10 +° 10 ÷° 10 *° 10 ' , 10 +1 10 +' I0 ÷' 10*' 10 *~ 10 +1 10 *L 10 +~ 10 +1°

0.33333 x 0.33333 x 0.66667 x 0.66667 x 0.66667 x 0.66667 x 0.66667 x 0.66667 x O. 13333 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.20000 x 0.30000 x

10 ÷° 10 +° 10 *o 10 *° 10 +° 10 +° 10 *° 10 +° 10 +, 10 +1 10 +l I0 +' 10 *l 10 ÷~ 10 *1 10 +' 10 +~ 10 +1°

J.C. Simo et al., Improved versions of assumed enhanced strain

379

Table 4 Eigenvalues for the distorted bi-unit cube (K/~ = 1 x 10~): elements Q1, QM1 and Q1/E9 are integrated with the standard 8-point quadrature, while the elements Q 1 / E l 2 and QM1/E12 employ the special 9-point quadrature rule Q1 0.35301 0.61467 0.93135 0.10137 0.12260 0.16827 0.17109 0.20302 0.21661 0.21751 0.13725 0.94870 0.10647 0.14769 0.59706 0.68230 0.75728 0.29855

QM1 x x x x x x x x x x x x

10 '° 10 ~° 10 +° 10 +t 10 ÷t 10+ ~ 10 'j 10 ÷~ 10 ,l 10 +t 10 +6 10 +~

x | 0 +9

x x x x ×

10 ~'~ 10 *'~ 10 +9 10 ~" 10" 1o

0.27845 x 0.39390 x 0.91662 x 0.99054 x 0.10415 x 0.12639 × 0.17118 x o. 17453 x 0.20457 x 0.22116 x 0.22257 x 0.94030 x 0.10434 x 0.14511 x 0.59036 × 0.71158 x 0.74401 x 0.30150 x

10 *° 10 +° 10 ~° 10 +t 10 *t 10 +t 10 .1 10 +i 10 *t 10 +l 10 ÷t 10 +a I 10 .9 [ 10 +'~ 10 ~9 10 *~ I0 +~ 10+ io

Q1 / E9

q l / El2

0.31724 x 10 ~° 0.44752 x 10 *° 0.59978 × 10 +° 0.64130 x 10 +° 0.10285 × 10 ÷1 0.11544 x 10 +t 0.15948 × 10 +t o. 16798 x 10 +z 0.17912 x 10 +~ 0.19402 x 10 ~l 0.20749 x 10 +l 0.22097 × 10 ÷t 0.24140 x 10 *~ 0.15130 x 01 +6 0.95983 x 10 +s 0.11073 x 10 +~ 0.15440 x 10 *~ 0.29832 x 10' io

0.35511 x 0.55361 x 0.62136 x 0.69599 × 0.70204 x 0.75893 x 0.11192 x o. 11678 x 0.17001 x 0.17569 x 0.19207 x 0.20323 x 0.22323 x 0.23151 × 0.23548 x 0.36086 x 0.69272 x 0.29812 ×

[ [

QM1 / E12 10 ÷° . 10 ~0 i0 +U 10 +° 10 +° 10 +° 10 +l 10" t 10 +1 10 ~J 10 +~ 10 +~ 10 ÷~ 10 *l 10 ÷j 10 ÷6 10 ÷7 I0 +to

0.24352 0.37683 0.51139 0.57520 0.61343 0.65762 0.71882 0.80084 0.12555 0.15870 0.17231 0.18349 0.18410 0.20342 0.21403 0.24072 0.26168 0.30075

x x x x x x x x x x x x x x x x x x

10 ÷0 10 ~° 10 +° 10 +° 10 ÷° 10 ÷° 10 ÷° 10 +° 10 +t 10 ~l 10 *1 10 *~ 10 .1 10 *~ 10 +t 10 ÷l 10 ~l 10 ÷Io

i n f i n i t y f o r t h e Q 1 d i s p l a c e m e n t e l e m e n t , in a g r e e m e n t w i t h t h e w e l l - k n o w n s e v e r e l o c k i n g r e s p o n s e to b e e x p e c t e d f r o m this e l e m e n t . A s e x p e c t e d , t h e p r o p o s e d m o d i f i c a t i o n o f t h e g r a d i e n t s o f t h e s h a p e f u n c t i o n s , l e a d i n g to t h e e l e m e n t l a b e l e d as Q M 1 , has n o e f f e c t in e i t h e r t h e n u m b e r o r m a g n i t u d e o f t h e e i g e n v a l u e s c o m p u t e d in u n d i s t o r t e d c o n f i g u r a t i o n s . T h e e n h a n c e d e l e m e n t Q 1 / E 9 e x h i b i t s f o u r e i g e n v a l u e s g o i n g to infinity, t h u s l e a d i n g to a f o r m u l a t i o n t h a t still e x h i b i t s a ' m i l d ' f o r m o f l o c k i n g . T h e a d d e d n i n e e n h a n c e d m o d e s in Q 1 / E 9 a r e o n l y a b l e to e l i m i n a t e t h e l o c k i n g a s s o c i a t e d w i t h t h e t h r e e b e n d i n g m o d e s o f t h e o r i g i n a l Q 1 , l e a d i n g to an e l e m e n t o p t i m a l in b e n d i n g . T h e r e m a i n i n g t h r e e v e r y l a r g e e i g e n v a l u e s a r e a s s o c i a t e d w i t h m o d e s a r i s i n g f r o m t h e t e r m ~t ~2~3 in t h e d i s p l a c e m e n t field. F r o m t h e r e s u l t s in T a b l e 3 w e o b s e r v e t h a t t h e s e t h r e e l a r g e e i g e n v a l u e s a r e n o l o n g e r p r e s e n t in t h e Q1/E12 and the QM1/E12 elements, which incorporate the three additional enhanced modes d e s c r i b e d in S e c t i o n 3.4. T h e o n l y r e m a i n i n g l a r g e e i g e n v a l u e in b o t h e l e m e n t s is a s s o c i a t e d w i t h p u r e v o l u m e t r i c r e s p o n s e . T h e n u m e r i c a l r e s u l t s also d e m o n s t r a t e t h a t t h e s e t w o e l e m e n t s a r e i d e n t i c a l in u n d i s t o r t e d c o n f i g u r a t i o n s , a r e s u l t to b e e x p e c t e d s i n c e t h e m o d i f i e d f o r m u l a ( 2 . 1 4 ) f o r t h e g r a d i e n t o f the isoparametric shape functions becomes exact for undistorted configurations. O n t h e o t h e r h a n d , i n s p e c t i o n o f t h e r e s u l t s in T a b l e 4 r e v e a l s s i g n i f i c a n t d i f f e r e n c e s f o r t h e d i s t o r t e d e l e m e n t . W e o b s e r v e t h a t t h e Q 1 e l e m e n t e x h i b i t s an a d d i t i o n a l l a r g e e i g e n v a l u e , w i t h o r d e r o f m a g n i t u d e l o w e r t h a n t h e o t h e r s e v e n v e r y l a r g e e i g e n v a l u e s . T h e e l e m e n t Q 1 / E 9 also e x h i b i t s a n e x t r a very large eigenvalue while the Q1/E12 element now exhibits two additional very large eigenvalues. T h i s p h e n o m e n o n is a d i r e c t c o n s e q u e n c e o f t h e e x a c t c o m p u t a t i o n o f t h e g r a d i e n t o f t h e i s o p a r a m e t r i c h o u r g l a s s m o d e s . A s w a s d i s c u s s e d in S e c t i o n 2.3, t h e c o n v e n t i o n a l c o m p u t a t i o n r e s u l t s in a d d i t i o n a l n o n - c o n s t a n t t e r m s in d i s t o r t e d e l e m e n t s t h a t c a u s e this ' m i l d f o r m ' o f l o c k i n g r e s p o n s e , In f a c t , t h e a p p e a r a n c e o f t w o e x t r a m o d e s in t h e Q 1 / E l 2 e l e m e n t i n s t e a d o f j u s t o n e as in t h e Q 1 a n d Q 1 / E 9 is a c o n s e q u e n c e o f t h e m o r e e x a c t i n t e g r a t i o n o f t h e s e n o n - c o n s t a n t t e r m s by t h e 9 - p o i n t q u a d r a t u r e r u l e e m p l o y e d in t h e f o r m e r e l e m e n t . T h e r e s u l t s f o r t h e Q M 1 a n d Q M 1 / E l 2 elements, which employ the m o d i f i e d f o r m u l a ( 2 . 1 4 ) , s u p p o r t this c o n c l u s i o n a n d s u b s t a n t i a t e t h e a r g u m e n t s p r e s e n t e d in S e c t i o n 2.3. T h e Q M 1 e l e m e n t n o l o n g e r e x h i b i t s t h e a d d i t i o n a l l a r g e e i g e n v a l u e p r e s e n t in t h e Q 1 e l e m e n t in distorted configurations, while the QM1/E12 element does not exhibit the two additional large e i g e n v a l u e s p r e s e n t in t h e Q 1 / E 1 2 e l e m e n t .

5.2. Upsetting o f a cylindrical billet T h e g o a l o f this t e s t p r o b l e m is to assess n u m e r i c a l l y t h e p e r f o r m a n c e o f t h e e n h a n c e d e l e m e n t s u n d e r i n v e s t i g a t i o n in h i g h l y s t r a i n e d s i t u a t i o n s a r i s i n g in c o m p r e s s i o n p r o b l e m s . T h e 3 D s i m u l a t i o n

J.C. Simo et al.. Improved versions of assumed enhanced strain

380

"~=0

fi=l

700

600 Q1----

500

QI/P0 - - - - Q1/E9 + 9 - p o i n t - - qM1/m2

400 fi=2

i~=3

Q

"~ ¢~

/ /

300 /

200

/

/

/

100 0

~ 1

2

I 3

I 4

5

Displacement, [ram]

fi=4

Fig. 3. Upsetting of a cylindrical billet. Initial and deformed

Fig. 4. Upsetting of a cylindrical billet. Load versus imposed

configurations at different imposed vertical displacements at the top.

top displacement curves obtained with different elements.

described below corresponds to the upsetting of a cylindrical billet and was suggested to us by Nagtegaal [25]. A perfectly axisymmetric specimen, with initial radius and height of 6.35 mm, is subjected to compressive loading obtained by prescribing the vertical displacement at the top face to a final value of 5.0 mm applied in 32 increments. Due to obvious symmetry considerations, only a quarter of the specimen need be considered. The initial finite element mesh, shown in Fig. 3, consists of a total of 384 brick elements. The contact condition that takes place at high compressive strains between the plane defined by the lower base and the lateral surface of the specimen is enforced via a conventional nodal penalty formulation, with a penalty parameter of 1 × 10 6. The finite deformation (uncoupled) elastoplastic constitutive model employed in the simulation is summarized in Box 1, with material properties recorded in Table 2. We note that the isochoric character of the plastic flow together with the very large plastic deformations result in a nearly incompressible response of the specimen. Figure 3 depicts the solution obtained with the proposed QM1/E12 element. The results demonstrate the good performance to be expected from this element in highly strained compressive situations. Figure 4 shows the load-displacement curves computed with different elements. Inspection of these results immediately reveals the significantly stiffer response exhibited by the Ol element and the slightly softer response of the new Q M 1 / E l 2 over the traditional Q1/P0 element. That the original enhanced element Q1/E9 with the standard 8-point quadrature formula is unable to reach a meaningful final

Q 1 / E 9 q--

QI/E9 +

8-point quadrature

g-polnt quadrature

Fig. 5. Upsetting of a cylindrical billet. Solution obtained with the Q I / E 9 element with the standard Gauss 8-point quadrature rule (6 = 3.25) and the same O l / E 9 enhanced element but integrated with the special 9-point quadrature rule (12 = 7.0).

J.C. Simo et al., Improved versions of assumed enhanced strain

381

solution in this type of problems is dramatically illustrated by the result depicted in Fig. 5(a). This disappointing performance of the Q 1 / E 9 element has been independently verified by Nagtegaal [25] and also observed by the authors in 2-D plane strain and axisymmetric problems in the presence of high compressive strains. We emphasize that the pathology illustrated in Fig. 5(a) should be attributed to under-integration resulting from the conventional 8-point quadrature formula and not to the mild form of locking exhibited by this element. In fact, use of a higher order quadrature rule (e.g., 3 x 3 × 3 Gaussian quadrature) completely circumvents this pathology while preserving the excellent response of the element in localization problems. This point is demonstrated in Fig. 5(b) where the final solution obtained with the same Q 1 / E 9 element, but with the 8-point formula replaced by the proposed 9-point quadrature rule, no longer exhibit the pathology present in Fig. 5(a). Figure 4 includes the loaddisplacement curve for this case and shows a response of this element very close to the QM1/E12 element for this problem. Computationaily, the remarkable result is that the proposed 9-point quadrature formula produces results essentially identical for all practical purposes to the far more expensive 27-point quadrature formula. The performance of this higher order quadrature in plane strain localization problems is investigated in the next section. 5.3. Plane strain localization p r o b l e m

The 2-D counterpart of the enhanced Q1/E9 element has 4 enhanced modes (Q1/E4), reduces to the QM6 element of [21] for the linearized problem and does not exhibit locking in the incompressible limit. However, as pointed out above, this element is also underintegrated by the standard 4-point quadrature leading to severe hourglassing in certain highly strained situations. The replacement of the conventional 4-point formula by the 5-point quadrature described in Section 3.5 provides a computationally effective procedure for circumventing these difficulties, without resorting to the more expensive 9-point quadrature formula (in 2-D, or 27-point formula in 3-D).

EQ. PI AST. STN. ¢ < 4.000E-01 ¢

F¢ ¢. ¢ ¢

> 8.000E-01

a

b

Fig. 6. Plane strain tension test. (a) Initial configuration (200 element mesh). (b) Deformed configuration at an imposed top displacement of 4.24, computed with the enhanced formulation (4 enhanced modes, Q1/E4) and the special 5-point quadrature rule.

382

J.C. Simo et al.. Improved versions of assumed enhanced strain

The goal of this problem is to investigate how the good performance exhibited by the Q 1 / E 4 element with the conventional 4-point quadrature formula is affected by higher order quadrature rules. For this purpose, the same plane strain tension test problem used in [13] is considered here, again within the framework of finite strain (uncoupled) plasticity. The specimen under investigation has a width of 12.826 mm and a length of 53.334 mm. The initial configuration is depicted in Fig. 6(a), together with the imposed boundary conditions. A quarter of the specimen is considered along with a finite element mesh consisting of 200 enhanced quads. In order to trigger the formation of the shear bands, which in this case appear at 45 ° of the load axis, a mechanical hardening-softening law is introduced together with a small geometric imperfection (the width at the center is reduced by 0.982%). The saturation hardening law in Box 1 is then replaced in this simulation by the following double-parabola hardeningsoftening law used in [26]: K , , ( e ) = O'o[C + B e - A e 2 ] ,

where A=a,

B=b,

C=1,

f o r e ~<6,

A=a/d,

B=b/d,

C=c,

fore1>6.

The following numerical values are used: o-o = 0.45 GPa, a = 14.5, b = 5.8, c = 7.41875, d --- 5.0625 and 6 =0.2. Figure 6(b) shows the solution obtained with the enhanced element and the 5-point quadrature rule along with the distribution of the equivalent plastic strain a at an imposed top displacement of 4.24 mm. In Fig. 7 this solution is compared with the one obtained employing the enhanced element with the standard 4-point quadrature and the Q1/P0 element. The load-displacement curves obtained with this elements are shown in Fig. 8. From these results, we conclude that although the 5-point quadrature formula slightly increases the stiffness of the original element, the new element retains essentially the same good performance in problems involving localization obtained with the conventional 4-point formula. The vastly superior performance of the element with the modified 5-point quadrature rule over that exhibited by the mean dilation approach should be noted.

IIIlUlIIIIIIIIIIII IIIIHIIIIIIIII IIIIHIIIIIIIII IIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIII IIIIHIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII

IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII

i i i i{i i i i!i.

IIIIIIIIIIIIIIIIlUl IIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIII IIIIIlUlIIIIIIIIII IIIIIIIIIIIIIIIIIII lllllllllllllllllll illllllllllllllllll i111111111111111111

Q1/E4 -+4-point quadrature

IIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIII IllllIIIIIIIIIlllll IIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIII IIIIHIIIIIIII] IIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIII [1111111111111111111 iii{!!i!~!!iiiii~i~i iiiiiiiiiiiiiiiiiiii iiiiiiiiiiiiiiiii~i

IIIIIIIIIIIIIIIIIIII IIIIIIIIIIlUlIIIIII IIIIIIIIIIlUlIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIIIII

Q1/E4 + 5-point quadrature

IIIIIIIIIIIIIII]1111 tllllllllllUlllllll IIIIIIIIIIIIIIIIIII IIIIIIIIIIMIII IIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIII IIIIIIIIIIIIIIIIII IIIlUlUlIIIIIIII IIIIIIHIIIIII IIIIIMIIIIIII

Q1/PO

Fig. 7. Plane strain localization problem. Solutions at an imposed top displacement u = 4.24 mm obtained via the enhanced element (QI/E4) with 4 or 5 point quadrature rule and the mixed element Ol/P0.

J.C. Simo et al., Improved versions of assumed enhanced strain

383

10

8

=-

6

Q

Q1/P0-.5-pointq u a d r .

Q1/E4 ---

4 - p o i n t quach'. Q 1 / E 4 i i i 1 2 3

: 4

Displacement, [ram] Fig. 8. Plane strain localizationproblem. Load versus imposed top displacementcurves obtained with differentelements. 5.4. T h e r m o m e c h a n i c a l n e c k i n g p r o b l e m

This final test problem is chosen to highlight the undesirable effects associated with the form of mild locking observed in the Q1/E9 element in the presence of finite deformations. We consider the full coupled thermomechanical problem associated with the necking of a cylindrical bar in simple tension. Both the physical coupling (arising from thermal stresses, self-heating induced by plastic dissipation, elastoplastic structural heating as well as thermal softening) and the geometrical coupling induced by finite deformations in the specimen are accounted for in the simulation. The integration algorithm employed for this problem is the unconditionally stable fractional step method first proposed in [27], which consists of a mechanical phase at constant entropy followed by a solution of the heat equation. The property of unconditional stability can be shown to hold even in the full non-linear regime of interest in this investigation. Moreover, the algorithm retains the key computational convenience of the conventional conditionally stable staggered schemes for thermomechanical problems; i.e., modularity of the implementation and symmetry of the sub-problems. It should be noted that the enhanced high-performance elements described in this work are used only in the solution of the mechanical (adiabatic) sub-problem. The conventional displacement formulation is used for the thermal subproblem. The bar has a total length of 53.34 mm and a radius of 6.35 mm. Due to the symmetry of the problem only one-eighth of the specimen is discretized and subjected to the obvious symmetry boundary conditions. Figure 9(a) shows the finite element mesh which consists of a total of 120 elements. As in the previous examples, the loading is applied by prescribing the displacement at the top. We carry the simulation until a top imposed displacement of 7.0 mm is achieved in a total number of 25 increments. The iterative scheme for the static condensation of the enhanced modes described in Section 4.5 is employed in this example. The bar is assumed thermally insulated on the lateral face and the temperature is assumed fixed to the reference value at the top face. No material or geometric imperfection is introduced to trigger the necking of the bar. We note again that the isochoric nature of the plastic flow effectively enforces the incompressibility constraint. Elements exhibiting severe locking, such as the conventional Q1 displacement element, produce simulation results in which not even mild necking of the specimen can be detected. Figures 9(b) and 9(c) respectively show the distribution of the equivalent plastic strain and temperature in the final deformed configuration computed with the QM1/E12 element. The solutions obtained with QM1/E12 and the original enhanced element Q1/E9 are compared in Fig. 10. The Q1/E9 element was shown (see Section 5.1) to exhibit a mild form of locking, the consequences of which are clearly demonstrated in Fig. 10. While the Q1/E9 element is able to capture the necking process in the bar, the simulation results exhibit a severe global hourglassing pattern. We remark that the use of a higher order quadrature rule does not improve the response of the element in this test problem (contrary to the situation found in Section 5.2; see Fig. 5). The cause of this hourglassing mode is to be found in the extra locking modes described in Section 5.1 which appear in 3-D and are removed via an additional enhancement in the QM1/E12 brick.

384

J.C. Simo et al.o Improved versions of assumed enhanced strain

/\

/

\

/ /

/ /

/

/ /

/

/ \

/ /

/ \

/ /

/ J J

....,,..

J

a

TEMPERATURE

EQ. PLAST. ST. Min = 0.00E+00 Max = 1.23E+00

Min- 0.00E+00 M a x - 7.23E+01

j 28.0E-01

1.00E+01 2.00E+01 3.00E+01 4,00E+01 5.00E+01 6.00E+01

~.ooE-o~

4.60E-01 6.40E-01 8.20E-01 1.00E+00

b

c

Fig. 9. Thermomechanical necking of a cylindrical bar. (a) Initial configuration. Distributions on the deformed configuration of the (b) equivalent plastic strain and (c) temperature at an imposed top displacement u = 7 (computed with the Q M 1 / E l 2 brick).

6. Concluding remarks A general framework for the construction of enhanced interpolations of the deformation gradients has been exploited in the design of low order finite elements which exhibit improved performance in the finite deformation regime over the conventional Galerkin finite element method. A main goal of this

J.C. Simo et al., Improved versions o f assumed enhanced strain

385

Q1/E9 QM1/E12 Fig. I0. Thermomechanical necking of a cylindrical bar. Solutions at an imposed top displacement u = 7 obtained with the QI/E9 and QMI/EI2 elements. work has been the design of a low-order 3-D brick that generalizes to the finite deformation regime incompatible mode elements based on Wilson's brick. The approach circumvents the locking response experienced by this element in the nearly incompressible limit while preserving its good performance in bending dominated and localization problems. The final element, however, can no longer be interpreted within the framework of a method of incompatible modes, even in the linear regime. It has been shown that the conventional 8-point quadrature formula in 3-D under-integrates both the original Wilson brick and the proposed enhanced element. The effects of this under-integration can be disastrous in practical calculations involving high compressive strains. These undesirable effects can be circumvented by a 9-point quadrature formula which produces results nearly identical to those obtained with the far more expensive 27-point formula. Finally, the mild form of locking arising in distorted elements can be bypassed by introducing a simple modification in the computation of the shape function derivatives which has no effect in undistorted configurations. A key advantage of the proposed methodology arises in the actual implementation of inelastic constitutive models, since the conventional 'strain-driven format' of current stress-point integration algorithms is preserved. We conclude by noting that the proposed methodology can be exploited in the design of sound one-point quadrature elements free from adjustable parameters. This topic will be pursued in a separate publication.

Acknowledgment We thank G. Krishnakumar for his comments and input on the subject of this manuscript. Support for this research was provided by N C E L (Naval Civil Engineering Laboratory), Port H u e n e m e , and C E N T R I C Engineering Systems Inc, Palo Alto, CA. This support is gratefully acknowledged.

References [1] L.R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials, Math. Modelling Numer. Anal. 9 (1985) 11-43. [21 F. Brezzi and M. Fortin, Mixed and Hibrid Finite Element Methods (Springer, Berlin, 1991). [3] J.C. Nagtegaal, D.M. Parks and J.R. Rice, On numerically accurate finite element solutions in the fully plastic range, Comput. Methods Appl. Mech. Engrg. 4 (1974) 153-177. [4] M. Ortiz, Y. Leroy and A. Needleman, A finite element method for localized failure analysis, Comput. Methods Appl. Mech. Engrg. 61 (1988) 189-214.

386

J.C. Simo et al., Improved versions o f assumed enhanced strain

[5[ T. Belytschko, J. Fish and B.E. Engelman, A finite element with embedded localization zones, Comput. Methods Appl. Mech. Engrg. 70 (1988) 59-89. [6[ J.C. Simo and T.J.R. Hughes, On the variational foundations of assumed strain methods, J. Appl. Mech. 53 (1986) 51-54. [7[ J.C. Simo and M.S. Rifai, A class of assumed strain methods and the method of incompatible modes, Internat. J. Numer. Methods Engrg. 29 (1990) 1595-1638. [8] E.L. Wilson. R.L. Taylor, W.P. Doherty and J. Ghaboussi, Incompatible displacement models, in: S.J. Fenves, N. Perrone, A.R. Robinson and W.C. Schnobrich, eds., Numerical and Computer Models in Structural Mechanics (Academic Press, New York. 1973). [9] G. Strang and G. Fix, An Analysis of the Finite Element Method (Prentice Hall, Englewood Cliffs, NJ, 1973). [10[ R.L. Taylor, J.C. Simo, O.C. Zienkiewicz and A.C. Chan, The patch test: A condition for assessing finite element convergence, Internat. J. Numer. Methods Engrg. 22 (1986) 39-62. [I 1] D. Kosloff and A. Frazier, Treatment of hourglass patterns in low order finite element codes, Internat. J. Numer. Methods Geomech. 2 (1978) 57-72. [12] U. Andelfinger, E. Ramm and D. Roehl, 2D- and 3D-enhanced assumed strain elements and their application in plasticity, in: Onate et al., eds., Proc. Complas IV, Barcelona, 1992. [13] J.C. Simo and F. Armero, Geometrically nonlinear enhanced strain mixed methods and the method of incompatible modes, Internat. J. Numer. Methods Engrg. 33 (1992) 1413-1449. [14] T. Belytschko, J.S.-J. Ong, W.K. Liu and J.M. Kennedy, Hourglass control in linear and nonlinear problems, Comput. Methods Appl. Mech. Engrg. 43 (1984) 251-276. 1151 P.G. Ciarlet, Mathematical Elasticity. Volume 1: Three-dimensional Elasticity (North-Holland, Amsterdam, 1988). [16] T.J.R. Hughes, The Finite Element Method (Prentice Hall, Englewood Cliffs, NJ, 1988). [17] O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, 4th Edition, Vol. 1 (McGraw-Hill, London, 1989). [18] T. Belytschko and L.P. Bindeman, Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for nonlinear problems, Comput. Methods Appl. Mech. Engrg. 88 (1991) 311-340. [19] B.C. Koh and N. Kikuchi, New improved hourglass control for bilinear and trilinear elements in anisotropic linear elasticity, Comput. Methods Appl. Mech. Engrg. 65 (1987) 1-46. [20] M. Mallet, C. Poirier and F. Shakib, A new finite element formulation for computational fluid dynamics: Development of an hourglass control operator for multidimensional advective-diffusive systems, Comput. Methods Appl. Mech. Engrg. 94 (1992) 429-442. [21] R.L. Taylor, P.J. Beresford and E.L. Wilson, A non-conforming element for stress analysis, lnternat. J. Numer. Methods Engrg. 10 (1976) 1211-1219. [22] B.D. Reddy and J.C. Simo. Stability and convergence of a class of enhanced strain methods, F R D / U C T Center for Research in Computational and Applied Mechanics, Report No. 174, University of Cape Town, South Africa, 1992. [23] P.G. Ciarlet, The Finite Element Method for Elliptic Problems (North-Holland, Amsterdam, 1978). [24] J.C. Simo, R.L. Taylor and K.S. Pister, Variational and projection methods for the volume constraint in finite deformation elasto-plasticity, Comput. Methods Appl. Mech. Engrg. 51 (1985) 177-208. [25] J.C. Nagtegaal, Private communication, 1992. [26] A. Nacar, A. Needleman and M. Ortiz, A finite element method for analyzing localization in rate dependent solids at finite strains, Comput. Methods Appl. Mech. Engrg. 73 (1989) 235-258. [27] F. Armero and J.C. Simo, A-priori stability estimates and unconditionally product formula algorithms for nonlinear coupled thermoplasticity, SUDAM Report no. 92-4, 1991, Stanford University; Internat. J. Plasticity, in press. [28] J.C. Simo and C. Miehe, Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation, Comput. Methods Appl. Mech. Engrg. 98 (1992) 41-104.