INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL.39,2337-2361 (1996) HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS SHMUEL L. WEISSMAN Symplectic Engineering Corporation, 1350 Solano Ave. #26, Albany, CA 94706, U.S.A. SUMMARY This paper presents a family of low-order high-accuracy three-dimensional brick elements.The elements are formulated via a three-field variational principle. The assumed (independent)strain field is constructed from two disjointed distributions. The first contains the lower-order distribution and its dimension is the minimum required to satisfy stability requirements. An energy constraint, which is enforced weakly at the element level, is used to relate the second distribution to the first. The stress field is chosen to a priori satisfy a similar energy constraint. As a result, internal constraints (e.g. incompressibility) are automatically satisfied by these fields, and locking behaviour is avoided. A J,-plasticity model illustrates the proposed elements' performance in nonlinear solids. The excellent performance of the proposed elements is demonstrated with numerous challenging examples, including many that are usually modelled by shell elements. KEY WORDS elastoplastic; high-accuracy; locking-free;low-order; nonlinear; three-field functional 1. INTRODUCTION This work presents a family of low-order high-accuracy eight-node brick elements that are suitable for elastoplastic analysis (including the case of perfect plasticity). The importance of this work stems from the following: 1. The proposed elements exhibit excellent performance even when applied to thin shell problems; therefore, they offer a viable alternative to shell elements especially in problems involving non-smooth shells. 2. None of the proposed elements locks at the nearly incompressible limit. This stands in sharp contrast to other recently proposed brick elements. For example, the nine- and fifteen-parameter elements of Andelfinger et al.' and the nine-parameter element of Simo and Armero' possess four eigenvalues that go to infinity as the incompressible limit is approached (five when the elements are distorted). A locking-free element should exhibit only one eigenvalue going to infinity as the nearly incompressible limit is approached. The development of low-order high-accuracy finite elements has traditionally been one of the main goals of finite element research, mainly because the computational effort associated with element array processing increases linearly with the number of elements. The effort associated with the solution of the global system of equations is proportional to the number of equations to the third power. Therefore, a significant reduction in the computational effort can be achieved if the number of global degrees-of-freedom is reduced, even at the expense of a somewhat larger computational'effort at the element level. CCC 0029-59811961142337-25 0 1996 by John Wiley & Sons, Ltd. Received 28 April 1995 Revised 15 November 1995 2338 S. L. WEISSMAN The design of low-order high-accuracy finite elements is based on mixed formulations. Two approaches to mixed finite elements are commonly used: (1) a two-field (Hellinger-Reissner) functional where the displacement and stress fields are the assumed independent variables (the works of Pian and Tong3 and Punch and Atluri4 are representative of this approach); and (2) a three-field (Hu-Washizu) functional, where the displacement, stress and strain fields are the assumed independent variables. Representative of this approach are the works of Nagtegaal et ~ l . Hughes,6 , ~ and Simo et aL7 The works of Malkus and Hughesa and Simo and Rifai’ have shown that mixed formulations encompass, as particular cases, reduced and selective reduced integration schemes (e.g. Reference lo), and the method of incompatible modes (e.g. Reference 11). The design of mixed elements requires special care in the selection of the assumed fields so that stability requirements are observed (e.g. References 12 and 13). Moreover, the use of mixed finite elements does not guarantee high-accuracy elements, a fact that was cast into a limitation principle by de Veubeke.I4 From a practical standpoint, it is easier to design the assumed stress field than the assumed strain field to satisfy the above restrictions. Additionally, elements derived from the two-field formulation are more computationally efficient. Consequently,elements based on the Hellinger-Reissner functional are more popular. However, the Hellinger-Reissner functional is based on the complementary energy function, which for general constitutive laws is not readily ascertainable. The Hu-Washizu functional (see e.g. Reference 15) is formulated in terms of the strain energy function, thus making it more attractive for non-linear materials. For example, integration schemes for elastoplasticity, such as the radial return algorithm for J2-plasticity of Wilkins,16 are founded on the premise that the elastoplastic problem is strain driven. Consequently, in elements based on the Hu-Washizu variational principle, these integration algorithms are formulated entirely in terms of the assumed total and plastic strain fields, and testing for plastic loading can be performed independently at each material (integration) point. Furthermore, elements based on the Hellinger-Reissner functional are not suitable for the modelling of perfect plasticity (see e.g. Reference 17). For these reasons, the more general Hu-Washizu functional is used in this work. An outline of the paper follows. Section 2 presents the strong (or local) form of the elastoplastic boundary value problem (the Appendix contains a summary of the radial return mapping used in the numerical simulations).Section 3 presents the weak form (includingthe construction procedure for the assumed stress and strain fields), derived from a three-field functional. Section 4 provides the assumed fields used to generate the elements proposed herein. In Section 5, numerical results showing the excellent performance of the proposed elements are tabulated. Finally, Section 6 presents concluding remarks. 2. THE BOUNDARY VALUE PROBLEM: STRONG FORM This section presents the strong, or local, form of the boundary value problem addressed in this work. The constitutive law is the rate-independent elastoplastic model, J2-plasticity with linear isotropic and kinematic hardening, used in the examples provided in Section 5. 2.1. Kinematics Let B c R3 denote a bounded open body whose closure has a (piece-wise)smooth boundary, aB; where aB is assumed to have a (piece-wise)continuous with outward unit normal, n, and to possess the following structure: a,BnaUB= 8 and a&%udvB = aB where avBis the part of dB on which displacements are prescribed as UIaug= 0 (given) (1) HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2339 and 8 9 is the part of 8B on which traction is specified as t = o*nla3 (given) (2) In the above equations, U(X) is the displacementfield associated with particles X c 1. Assuming small strains, the (total) strain field, E, is related to the displacement field, U,by E:= V'U (3) where V s is the symmetric part of the gradient operator. In component form, the total strain is given by Eij = +(Vi,j + Vj,i) (4) Also, the usual assumption that E:= where E' and E~ Ee + EP (5) are the elastic and plastic strain tensors, respectively, is adopted. 2.2. Balance of momentum Let b : g --* R3and p : W + R be the (given) body force and mass density, respectively. The local form of the balance of momentum (equilibrium)equations (neglectinginertia terms) are given by:* in (divo + pb = 0 dT = d 2.3. Constitution A constitutive law that relates the stress tensor, a(X), to the strain tensor, E(X),is necessary to complete the statement of the boundary value problem. To this end, a stored energy function W (X,E(X)- E ~ ( X ) x) :Y~ + R is introduced so that the second rank stress tensor (CE 9)is given by:' o(X):= a, W(X,E(X) - EP(X)) (7) and Y is the space of second rank symmetric tensors. The elastic tensor, C (a rank four tensor), is defined by c:=aEEw (8) and possesses the following symmetries: Cijkl= C k l i j = C i j l k = C j i l k . The strain energy function W is assumed to be quadratic so that W(Ec):=+EC:C:EC (9) and C is given in terms of the Lame coefficients, A and p, by c = 118 1 + 2pI The symmetry of the stress tensor is a result of the balance of angular momentum If only initially isotropic and homogeneous materials are considered (as in this work), the explicit dependence of the stress tensor on X can be dropped" 2340 S. L. WEISSMAN where 1 and I are the rank two and rank four identity tensors, respectively. The Lame coefficients are expressed in terms of Young’s modulus, E, and Poisson’s ratio, v, by I= (1 + vE v)(l - 2v) and p=- E 2(1 v) + To complete the constitutive law, the evolution of the plastic strain and hardening parameters , x R“‘+ R be the yield function,: where rn is the dimension must be described. Letf(a(8, E ~ ) q):Y of the hardening parameters, vector q. An admissible state is such that f (0, q) d 0 (12) The surfacef(a, q) = 0 is known as the ‘yield surface’ in stress space. The evolution of the plastic strain and hardening parameters are referred to as the ‘flow rule’ and ‘hardening law’, respectively. For the case of associative J,-plasticity, used in this work, the yield function, flow rule, and hardening law are given by where, following standard notations, the hardening vector, q, is taken as q:= { a , ij} where a is the ‘equivalent plastic strain’ that defines ‘isotropic hardening’ and q defines the centre or ‘kinematic hardening’ of the von Mises yield surface (see, e.g. Reference 19); and q:= dev[a] - q, tr[q]:= 0, and H’(a) and K ( a )are the kinematic and isotropic hardening moduli, respectively. In this work, the following linear isotropic and kinematic hardening laws are adopted for H’(a) and K(a): + H’(a):= (1 - B)Ra and K(a):= ay BRa with BE[O, 13 (14) Finally, i, is referred to as the ‘consistency parameter’, and is assumed to obey the following Kuhn-Tacker complementary conditions: I i, 2 O,f(a, q) d 0, and i,f(a, q) = 0 I In addition, f is constrained to satisfy the ‘consistency requirement’: The Appendix provides an outline of the radial return mapping used to update the plastic strain and hardening parameters along with the form of the algorithmic tangent. For further details regarding these issues, see e.g. Reference 19. 3. THE WEAK FORM AND FINITE ELEMENT APPROXIMATION This section presents a mixed finite element approximation to the elastoplastic boundary value problem considered in Section 2. The formulation is based on a three-field (Hu-Washizu) ‘The explicit dependence of the stress on the total and plastic strains is included as a reminder that the total and plastic strain tensors and hardening parameters are the independent variables in this development HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2341 functional. For highly accurate elements in coarse meshes the approach advocated in Reference 20 and exploited in the context of two-dimensional elements in Reference 21 is adopted. Accordingly, the assumed stress and strain fields are enlarged (over the minimum assumed fields that satisfy stability conditions) and constrained so that the additional terms do not contribute (weakly) to the internal energy, thus obtaining additional freedom to satisfy internal constraints. 3.1. Boundary value problem: weak form Let the total free energy be given by the following three field functional: I ~ ( o E,, U):= I [ W ( E- 9 )+ 0 : {V’U - E}] dY - nExT(U) (17) where the independent variables are the stress, strain and displacement fields. The plastic strain, cP, is given by rate equations such as (13)2 (i.e. E~ is viewed as an internal variable; see e.g. (U) is the potential energy of the external loading, which under the Reference 22), and nExT assumption of dead loading is given by The weak form of the boundary value problem is obtained by looking at stationary points of the functional (note that this is a constrained optimization problem), and is given by Dn = ss + r:[V”U - E ] ~ Y ( : [ a , W - o]dY j9 M kinematics + I V’q:adV - D&xT.q = 0 M constitution (19) V balance of momentum where r, ( and q are arbitrary admissible variations of the stress, strain and displacement fields, respectively (note that q is required to satisfy the displacement boundary conditions). Since r, 5, and q are arbitrary and independent of each other, it follows that equation (19) is composed of three independent parts as indicated above. Furthermore, by the fundamental lemma of calculus, the strong form of the strain-displacement,constitutive, and (after integration by parts) balance of momentum equations are recovered. 3.2. Finite element approximation In this work, a finite element approximation in which the assumed stress and strain fields are approximated independently in each element, and thus are discontinuous across element boundaries, is presented.” A typical element’s domain .@:’ c a is obtained by a discretization of the body, B, into finite dimensional regions such that a z up: .@:I, where ‘nel’ is the number of elements in the mesh. Consequently, the total free energy is approximated by nel el $By equation (19), only C-,continuity across the element boundary of the assumed stress and strain fields is strictly required 2342 S. L. WEISSMAN nr where (&, cC1, Ucl)is given by equation (17) with all fields and domain understood as element variables and domain, respectively.The remainder of this work is carried out at the element level. Consequently, the superscript 'el' denoting the element level can be omitted without any confusion. Finite element approximations can be constructed directly from equation (20) by introducing explicit approximations to the assumed fields. Unfortunately, selecting the assumed fields to achieve high accuracy in coarse meshes is not a trivial task in the presence of internal constraints. To overcome this difficulty, the approach proposed by Weissman" is followed. Accordingly, the assumed strain distribution is enlarged, over the minimum required for stability, and is constrained so that the added terms do not contribute to an energy-like quantity, consequently providing the additional freedom required to satisfy internal constraints (e.g. incompressibility). This constraint is appended to the free energy functional as follows:7 n(a,E, U):= J', [ W ( E- EP) + cr:{VSU- + 3,:&W]dV - IIExT(U) E} (21) Let the assumed stress, strain and displacement fields be approximated by a:= ss E:= E~ (22) + c2 = Elel + Ezez (23) U:= Nd (24) and respectively, E~ is the added strain distribution. In the above equations, the assumed stress and strain interpolations are constructed so that S and El contain the minimum distribution required for stability, E l n E z = 8, and the e's and s's are the vectors of stress and strain unknowns, respectively; the standard isoparametric interpolation is used for the displacement field where N is the interpolation matrix and d is the vector of nodal displacements. Also let the Lagrange multipliers be given by A:= E'z (25) where E' is the interpolation matrix and z is a vector of independent parameters." The added Euler-Lagrange equations due to 3, are used to relate Aez (the increment in ez) to Ael by Aezn= - (laEiTC:'gEzdV)-l[ + jaEiT&W,,dV EiTC:lgE1dVAel, Recall that the functional in equation (21) represents the approximation to the true functional within the element's domain, and all fields are finite dimensional. Consequently, the added constraint serves as an orthogonality restriction on the acceptable distribution, rather than implying that d, W = 0. An additional constraint can be added explicitly in which d, W is replaced by S (see Reference 20). This is not explicitly included in the present work because S and 1. are selected to be orthogonal and thus, this added constraint is automatically satisfied "Notethat dim(E') = dim(E,) and that the added constraint replaces a natural Euler-Lagrange equation due to the added strain terms. This choice will be explained in Section 4 below 2343 HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS respectively, where the subscript 'n' denotes the iteration number within the Newton iteration used in solving the nonlinear problem, and C;lg is the algorithmic elastoplastic tangent (see the Appendix). Let the following definitions be introduced: STIElel, r2,:= + E2ezn+ E2V, - Bd,] dV and r3,,:= J apbdV +J t dT - GTsln (32) where B is given by B:= V"N. (33) Let S be chosen so that it is orthogonal to A,** the Lagrangian multiplier, and note that as a result of equation (27) En too is orthogonal to A so that The remaining (linearized) Euler-Lagrange equations can be cast into a matrix form as follows: Eliminating the stress and strain coefficients yields the standard form: GT(A,H;'AT)-'GAd, L--y----J = r3, + GT(A.H;lA;f)-l(H,,A~lrl. + rZn) L " (36) J R" K. In the present work, A,, is constructed to be a fully ranked square matrix so that the stiffness matrix, K,, and residual, R,, are given by K,:= GTAiTH,A;'G and R,:= r3, + GTA,T(rl,,+ H,,Al1r2,,) (37) respectively. For further details, including analysis of algorithmic stability, see Reference 20. **This choice of S is motivated by numerical efficiency 2344 S. L. WEISSMAN 4. ASSUMED FIELDS SELECTION AND PROPOSED ELEMENTS This section presents the assumed independent fields required for the generation of a family of three-dimensional eight-node elements. These elements are constructed following the procedure outlined in Section 3. Section 5 will show that these elements achieve the desired high accuracy in coarse meshes while satisfying stability requirements. 4. I . Displacement interpolation The displacement field is interpolated using the standard isoparametric shape functions (e.g. Reference 23). Consequently, the displacement field is approximated by 8 where dI and NI are the vector of nodal displacement and shape function associated with node I, respectively, and NIis given as a function of the element’s natural coordinates by + &O(l+ qrq)(l + <rO NAt, q, (39) 4.2. Stress interpolation Following Pian and S~mihara,’~ the assumed stress field is first constructed in the element’s natural coordinates (5, q, <) as (40) [9ii 9 2 2 y 3 3 9 1 2 9 2 3 9 3 i I T ’ % where r 1 0 0 Y=l 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 q 0 0 0 0 0 ~ q < 0 0 0 0 0 0 0 0 0 ~ ~ ~ ~ 0 0 0 0 0 ~ q ~ 0 0 0 0 0 0 0 0 < 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 q 0 0 0 11 1 0 0 0 0 ~ q 0 0 0 0 0 0 (41) 0 is the vector of stress parameters associated with 9and thus its dimension is 18 (i.e. it contains 18 parameters).The assumed stress field is transformed from the (5, q, [) space into the (xl, x2, x3) space by means of the following transformation: aij = YiIRfJXIJ (42) where the indices i a n d j take the values xl, x2 and x3; and the indices I and J take the values 5 q I; and the ‘deformation gradient’ YiIis given by Yx,e = axl/a5, etc. To preserve the ability to pass the constant strain patch test, the transformation (42) is based on values .Fir at the element’s center.24Consequently, the transformation (in matrix notation) is given by S = Fox (43) 0 2345 HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS where 4.3. Strain interpolation The stress interpolation in the element's natural coordinates, equation (41), is also used as the strain g1interpolation (i.e. B1= Y ) ,while g2is given by < < q e g o o 0 0 0 0 0 0 q q c g q o o 12 = 0 0 0 0 0 0 0 0 i 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 V 9 parameter 0 0 rcrq1: 0 0 0 0 0 0 0 0 0 J 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - g q o o o o o o q ~ o o o o o o < q 15 parameter The assumed strain interpolation in the (xl, x2,x3) space is obtained by using the inverse of equation (43), and is given (in matrix form) by E = F i T 1= FiT[&l 8 2 1 where C is the assumed strain interpolation in the element's natural coordinates. (46) 2346 S. L. WEISSMAN 4.4. Lagrange multipliers interpolation Analysis of equation (26) shows that stability is obtained for E’= E2 (47) (see Reference 20). This choice is also in agreement with the requirement that the enhancement to the assumed fields does not contribute to the internal energy. Moreover, it naturally gives rise to an Euler-Lagrange equation identical to the one due to the constraint introduced in equation (21). However, this choice fails to satisfy equation (34), for elements with a non-constant transformation Jacobian J = det(9). For this reason, the following modified interpolation is used for the Lagrange multipliers:tt 4.5. Proposed elements The assumed fields presented above are used to generate three-dimensional eight-node brick elements. Three elements are presented that differ by the enhancement’s dimension. The three options are 9, 15 and 24 parameters, as shown in equation (45), and accordingly are denoted by B8-9P (8-node brick element with 9 parameters), B8-15P and B8-24P, respectively. The motivation for these selections is as follows: To overcome locking behaviour only the normal stresses/ strains must be enhanced; thus the selection of nine-parameters’ enhancement, which provides the same approximation in the normal strains in all directions. In this case, however, the shear terms are approximated by an incomplete linear distribution (in the element’s natural coordinates). A better approximation is obtained if a complete linear distribution is adopted (hence the 15-parameter element). The third element is constructed to have the same distributions for both the normal and shear terms, and thus the choice of the 24-parameter element. Finally, the stress field can be used to formulate an element based on the Hellinger-Reissner formulation. This element is identical to that proposed by Pian and Tong.3 4.6. Remarks on implementation The procedure presented above is best suited for conveying the proposed ideas. From an implementational standpoint, however, it is best to operate using the sparse matrices in the element’s natural co-ordinates. Thus, for efficiency, the assumed field should not be pushed forward to the physical space; rather, the tangent matrix and B-matrix should be pulled back to the element’s natural co-ordinates (e.g B = FtB). As will be shown in the next section, the last enhancement (i.e. extension to 24 parameters) does not significantly improve the results over the B8-15P elements and therefore, should not be used in practice. It is included precisely to demonstrate that adding modes does not necessarily improve element performance. ++Notethat in the case of constant Jacobian elements the Lagrangian multiplier’s interpolation provided by equation (48) is equivalent to that of equation (47) HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2347 5. NUMERICAL EXAMPLES The performance of the proposed elements is now evaluated with a number of challenging problems selected from the literature. First, in Section 5.1, the elements are shown to pass the constant strain patch test. Second, in Section 5.2, the eigenvalues of a single element's stiffness matrix are examined to demonstrate that the proposed elements do not lock at the nearly incompressible limit. Third, the higher-order patch test is used to evaluate the proposed elements' performance in bending-dominated problems. The next three examples, in Sections 5.4, 5.5 and 5.6 are used to illustrate the excellent performance of the proposed elements in elastoplastic problems (including the case of perfect plasticity). The following three examples, Sections 5.7,5.8 and 5.9, along with Section 5.6, are used to present the excellent performance of the proposed three-dimensional brick elements in shell-like problems. Finally, Section 5.10 illustrates the performance of the proposed elements in the context of finite strain elasticity. 5.1. Constant stress/strain patch test A common objective in element design is that regardless of how skewed the element is, it should model a state of constant strain exactly. This characteristic is demonstrated with the mesh presented in Reference 25. The displacements at the vertices are specified to coincide with a homogeneous state of strain (see Reference 25 for specified values), while the interior nodes' displacements are left free to be determined. The material constants used are: Young's modulus, E = 1.OE + 6, and Poisson's ratio, v = 0.25. All three proposed elements reproduce the exact nodal displacements as specified in MacNeal and Harderz5 and thus, pass this test. The same mesh is also used with only a minimum number of constrained displacements (to prevent rigid body modes) and applied external loads that yield a constant state of stresslstrain. The proposed elements pass this version of the patch test, too. 5.2. Eigenvalue analysis of a single element To test the elements' performance at the nearly incompressible limit, and to show that they are free of spurious zero energy modes, an eigenvalue analysis of a single element mesh is performed. First, the eigenvalues of an undistorted element of side length L = 1 are presented in Table I.$$ Second, the eigenvalues of a rhombic-shaped element of side length L = 1 and distortion angle @ = 45" are given in Table 11. A perspective view of the elements' geometry is shown in Figure 1. The material properties used in both cases are: E = 1 and v = 0.4999. It is evident from the tables that all elements possess the correct non-zero eigenvalues and that all three are free of spurious locking behaviour at the nearly incompressible limit. The B8-9P element, as expected, exhibits a higher sensitivity to mesh distortion, as manifested by the higher eigenvalues (2-18). This result stands in sharp contrast to the nine-parameter elements of Andelfinger et al.' and Simo and Armero.2 These latter elements possess four large eigenvalues, three of which are undesired, indicating possible volumetric locking in some cases. Moreover, the ~ 15-parameter element of Andelfinger et al.' and the 12-parameter element of Simo et ~ 1 . ' possess four and three large eigenvalues, respectively. Locking behaviour of the Simo et a1.26element can be avoided, however, through a modification of the isoparametric shape function derivatives (see "The eigenvalues 19-24 are zero for both cases considered (i.e. the six rigid body modes) 2348 S. L. WEISSMAN Table I. Eigenvalues for an undistorted element of side length L = 1 Eigenvalue All elements 1 2500E +3 2 3.333E - 1 3 3.333E - 1 4 3.333E - 1 5 3.333E - 1 6. 3.333E - 1 Eigenvalue All elements 7 3.333E - 1 8 3.333E - 1 9 3.333E - 1 10 2.222E - 1 11 1.111E - 1 12 1.111E - 1 Eigenvalue All elements 13 l.111E - 1 14 5.555E - 2 15 5.5558 - 2 16 5-555E - 2 17 5-555E - 2 18 5555E - 2 Table 11. Eigenvalues for a rhombic - shaped element of side length L = 1 and 4 = 45" Eigenvalue B8-24/15P B8-9P 1 3.333E + 3 3.333E 3 2 6.372E - 1 7.690E - 1 3 5695E - 1 6.724E - 1 4 5WOE - 1 6.372E - 1 5 6 5'OOOE - 1 4303E - 1 5'000E - 1 5.000E - 1 Eigenvalue B8-24/15P B8-9P 7 1.962E - 1 3.934E - 1 . 8 1.897E - 1 3.545E - 1 1-603E - 1 3'150E - 1 9 10 1.3278 - 1 2'385E - 1 11 8.334E - 2 1'962E - 1 12 8.012E - 2 1.766E - 1 Eigenvalue B8-24/15P B8-9P 13 5.9658 - 2 8.334E - 2 14 15 2W0E - 2 6.647E - 2 1'069E - 2 4.755E - 2 16 9.849E - 3 2.508 - 2 17 4.277E - 3 2.387E - 2 18 3.736E - 3 1.607E - 2 + Figure 1. Perspective views of an undistorted and a rhombic-shaped element the eigenvalue analysis in References 1 and 26). A less expected result is that the BS-24P and B8-15P elements yield identical results even for the distorted case. Finally, it is noted that the eigenvalues obtained for the Pian and Tong3 element are identical to those obtained for the B8-24P and B8-15P elements for both cases. 5.3. Bending of a cantilever beam: Higher order patch test This is a standard benchmark problem used to assess element sensitivity to mesh distortion. In this problem, a cantilever beam of length L = 10, height H = 2, and width W = 1 is loaded by an end unit moment. To highlight the proposed elements' performance in the presence of internal 2349 HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONALBRICK ELEMENTS 10.0 I, 0.5 0.5 Figure 2. Higher-order patch test mesh, plan view c 5 1.0 g 0 9 0.8 0 a BN ._ 0.6 ‘;ii 0.4 E b 0.2 I 0.0 I 0 I 1 . I I I 1 2 3 4 5 A Figure 3. Higher order patch test, v = 0.49999; vertical displacement of the bottom edge, normalized with beam solution constraints, the lateral displacements are fixed (i.e. a state of plane strain is simulated) and the material constants are E = 1 and v = 0.4999 (i.e. nearly incompressible material). The beam is modelled by two elements, as shown in Figure 2. Figure 3 shows the lower edge vertical displacement, normalized with the exact beam theory solution of 56.2565, as A is increased from 0 to 5 (see Figure 2). To further emphasize that the deterioration in the elements’performance is due to geometrical distortion, the same problem is repeated for the case of v = 0. The lower edge vertical displacement (normalized with the exact beam solution of 75.0)is presented in Figure 4. It is evident from comparing Figures 3 and 4 that the deterioration in performance is indeed due to the geometrical distortion. In fact, the performance of the B8-9P element is slightly better for the nearly incompressiblecase. The B8-15PYB8-24P and the Pian and Tong3 elements yield identical results, and exhibit excellent insensitivity to mesh distortion. This identity is, however, restricted to linear elastic constitutive law. Finally, the B9-9P element reveals significant sensitivity to excessive mesh distortion and thus, as expected, it is shown to be inferior relative to the other proposed elements. 5.4. Conjined plastic flow This problem was suggested by Andelfinger et al.’ as a challenging problem for elements’ performance under confined plastic flow. In this problem, a ‘regular brick’ (termed in Andelfinger 2350 S. L. WEISSMAN I 0.04 0 . I I 1 2 I . 3 4 5 A Figure 4. Higher order patch test, v = 0; vertical displacement of the bottom edge, normalized with beam solution Figure 5. A regular mesh of 5 x 5 x 5 elementsis used to model (taking advantage of symmetry conditions)a regular brick of side length L = 100 and height h = 50, loaded by a uniform pressure load of 4 = 250 acting on its center area of 20 x 20 et al. as ‘regular block‘) of side length L = 100 and height h = 50 is fixed at its bottom, and a uniform pressure of q = 250 is applied to the center 20 x 20 area of its top surface. Using symmetry conditions, only one quadrant is discretized by a uniform 5 x 5 x 5 mesh, as shown in Figure 5. The material properties are: E = 210000-0, v = 0-3, c,,= 250-0, B = 1.0 (i.e. isotropic hardening), and R = 1/210. Figure 6 presents the load-deflection curves obtained for the proposed elements and for the displacement model. While the three proposed elements yield slightly different curves, they all have the same slope at the end (this difference diminishes under mesh refinement).This is in sharp contrast to load-deflection curve for the displacement model, which exhibits a higher slope. The failure of the displacement model to achieve the correct slope is in agreement with the findings of Nagtegaal et al.’ for the double-edge-notched tensile specimen, where this behaviour HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2351 L 0 CI 0 rz I............. DlSP Figure 6. Load-deflection curves for simulations of the regular brick problem using various elements I 48 1 I Figure 7. Cook's membrane problem (plan view) was observed for very fine meshes. It should also be noted that the nine-parameter element of Andelfinger et al.' also overestimates the load carrying capacity of the brick. However, the twenty-one- and thirty-mode elements of Andelfinger et al.' yield results similar to the B8-15P and B8-24P elements. 5.5. Cook's membrane problem: elastoplastic A tapered panel (of unit thickness) is clamped at one end and loaded by an in-plane bending load at the other end, as shown (plan view) in Figure 7. This benchmark problem demonstrates the proposed elements' superior coarse mesh accuracy in bending-dominated elastoplastic problems. The material properties are: E = 70, v = 3, o,, = 0.243, = 0.1,and ff = 0-15. A plane strain simulation is imposed by fixing the out-of-plane displacements. Figure 8 presents the convergence (of the centre node on the loaded edge) under uniform mesh refinement. For comparison, the results obtained for the plane strain implementation of the mean dilation (B-bar) element of 2352 S. L. WEISSMAN 2 4 8 16 32 NELlSide Figure 8. Cook‘s membrane problem+lastoplastic plane strain simulation.Convergencefor the proposed elements and the B-bar (plane strain) element Nagtegaal et aL5 are provided.@The superiority of the proposed elements is evident from the figures. It is interesting to note that even the B8-9P element exhibits superior performance over the mean dilation element. 5.6. Clamped square plate: elastoplastic This example, typically modelled by plate and shell elements, is introduced to demonstrate the excellent performance of the proposed elements even when the element’s geometrical aspect ratio is large. The plate’s characteristics are: side length L = 10 and thickness t = 0 1 (a plan view is shown in Figure 9); the material properties are E = 10920.0, v = 0.3 and by = 1000.0 (i.e. an elastic-perfectly plastic material); and the plate is loaded by a uniformly distributed pressure. First, the convergence under uniform mesh refinement (with only one element through the thickness) for the elastic case is considered. The transverse displacement of the mid-surface at the plate’s centre is presented in Figure 10 (normalized with the Kirchhoff solution of w = 12.6). The results obtained for the T1 plate element” are given for comparison. It is evident that the three-dimensional brick elements compare favourably with the plate element results. Second, the load-deflection curves for the mesh of the 16 x 16 (in-plane)elements with two and four elements through the thickness are presented in Figure 11(stress control simulation). Comparison is made with the upper Kirchhoff limit load,22given by Fuppcr = (45-05/& h2ay (49) where h is the plate thickness. The results for the mesh with only one element through the thickness are not reported since it does not have enough integration points through the thickness. As demonstrated by Figure 10, the fact that there are only two integration points through the thickness does not affect the accuracy in the (linear)elastic case. However, during the elastoplastic #Recently a number of elements have been proposed that outperform the mean dilation element. However, the mean dilation element remains the standard benchmark in the literature. Therefore, the mean dilation element is used here to provide a measuring stick for the performance of the proposed elements HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2353 Figure 9. Clamped square plate of side length L = 10, thickness h = 0.1. Due to symmetry only one quadrant is discretized (plan view) 0.950 4 256 64 16 1024 NEL Figure 10. Uniformlyloaded square clamped plate; convergenceof the centre transverse.displacementwith uniform mesh refinement B8-15P14 ----o---. om+ 0 1 5 . I . . I . a - - - . - . 1 I I 1 - I -1 10 15 20 25 30 35 40 45 50 55 60 65 W Figure 1 1. Load-displacement curves for the clamped square plate-dastoplastic 2354 S.L. WEISSMAN (a) (b) Figure 12. Circular plate meshes-plan view I3 .--B Q 1.00- 0.75 - U - 1 0.50 1 Y L - U 41 .- 0.25 // I____----- - Q E c p 1 o.m! 1 I I 100 Number of elements in the mesh 10 lo00 Figure 13. Clamped thin circular plate--convergence of the transverse displacement at R refinement =0 with uniform mesh phase, more integration points through the thickness are required. In particular, when only one element is used through the thickness, the integration points are relatively far from the upper and lower surfaces; thus, the model predicts elastic, and consequently stiffer, response even when the lower and upper surfaces have plastified. 5.7. Clamped thin circular plate In this example a thin clamped circular plate of radius R = 5 and thickness t = 0.1 is subject to a uniformly distributed unit load. This classical benchmark problem is used to demonstrate that plate and shell elements are free from shear locking when irregular shaped elements are used. Taking advantage of symmetry conditions, the plate is meshed using 3, 12,48 and 768 elements (only one element through the thickness). Two mesh topology options, illustrated in Figure 12, are examined. The linear elastic material properties are E = 10920.0 and v = 03. The convergence of the transverse displacement at the plate’s centre (R = 0) with uniform mesh refinement is shown in Figure 13 (normalized with the Kirchhoff solution of w = 9.78348).The HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2355 t’ I F = 1.0 Figure 14. Pinched cylinder with end diaphragms, perspective view results obtained for the T1 plate element2’ are given for comparison. The plate element outperforms the brick elements. This result is expected as plate elements are specifically designed for this class of problems. However, the proposed brick elements produce very good results. With only twelve elements, 92 and 78 per cent of the Kirchhoff solution is obtained when using mesh a and mesh b, respectively. Note that all the proposed elements yield virtually identical results. Furthermore, observe the significant improvement in course meshes obtained when using meshes of type a over meshes of type b. Finally, it should be reiterated that for this linear elastic problem, the element of Pian and Tong3 yields identical results to these reported for the proposed elements. 5.8. Pinched cylinder with end diaphragms A cylinder of length L = 600, radius R = 300, and thickness t = 3, with two end rigid diaphragms, is loaded by pinching vertical forces at the middle section. The material properties are E = 3.0 E 6 and v = 0.3. The objectives of this problem are: (1) to test the proposed elements’ performance in a curved geometry, and (2) to test the importance of the geometrical approximation. To this end, using symmetry conditions, only one octant of the cylinder is discretized by uniform meshes, as shown in Figure 14. First, the convergence under uniform mesh refinement is examined (only one element is used through the thickness). Figure 15 presents the displacement under the load (normalized with an analytical solution of 1432488 E - 5). The results obtained for the displacement model eight-node brick element and for the shell element of Simo et aL2* (the mixed approximation) are provided . ~ outperforms ~ the proposed brick for comparison. The shell element of Simo et ~ 1 clearly elements. However, roughly the same number of elements is required to achieve a low level of error (e.g. 1 per cent). It should also be noted that the rate of convergence with mesh refinement is about the same for the shell element and the proposed brick elements. Thus, the proposed elements provide a viable alternative to shell elements. This alternative is particularly attractive in problems involving non-smooth geometry, such as the one described below. Second, the geometry is fixed to that obtained by the 4 x 4 mesh, and the convergence under uniform mesh refinement of the displacement under the load is presented in Figure 16 + 2356 S. L. WEISSMAN 16 NEWSIDE 8 4 Figure 15. Pinched cylinder with end diaphragms. Convergence of the displacement under the load with uniform mesh refinement - -- - - - B8-24/15/9P - ' 0.0 BE-DISP CONVERGED I I I I 4 8 16 32 NEWSIDE Figure 16. Pinched cylinder with end diaphragms-geometry fixed to that obtained by the 4 x 4 mesh. Convergence of the displacement under the load with uniform mesh refinement (normalized with a converged solution of 8-95 E - 6). The converged displacement in this case is only 49 per cent of the analytical solution obtained for the smooth cylinder. This result re-affirms the known sensitivity of curved shell-like structures to geometrical perturbations. Moreover, it demonstrates the need for the proposed elements since such geometry is outside the scope of classical shell theory (see e.g. Reference 29). 5.9. Pinched hemispherical shell with 18" hole The pinched hemispherical shell with an 18" hole, proposed in Reference 25, is used to evaluate the performance of the proposed elements in doubly curved structures. The shell's dimensions are: radius R = 10 and thickness t = 0.04.The material properties are E = 6.875 E + 7 and v = 0-3. The shell is subjected to alternating inward and outward forces (F = 2.0) 90" apart. Using symmetry boundary conditions, only one quadrant of the shell is modelled by uniformed meshes, as shown in Figure 17. HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2357 t= Figure 17. Hemispherical shell with 18" hole, perspective view -g 1.0- -.-EB 0.8n 0.6 FN t - 0.4 - zb 0.2 4 16 32 NEL/SIDE Figure 18. Hemispherical shell with 18" hole. Convergence of the displacement under the load with uniform mesh refinement Figure 18 presents the convergence of the displacement under the loads (normalizedwith 0.093, a value quoted in Reference 28 as obtained by asymptotic expansiony for the proposed elements. The results reported in Reference 28 for the mixed shell approximation and those obtained for the displacement model eight-node brick elements are also provided for comparison. Once again the shell element exhibits superior performance. However, as in the former example, roughly the same number of elements is required to achieve a low level of error. 5.10. Finite strain elasticity: Cook's problem The extension of the formulation presented herein to finite strain problems will be addressed in a subsequent paper. However, in the context of finite strain elasticity, formulated in the reference w MacNeal and Harderz5propose the value of 0.094 for the displacement under the load 2358 S.L. WEISSMAN configuration, this extension is straightforward. This example is intended to illustrate the performance of the proposed elements in this context. The changes required to transform the small strain formulation presented herein into a finite strain one is to add the geometrical stiffness,and to interpret the stress as the second Piola-Kirchhoff stress and the strain as the Lagrangian strain. To illustrate the performance of the proposed elements, consider once again Cook's problem of Example 5.5, now with a compressible neo-Hook model characterized by the following strain energy function: K w =Z ( t ( -~ 1)2 - In(J)) + P (tr(C) - 3) Figure 19. Finite deformation analysis of Cook's problem; deformed mesh and outline of the undeformed mesh (dashed line) 4 2 4 8 16 32 NEL./Side Figure 20. Finite deformation analysis of Cook's problem; Convergence of the finite element solution HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2359 where C is the Cauchy-Green strain tensor, J’ is the determinant of C, and K and p are the bulk and shear moduli, respectively. As in Example 5.5, a plane strain condition is imposed by fixing the out-of-plane displacements. The material properties are: K = 400942.0 and p = 80.1938 (i.e. a nearly incompressiblematerial). The simulation is carried out using force control with steps of AF = 5.0 and F E [0.0,250*0]. The deformed shape for the 64-element mesh is shown in Figure 19, and the convergence of the top-corner vertical displacement with mesh refinement is given in Figure 20. 6. CONCLUDING REMARKS A family of three-dimensional eight-node brick elements that exhibit high accuracy in coarse meshes and are suitable for the modelling of elastoplastic (including elastic-perfectly plastic) problems has been presented. The elements are formulated based on a three-field functional of the Hu-Washizu type; consequently, the classical strain-driven algorithms of plasticity are directly applicable. This is in contrast with other highly accurate elements that are based on the Hellinger-Reissner functional and require the development of special algorithms (see, e.g. Reference 17). The latter are not suitable for the modelling of elastic-perfectly plastic problems. Also, the elements derived via the Hellinger-Reissner functional are based on the complementary energy, which is not always readily ascertainable. The elements presented pass the constant stress/strain patch test, and exhibit good distortion insensitivity. In addition, these elements exhibit good performance even when the element’s geometrical aspect ratio is very large, as shown by numerous problems commonly used to evaluate the performance of shell elements. Moreover, unlike the recently proposed elements of Andelfinger et al.’ and Simo and Armero’ the nine-parameter element proposed herein does not lock at the nearly incompressible limit. Finally, the B8-15P and B8-24P elements yield almost identical results. Therefore, since no advantage is to be gained by adding the extra nine modes, it is recommended that the B8-15P be chosen for practical computations. (The B8-24P element was presented here to demonstrate that adding the extra terms has only a marginal effect.) APPENDIX: RADIAL RETURN MAPPING AND THE ALGORITHMIC TANGENT The radial return mapping used to update the plastic strain and hardening parameters is summarized in this Appendix. Also provided is the algorithmic elastoplastic tangent. For a detailed discussion of integration algorithms for elastoplasticity and the algorithmic tangent, see Simo and Hughes.’’ A.I. Return mapping Step 1: Given L-:, a,, qn and evaluate the trial deviator stress and q trial - devCaln+1- 2pdev[En+1 - &:I and q::; Step 2. Evaluate the yield function = dev[a:?; - q,,] 2360 S. L. WEISSMAN Step 3. Check for plastic flow < 0)O then if(f:i set (0). + = (o):?: go to step 5 else go to step 4 end if Step 4. Compute Ay.+ = y.+ - y. and update plastic strain and hardening parameters where (58) Step 5. Compute the stress and algorithmic elastoplastic tangent Q,+~ = k trace(&,+ c;y1 where k = { + dev[o]fp\ (59) - 2pn::\ k 1 @ 1 + 2 , ~ ( 1 - +1@1) k l @ l+2~16~(1-31@1)-2p6,n:Y\@n::~ if:::f iff;?\ <0 >O (60) + $,u is the bulk modulus, and the constants b1 and d2 are given by =1 REFERENCES 1. U. Andelfinger, E. Ramm and D. Roehl, ‘2D- and 3D-enhanced assumed strain elements and their application in plasticity’, Proc. COMPLAST 111, Barcelona, Spain, 1992. 2. J. C. Simo and F. Armero, ‘Geometrically nonlinear enhanced strain mixed methods and the method of incompatible modes’, Int. j. numer. methods eng., 33, 1413-1449 (1992). 3. T. H. H. Pian and P. Tong, ‘Relations between incompatible displacement model and hybrid stress model’, Int. j. numer. methods eng., 22, 173-181 (1986). 4. E. F. Punch and S. N. Atluri, ‘Development and testing of stable, invariant, isoparametric curvilinear 2- and 3-D hybrid-stress elements’, Comp. Methods Appl. Mech. Eng., 47, 331-356 (1984). 5. J. C. Nagtegaal, D. M. Parks and J. C. Rice, ‘On numerically accurate finite element solutions in the fully plastic range’, Comp. Methods Appl. Mech. Eng., 4, 153-177 (1974). 6. T. J. R. Hughes, ‘Generalization of selective integration procedures to anisotropic and nonlinear media’, Int. j. numer. methods eng., 15, 1413-1418 (1980). 7. J. C. Simo, R. L. Taylor and K. S. Pister, ‘Variational and projection methods for volume constraint in finite deformation elastoplasticity’, Comp. Methods Appl. Mech. Eng., 51, 177-208 (1985). 8. D. S. Malkus and T. J. R. Hughes, ‘Mixed finite element methods-reduced and selective integration techniques: A unification of concepts’, Comp. Methods Appl. Mech. Eng., 15, 63-81 (1978). 9. J. C. Simo and M. S. Rifai, ‘A class of mixed assumed strain methods and the method of incompatible modes’, Int. j . numer. methods eng., 29, 1595-1638 (1990). HIGH-ACCURACY LOW-ORDER THREE-DIMENSIONAL BRICK ELEMENTS 2361 10. N. Bicanic and E. Hinton, ‘Spurious modes in two-dimensional isoparametric elements’, Int. j. numer. methods eng., 14, 1545-1557 (1979). 11. E. L. Wilson, R. L. F. 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. 12. I. Babushka, ‘Error bounds for finite element method‘, Numer. Math., 16, 322-333 (1971). 13. F. Brezzi, ‘On the existence, uniqueness and approximation of saddle-point problems arising from lagrange multipliers, Rev. Francaise d’dutomati. Informat. Recherche Operationnelle, &R2, 129-1 51 (1974). 14. B. F. de Veubeke, ‘Displacement and equilibrium models in the finite element method’, in 0.C. Zienkiewicz and G. S. Holister (eds.), Stress Analysis, John Wiley, London, 1965. 15. K. Washizu, Variational Methods in Elasticity and Plasticity, 3rd edn., Pergamon, Elmsford, NY, 1982. 16. M. L. Wilkins, ‘Calculation of Elastic-Plastic Flow’, in B. Alder, et al. (eds.), Methods of Computational Physics 3, Academic Press, New York, 1964. 17. J. C. Simo, J. G. Kennedy and R. L. Taylor, ‘Complementary mixed finite element formulations of elastoplasticity’, Comp. Methods Appl. Mech., Eng., 74, 177-206 (1989). 18. J. E. Marsden and T. J. R. Hughes, Mathematical Foundations of Elasticity, Englewood Cliffs, Prentice-Hall, N.J., 1983. 19. J. C. Simo and T. J. R. Hughes, Elastoplasticity and Viscoplasticity: Computational Aspects, Springer, Berlin, to appear. 20. S. L. Weissman, ‘A mixed formulation of nonlinear-elastic problems’, Comput. Struct., 44, 813-822 (1992). 21. S. L. Weissman and M. Jamjian, ‘Two-dimensional elastoplasticity: Approximation by mixed finite elements’, Int. j . numer. methods eng., 36, 3703-3727 (1993). 22. J. Lubliner, Plasticity Theory, Macmillan, New York, 1990. 23. T. J. R. Hughes, The Finite Element Method, Englewood Cliffs, Prentice-Hall, N.J., 1983. 24. T. H. H.Pian and K. Sumihara, ‘Rational approach for assumed stress finite elements’, Int. j . numer. methods eng., 20, 1685-1695 (1984). 25. R. H. MacNeal and R. L. Harder, ‘A proposed standard set of problems to test finite element accuracy’, Finite Elements Anal. Des., 1, 3-20 (1985). 26. J. C . Simo, F. Armero and R. L. Taylor, ‘Improved versions of assumed enhanced strain tri-linear elements for 3D-finite deformation problems’, Comp. Methods Appl. Mech. Eng., 110, 359-386 (1993). 27. T. J. R. Hughes and T. E. Tezduyar, ‘Finite elements based upon mindlin plate theory with particular reference to the four node bilinear isoparametric element’, J . Appl. Mech., 46, 587-596 (1981). 28. J. C. Simo, D. D. Fox and M. S. Rifai, ‘On a stress resultant geometrically exact shell model. Part 11: The linear theory; computational aspects’, Comp. Methods Appl. Mech. Eng., 73, 53-92 (1989). 29. P. M. Naghdi, ‘Theory of shells’, in C. Truesdell (ed.), Handbuck der Physik, Vol. Via/2, Mechanics of Solids 11, Springer, Berlin, 1983.

1/--страниц