INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 44, 615—639 (1999) A SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS FOR LINEAR AND NON-LINEAR MIXED-ENHANCED FINITE ELEMENTS FOR PLANE ELASTICITY PROBLEMS R. PILTNER* AND R. L. TAYLOR Department of Engineering Mechanics, ¼317.2 Nebraska Hall, ºniversity of Nebraska-¸incoln, NE 68588-0526, º.S.A. Department of Civil Engineering, SEMM, ºniversity of California at Berkeley, Berkeley, CA 94270, º.S.A. SUMMARY In a previous paper a modified Hu—Washizu variational formulation has been used to derive an accurate four node plane strain/stress finite element denoted QE2. For the mixed element QE2 two enhanced strain terms are used and the assumed stresses satisfy the equilibrium equations a priori for the linear elastic case. In this paper an alternative approach is discussed. The new formulation leads to the same accuracy for linear elastic problems as the QE2 element; however it turns out to be more efficient in numerical simulations, especially for large deformation problems. Using orthogonal stress and strain functions we derive B functions which avoid numerical inversion of matrices. The B -strain matrix is sparse and has the same structure as the strain matrix B obtained from a compatible displacement field. The implementation of the derived mixed element is basically the same as the one for a compatible displacement element. The only difference is that we have to compute a B -strain matrix instead of the standard B-matrix. Accordingly, existing subroutines for a compatible displacement element can be easily changed to obtain the mixed-enhanced finite element which yields a higher accuracy than the Q4 and QM6 elements. Copyright 1999 John Wiley & Sons, Ltd. KEY WORDS: mixed finite element; quadrilateral enhanced strain element 1. INTRODUCTION Several methods have been developed to improve the performance of the standard four-node compatible displacement element which yields poor results for problems with bending and, for plane strain problems, at the nearly incompressible limit. Among the improvements we like to cite the following methods: (i) Selected reduced integration and the so-called B-bar method of Hughes. (ii) The QM6-element of Taylor/Beresford/Wilson, who used four incompatible quadratic displacement functions to calculate an approximated displacement gradient. (iii) The method of ‘enhanced strains’ introduced by Simo and Rifai. Simo and Rifai also demonstrated that the QM6 element can be viewed as an enhanced strain element with four enhanced strain terms. * Correspondence to: R. Piltner, Department of Engineering Mechanics, W317.2 Nebraska Hall, University of NebraskaLincoln, Lincoln, NE 68588-0526, U.S.A. E-mail: [email protected] CCC 0029—5981/99/050615—25$17.50 Copyright 1999 John Wiley & Sons, Ltd. Received 8 October 1997 Revised 7 May 1998 616 R. PILTNER AND R. L. TAYLOR (iv) The hybrid stress elements of Pian and Sumihara and Yuan/Huang/Pian. (v) The mixed finite element QE2 of Piltner and Taylor using a modified Hu—Washizu variational formulation with bilinear displacement interpolations, seven strain and stress terms in Cartesian co-ordinates, and two enhanced strain modes. The enhanced concept became quite popular in recent years and has been used for both linear and non-linear problems by several researchers, e.g. Simo/Armero, Simo/Armero/Taylor, Andelfinger/Ramm, Crisfield/Moita/Jelenic/Lyons, Freischläger/Schweizerhof, Glaser/ Armero, Korelc/Wriggers, Nagtegaal/Fox, Roehl/Ramm, de Sausa Neto/Peric/ Huang/Owen, Wriggers/Reese. In Reference 1, the stress functions for the QE2 element were chosen such that the homogeneous equilibrium equations are satisfied a priori. Therefore the linear version of the QE2 element can also be considered as a Trefftz-type element (e.g. References 21—28). For the QE2 element a matrix H has to be computed and inverted. The matrix H has the following block diagonal structure: H 0 (1) 0 H where H is a 3;3 diagonal matrix. So with the assumed stresses of Reference 1 we are left at the element level to invert the 4;4 symmetric sub-matrix H . In this paper we consider a method to avoid numerical inversion (except the one needed for the static condensation of the unknowns associated to the enhanced strain functions). In addition to the improvement in efficiency of the mixed four-node finite element we address the subject of non-physical instabilities of non-linear versions of enhanced finite elements which have been observed by Wriggers and Reese and discussed by Crisfield et al. and De Sousa Neto et al. A mixed finite element with four enhanced strain terms has been tested in a series of numerical examples. The element is denoted B -QE4. In several numerical tests for linear problems it is shown that the elements B -QE4 and QE2 give identical results. The possibility of coding the B -QE4 element like a displacement element and its excellent performance make it attractive for non-linear applications. H" 2. FINITE ELEMENT FORMULATION FOR THE LINEAR ELASTIC CASE As in Reference 1 we use the following modified Hu—Washizu variational formulation: 4 2 e2Ee d»!4 u2f d»!1 u2T dS!4 r2(e!Du!eG) d» %(u, eG, e, r)" 1 (2) where eG is the enhanced strain field satisfying an orthogonality condition with respect to properly chosen reference stresses. Without eG we have the classical Hu—Washizu formulation. Details about the choice of enhanced strain terms are given in Reference 1. The motivation for using a mixed variational formulation is the following: The good performance of the element QE2 showed that there can be an advantage of assuming stresses in Cartesian co-ordinates. Using both assumed stresses and assumed strains we will be able to derive a formulation which allows a fast implementation for both linear and non-linear problems. It will be shown that this implementation of the mixed element will be very similar to an implementation of a pure displacement element. Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 617 Carrying out the variation in (2) we obtain the following equations: 4 dr2[e!Du!eG] d»"0 (3) 4 de2[Ee!r] d»"0 (4) ! 4 d(Du)2r d»"4 du2f d»#1 du2T dS 4 (deG)2r d»"0 (5) (6) The displacement, strain and stress fields are chosen in the following form: u"Nq, e"Sa, r"Sb, eG"BGk (7) where N"N(m, g) is the matrix of compatible shape functions and q contains the nodal displacements of the finite element. The vectors a, b, k are strain, stress and enhanced strain parameters, respectively. The strain field obtained from the compatible displacement field u can be written as Du"Bq (8) where D is a linear differential operator matrix. The discretization of equations (3)—(6) with arbitrary variations db, da, dq and dk leads to the following system of equations at the element level: 0 !H L LG b 0 a 0 L2 H 2 0 0 0 q LG2 0 0 0 k !H 0 0 " f 0 (9) where 4 S2S d», H" 4 S2ES d», H " 2 4 S2B d», L" 4 S2BG d» LG" 4 N2f d»#1 N2T dS f " (10) From the first two equations of (9) we obtain the strain and stress parameters as a"H\Lq#H\LGk (11) b"H\H a"H\H H\Lq#H\H H\LGk 2 2 2 (12) Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 618 R. PILTNER AND R. L. TAYLOR Substitution of (11) and (12) into (9) gives us the following system of equations: L2H\H H\L L2H\H H\LG 2 2 LG2H\H H\L LG2H\H H\LG 2 2 or K C2 C Q q f " k 0 q f " k 0 (13) (14) where K"L2H\H H\L, C"LG2H\H H\L, Q"LG2H\H H\LG 2 2 2 (15) Using a static condensation process with respect to the enhanced strain parameters k we finally get the element stiffness matrix in the form k"K!C2Q\C (16) In order to avoid calculating the matrices L, H, H , L, LG as well as performing the numerical 2 inversion of H and matrix multiplications in (15), we seek a more efficient way to compute the matrices K, C and Q. The key for an efficient implementation of the mixed finite element is to exploit the structures of the matrices H, L, LG for a proper choice of stress and strain functions collected in the matrix S. In order to achieve high accuracy with our finite element we choose the strain and stress functions in Cartesian co-ordinates. This choice is motivated by the good numerical results of the element QE2 and the Trefftz-type finite elements from References 21—28 for which stresses in Cartesian co-ordinates have been chosen such that the equilibrium equations are satisfied a priori. The disadvantage of the set of assumed stresses for element QE2 (see Equation (77) in Reference 1) is that the matrix of assumed stresses is full and did not lead to a diagonal matrix H. In Reference 1 a block diagonal matrix with one 3;3 diagonal sub-matrix and a 4;4 fully populated sub-matrix was constructed. For the present element we do not require that the assumed stresses satisfy equilibrium at the element level a priori for arbitrary stress parameters b, as was the case in Reference 1. Instead we require only that the stresses can satisfy the equilibrium equations for proper relations between the stress parameters b collected in the vector b. This implies that we now need more stress terms H than in Reference 1, where seven stress terms have been used. For the mixed two-dimensional element we can assume the stress field in the form p S S S 0 0 0 0 0 0 VV r" p " 0 0 0 S S S 0 0 0 [b]"Sb WW q 0 0 0 0 0 0 S S S VW (17) where the S are linearly independent functions in Cartesian co-ordinates. For the strain field we H use the same matrix S as for the stress trial field. Unlike the assumed stresses for element QE2 (see Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 619 equation (77) in Reference 1) the assumed stress components of (17) are not coupled. The resulting matrix H has the following block diagonal structure: Hª 0 0 H" 0 Hª 0 0 0 Hª (18) If our linearly independent trial functions S ( j"1, 2, 3) do not satisfy the orthogonality H condition 4 SGSH d» O0 "0 for iOj for i"j (19) we will get a fully populated 3;3 matrix Hª . Since the matrix H has to be inverted it is better to orthogonalize a set of linearly independent functions so that Hª and therefore H become diagonal. Possible choices for our initial (non-orthogonal) set of linearly independent trial functions S H are S "1, S "xN , S "yN (20) and S "1, S "m, S "g (21) where xN , yN are Cartesian co-ordinates with origin at the centre of the quadrilateral finite element. The co-ordinates m, g are obtained as a linear combination of the Cartesian co-ordinates xN , yN from the relationship m b !a xN 1 " (22) g J !b a yN The relationship between the global Cartesian co-ordinates x, y and the co-ordinates xN , yN , shifted to the center of the element, is given as xN "x!a , yN "y!b (23) where a " (x #x #x #x ), b " (y #y #y #y ) The coefficients of the linear mapping (22) are given as (24) a " (!x #x #x !x ), b " (!y #y #y !y ) a " (!x !x #x #x ), b " (!y !y #y #y ) (25) J "a b !a b Finally, we can express the element coordinates xN , yN and m, g through the natural element co-ordinates m, g via xN "a m#a g#a mg yN "b m#b g#b mg Copyright 1999 John Wiley & Sons, Ltd. (26) Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 620 R. PILTNER AND R. L. TAYLOR and J m"m# mg J J g "g# mg J (27) where a " (x !x #x !x ), b " (y !y #y !y ), (28) and J "a b !b a , J "a b !a b (29) Since the expressions (27) are simpler than (26) we choose [SK "1, SK "m, SK "g ] as initial set of linearly independent trial functions to construct an orthogonal set of linearly independent functions [S , S , S ]. This can be achieved by a Gram—Schmidt-orthogonalization process through the formula SK S d» H I H\ 4 S "SK ! S H H I I S S d» I I 4 (30) The resulting set of orthogonal functions is S "1 J J J 1 J S "m! "m# mg! " b xN !a yN ! (31) 3J J 3J J 3 J J J 1 J S "g ! !f S "g# mg! !f S " !b xN #a yN ! !f S J J 3J 3 3J where 2J J f" (32) 3J!J#3J It should be noted that the functions S and S used for the stress/strain approximations can be expressed as polynomials in xN , yN , m, g or m, g-co-ordinates, whereas in a pure displacement formulation the stresses and strains are rational functions in natural element co-ordinates for a general quadrilateral element. Because the functions S can be expressed as a complete set of H linear functions in Cartesian co-ordinates, the special case of the seven parameter stress field of the QE2 element satisfying the equilibrium equations a priori is included in the stress field (17). With the choice of functions [S , S , S ] from equation (31) the entries of the matrix S for the stress and strain fields are defined. Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 621 SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS The compatible displacement field u is interpolated with standard shape functions according u 1 u u " (1#m m)(1#g g) G " N G "Nq (33) G G G v v v 4 G G G G An admissible enhanced strain field with four terms can be obtained from the incompatible displacement field u" uG" N 0 N 0 0 N 0 N j j j j (34) where N"1!m, N"1!g (35) Because the constant strain patch test cannot be satisfied when using the exact gradient of the incompatible shape functions * y (m, g) !y (m, g) *x 1 E K N" N" H H * J(m, g) !x (m, g) x (m, g) E K *y * *m N"J\(m, g) H * *g * *m N H * *g (36) the approximated gradient * *x 1 y (0, 0) !y (0, 0) E K N" N" H H J(m, g) !x (0, 0) * x (0, 0) E K *y * *m J N" J\ H J(m, g) * *g * *m N H * *g (37) as proposed by Taylor/Beresford/Wilson in Reference 4 is used to get the following admissible enhanced strain field: eG !b m 0 b g 0 VV 2 eG" eG "BGk" 0 a m 0 !a g WW J(m, g) cG a m !b m !a g b g VW j j j j (38) where the determinant J of the Jacobian matrix can be expressed as J(m, g)"J #J m#J g (39) With the assumed displacement, strain and stress fields, all integrals in the element formulation can be calculated exactly with 2;2 Gauss-integration since the expression for the matrices H, H , 2 L, and LG contain only polynomials in m and g. This is in contrast to the situation in the Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 622 R. PILTNER AND R. L. TAYLOR displacement element formulation where we cannot get exact integral results for an arbitrary quadrilateral element shape. This is due to the factor 1/J(m, g) in all integrals for the stiffness matrix. The diagonal matrix H is defined through the entries of the 3;3 matrix Hª . Through analytical integration we get the diagonal coefficients of the matrix Hª as hK "4J Kh " J (3j!j#3) 4(3#2j!j#2j#2j j!j) J hK " 3(3j!j#3) where j "J /J and j "J /J . The strain matrices B and BG have the structure N 0 N 0 N 0 N 0 V V V V N 0 N 0 N 0 N B" 0 W W W W N N N N N N N N W V W V W V W V and * * N 0 N 0 *x *x BG" 0 * N *y 0 * N *y * N *y * N *x * N *y * N *x (40) (41) (42) where */*x and */*y indicate that the derivatives of the incompatible shape functions are an approximation given by equation (37). It is important to notice that due to the sparse structure of S, B, and BG and the fact that these matrices have repeated entries for their coefficients the matrices L and LG have only a few non-vanishing different coefficients. The non-vanishing coefficients are denoted by f GH , f GH , f GI , V W V f GI , where (i"1, 2, 3), ( j"1, 2, 3, 4), (k"1, 2) and given by W * f GH" S (m, g) N (m, g) d» V G *x H 4 4 * S (m, g) N (m, g) d» G *y H 4 * S (m, g) N (m, g) d» G *x I f GH" W f GI " V (43) * 4 SG (m, g) *y NI (m, g) d» f GI " W Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 623 There is no particular advantage to give analytical expressions for the results of the integrals in equations (43). The exact values for f GH , f GH , f GI , and f GI can be obtained with 2;2 Gaussian V W V W integration. Having calculated the exact coefficients of the diagonal matrix H and the sparse matrices L and LG we can calculate the strain parameters a from relationship (11). Now we can write the assumed strain field in the form: e"B q#B Gk (44) NM 0 NM 0 NM 0 NM 0 V V V V B " 0 NM 0 NM 0 NM 0 NM W W W W NM NM NM NM NM NM NM NM V W V W V W V W (45) where and NM V B G" 0 0 NM W NM V NM W NM V 0 NM W 0 NM W NM V (46) The matrices B , B G have the same sparse structure as the matrices B, BG. Taking into account that S "1, the NM functions can be written as f H f H f H NM H (m, g)" V # V S (m, g)# V S (m, g) V hK hK hK f H f H f H NM H (m, g)" W # W S (m, g)# W S (m, g) W hK hK hK f I f I f I NM I (m, g)" V# V S (m, g)# V S (m, g) V hK hK hK f I f I f I NM I (m, g)" W# W S (m, g)# W S (m, g) W hK hK hK (47) where ( j"1, 2, 3, 4) and (k"1, 2). Note that the matrices B and B G contain polynomials in m, g whereas the matrices B and BG contain rational functions of the natural element co-ordinates m and g. Constant stress and strain modes can be represented exactly with the proposed formulation involving the matrices B and B G. Instead of computing the matrices K, C and Q from equations (15), which involves time consuming matrix multiplications of H, L, LG, H we can compute them faster by using the sparse 2 strain matrices B and B G: 4 B 2EB d», K" Copyright 1999 John Wiley & Sons, Ltd. 4 B G2EB d», C" 4 B G2EB G d» Q" (48) Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 624 R. PILTNER AND R. L. TAYLOR The external load vector involves compatible shape functions from equation (33) and is given by 4 N2f d»#1 N2T dS f " (49) 3. FORMULATION FOR LARGE DEFORMATION OF HYPERELASTIC MATERIALS In this section we will first use tensor notation and later switch to matrix and vector notation. For an overview on non-linear elasticity we refer to Beatty’s review articles. For the large deformation case we will use the first (un-symmetric) Piola—Kirchhoff tensor P which is related to the second (symmetric) Piola—Kirchhoff stress tensor S , the Kirchhoff stresses s and the Cauchy stresses r through the relation P"FS "sF\2"JrF\2 (50) where F is the deformation gradient *x F" "I#Grad u *X (51) In this notation, x is the position vector in the current configuration whereas X denotes the position vector in the reference configuration. Using u as the displacement vector, the position vectors are related through the equation x"X#u (52) The displacement gradient with respect to the reference configuration is given for our twodimensional example by *u *u *X *½ Grad u" (53) *v *v *X *½ The relations for the symmetric Kirchhoff stress tensor s and the Cauchy stress tensor r are given by s"PF2"FS F2 1 r" s J (54) Here we want to consider the case of a hyperelastic material which is characterized through a frame invariant stored energy function ¼(X, F). For objectivity the strain energy function must satisfy the relation W(X, QF)"¼(X, F) (55) for all proper orthogonal Q. For a hyperelastic material the un-symmetric Piola—Kirchhoff stress tensor is given by *¼ P" *F Copyright 1999 John Wiley & Sons, Ltd. (56) Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 625 In the finite element implementation we will make use of the incremental relation *P"A : *F (57) *P "A *F G( G(I* I* (58) or in index notation The fourth order tensor A, which is called the first elasticity tensor, is given through *¼ *¼ or A " A" G(I* *F *F *F*F G( I* (59) Using an assumed displacement gradient G and an enhanced displacement gradient Grad u in the reference configuration the modified Hu—Washizu formulation for large deformations can be given as 4 Grad du : P d»#d% . d%" (60) 4 dP : [G!Grad u!Grad u ] d» ! 4 dG : *F > *¼ # F I !P d» G 4 d[Grad u ] : P d» # where dF"Grad du and *u *u *X *½ uJ uJ G " 6 7 , Grad u " vJ vJ *v *v 6 7 *X *½ P" P P P P (61) (62) have been used. For a discretization of (60) we switch to matrix and vector notation. Using the vectors P2"[P , P , P , P ] u2"[u , v , u , v ] 6 7 7 6 g2"[uJ , vJ , uJ , vJ ] 6 7 7 6 u2 "[u , v , u , v] 6 7 7 6 Copyright 1999 John Wiley & Sons, Ltd. (63) Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 626 R. PILTNER AND R. L. TAYLOR The variational formulation can be expressed as 4 du2P d»#d% . d%" (64) 4 dP2 [g2! u! u ] d» ! 4 dg [P (g)!P] d» # 4 d u2P d» # where P denotes the vector of Piola—Kirchhoff stresses obtained from *¼/*F by using the assumed displacement gradient G . An increment for the Piola—Kirchhoff stress depending on an increment *g can be calculated from *P "A *g 2 (65) An expression for the tangent moduli A for one of the chosen constitutive models is given in 2 Section 4. For the increments of stresses and displacement gradients we assume the functions *P"S*b, *g "S*a, (*u)"B*q, (*u )"B*k (66) S S S 0 0 0 0 0 0 0 0 0 0 0 0 S S S 0 0 0 0 0 0 S" 0 0 0 0 0 0 S S S 0 0 0 0 0 0 0 0 0 0 0 0 S S S (67) where N B" 6 0 N 7 0 0 N 7 0 N 6 N 6 0 N 7 0 0 N 7 0 N 6 N 6 0 N 7 0 0 N 7 0 N 6 N 6 0 N 7 0 0 N 7 0 (68) N 6 B is the matrix for the enhanced displacement gradient. Unfortunately, for the large deformation case we cannot use a sparse matrix for the enhanced displacement gradient. Using a sparse matrix B yields an element which is not frame invariant. The following B has been chosen for Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 627 the proposed enhanced element: bm !b b m !b b g bg am !a a m !a a g ag 1 B" J(m, g) !a b m a b m a b g !a b g !a b m a b m a b g !a b g (69) The above matrix can be obtained by using the following form for the enhanced displacement gradient H: J F"H" J\H J\2 J (70) where the Jacobian at the element origin is given by a b J " a b and H is chosen as (71) j m j m (72) j g j g Recently instabilities have been detected in a non-linear version of the QM6 element. The QM6 element proposed by Taylor/Beresford/Wilson in 1976 can be viewed as an enhanced strain element with four enhanced strain terms. A non-linear version was proposed by Simo and Armero. The non-physical instabilities in enhanced strain finite element simulations can occur in a hyperelastic material under uniform compression, as initially observed by Wriggers and Reese. The problem has been investigated recently by several researchers, e.g. Crisfield et al., de Souza Neto et al., Glaser and Armero. In order to overcome the element instability problems, Glaser and Armero proposed to exchange the original basis for the enhanced deformation gradient H " H " j m j g j m j g (73) by either expression (72) or by j m j m#j g (74) j m#j g j g Glaser and Armero considered the following two transformations for constructing the enhanced part of the deformation gradient: H " Copyright 1999 John Wiley & Sons, Ltd. J (a) H" J H J\ J (75) J (b) H" J\2 H J\ J (76) Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 628 R. PILTNER AND R. L. TAYLOR The enhanced deformation gradient in Reference 13 is chosen in the form F"F H (77) as proposed by Simo/Armero and Taylor. Korelc and Wriggers proposed an enhanced deformation gradient of the form J F"H" J\2 H J\ J (78) where H is the same as in equation (72). However, Glaser and Armero pointed out that the use of (78) yields a non-objective element. For the proposed element H is chosen as in equation (70). Using (70) an objective element is obtained (see example 5.7). After elimination of the strain parameters the assumed displacement gradient for an increment takes the form *g "B *q#B *k (79) where B " NM 6 0 NM 7 0 NM 6 0 0 NM 7 0 NM 7 0 NM 6 0 NM 7 0 NM 6 NM 6 0 NM 7 0 0 NM 7 0 NM 6 NM 6 0 NM 7 0 0 NM 7 0 (80) NM 6 and bM bM B " bM bM bM bM bM bM bM bM bM bM bM bM bM bM (81) The functions bM (m, g) are defined as GH ¸ ¸ ¸ bM (m, g)" GH # GH S (m, g)# GH S (m, g) GH hK hK hK (82) and 4 SI (m, g)bGH (m, g) d» ¸I " GH (83) The functions b are the coefficients of the matrix B defined with equation (69). GH Because the assumed displacement gradient can be expressed in the form (79) involving the sparse matrix B , the element tangent stiffness matrix can be computed in a very economic manner. The element tangent stiffness matrix for the large deformation case is given as k "K !C2 Q\C 2 2 2 2 2 Copyright 1999 John Wiley & Sons, Ltd. (84) Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 629 where the matrices K , C , and Q can be computed in the following form: 2 2 2 4 B 2A2 B d» K " 2 4 B A2 B d» C " 2 (85) 4 B 2A2B d» Q " 2 The residual vectors for step n are obtained by computing 4 B 2P d», f " 4 B 2PL d» f " (86) 4. CHOSEN CONSTITUTIVE MODELS For the examples in the present paper we consider two constitutive models. 4.1. Model 1 The stored energy function for the first constitutive model is K (I , J)" k(I !3)!k ln J# j(ln J) ¼"¼(I , I , I )"¼ where I , I , I are the invariants of the right Cauchy—Green deformation tensor C"F2F (87) (88) and J"(I "(det C"det F The second Piola—Kirchhoff stress tensor is obtained from S "2 *¼ *I *¼ *I *¼ *I *¼ #2 #2 "2 *C *I *C *I *C *I *C (89) (90) where *I *I "I, "(det C)C\"JC\"JF\F\2 *C *C (91) For the above strain energy function ¼ we have *¼ 1 *¼ *¼ *¼ *J *¼ 1 " k, "0, " " *I 2 *I *I *J *I *J 2J Copyright 1999 John Wiley & Sons, Ltd. (92) Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 630 R. PILTNER AND R. L. TAYLOR and 1 *¼ k " #j ln J J J *J (93) so that we get the constitutive relationship for S in the form S "kI#[!k#jln J]F\F\2 (94) For the other stresses we obtain P"k[F!F\2]#jln JF\2, s"k[b!I]#jln JI, k j r" [b!I]# ln JI J J (95) where b"FF2 (96) The first elasticity tensor A can be obtained conveniently by using components of the unsymmetric Piola—Kirchhoff stress tensor P, given from (56) by 1 1 P "kF # [!k#jln J]F , P "kF # [!k#jln J]F J J 1 1 P "kF ! [!k#jln J]F , P "kF ! [!k#jln J]F J J (97) Using relationship (59) we get A 1 "k# [k#j!jln J]F F , J 1 A "k# [k#j!jln J]F F J 1 1 A " [!k#jln J]# [k#j!jln J]F F J J 1 1 "! [k#j!jln J]F F A "! [k#j!jln J]F F , A J J (98) 1 1 A "! [k#j!jln J]F F , A "! [k#j!jln J]F F J J 1 A "k# [k#j!jln J]F F , J 1 A "k# [k#j!jln J]F F J 1 1 A "! [!k#jln J]# [k#j!jln J]F F J J In matrix notation the incremental constitutive equation for a point under consideration can be written as *P"A *F 2 Copyright 1999 John Wiley & Sons, Ltd. (99) Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 631 where A A A A A A A A A " 2 A A A A A A A A *P2"[P , P , P , P ] (100) *F2" (*u2)"[*u , *v , *u , *v ] 6 7 7 6 It should be mentioned that although the symmetric matrix A is singular in the undeformed 2 configuration, the element stiffness matrix has the correct number of zero eigenvalues (i.e. three zero eigenvalues associated to the rigid body modes). 4.2. Model 2 The second constitutive model used in the examples is taken from a paper by Knowles and Sternberg. The strain energy function for this model is ¼(I , J)" k[I J\#2J!4] (101) Reference 31 gives a detailed discussion of this material model which is a special case of a Blatz—Ko material. 5. NUMERICAL EXAMPLES Several linear and non-linear problems have been selected to test the performance of the proposed four-node element denoted as B -QE4. For all problems 2;2 Gaussian integration is used. The element passes the patch test. 5.1. Beam bending A beam modeled with five elements is subjected to two load cases (Figure 1). Plane stress conditions are assumed in the model. The results of different elements for the maximum displacement at point A and the normal stress p at point B are given in Table I. The elements VV used in the examples are: (i) the bilinear isoparametric displacement element Q4 [32, 33, 34]. (ii) the enhanced strain element QM6 of Taylor/Beresford/Wilson, (iii) the hybrid stress element P—S of Pian/Sumihara, (iv) the enhanced mixed element QE2 by Piltner/Taylor, (v) the proposed element B -QE4. The elements B -QE4 and QE2 give identical results. 5.2. Mesh distortion test for beam bending In this test a beam under bending is analysed with only two plane stress elements (Figure 2). The degree of distortion of the element is measured with the distortion parameter *. The material parameters are E"1500 and l"0)25. From Table II we can see that the elements QE2 and B -QE4 show the least sensitivity to mesh distortion even for every severe distortions. Elements QE2 and B -QE4 provide identical results. Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 632 R. PILTNER AND R. L. TAYLOR Figure 1. Finite element mesh for cantilever beam problem Table I. Comparison of plane stress solutions obtained with four node elements for cantilever beam problems Case 1 Element v Q4 QM6 P-S QE2 B -QE4 Exact 45)65 96)07 96)18 96)5 96)5 100 Case 2 p VV !1604 !2497 !3001 !3004 !3004 !3000 v 50)96 97)98 98)05 98)26 98)26 102)6 p VV !2146 !3235 !3899 !3906 !3906 !4050 Figure 2. Cantilever beam for the mesh distortion test 5.3. Cook’s membrane problem The material parameters for the plane stress structure shown in Figure 3 are E"1 and l". Displacement and stress results are listed in Table III. Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 633 Table II. Displacenment v of cantilever beam (Figure 2) for different values of the mesh distortion parameter * * Q4 QM6 0 0)5 1 2 3 4 4.9 28)0 21)0 14)1 9)7 8)3 7)2 6)2 100)0 80)9 62)7 54)4 53)6 51)2 46)8 P-S QE2 Displacement v 100)0 100)0 81)0 81)2 62)9 63)4 55)0 56)5 54)7 57)5 53)1 57)9 49)8 56)9 B -QE4 Exact 100)0 81)2 63)4 56)5 57)5 57)9 56)9 100 100 100 100 100 100 100 Figure 3. Cook’s membrane problem: plane stress structure with unit load uniformly distributed along right edge (E"1, l") 5.4. Short beam under bending A short beam, supported at the left end and subjected to a couple at the right edge, is modeled with a mesh as shown in Figure 4. The purpose of this example is to test the influence of the mesh pattern on the stress distribution in the finite element solution. In the QM6 element the stress has the tendency to line up with the element edges of the element in the center of the mesh (Figure 5(a)). On the other hand, elements QE2 and B -QE4 preserve symmetries in the solution domain although a very coarse irregular mesh is used (Figure 5(b)). 5.5. Single element test for a Neo-Hookean material model Wriggers and Reese noticed that the enhanced element Q1/E4 can experience problems in the analysis of Neo-Hookean hyperelastic materials under uniform compression. They suggested the single element test shown in Figure 6. The material model (87) is used with the following material parameters: j"40 000, k"80)2. At 30)4 per cent compression the element tangent Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 634 R. PILTNER AND R. L. TAYLOR Table III. Results for the problem shown in Figure 3 Element N"2 N"4 N"16 Q4 QM6 P-S QE2 B -QE4 Displacement 11)85 21)05 21)13 21)35 21)35 v at C 18)30 23)02 23)02 23)04 23)04 23)43 23)88 23)88 23)88 23)88 Q4 QM6 P-S QE2 B -QE4 Maximum stress at A 0)1078 0)1814 0)1773 0)2225 0)1854 0)2241 0)1956 0)2261 0)1956 0)2261 0)2353 0)2364 0)2364 0)2364 0)2364 Figure 4. Finite element mesh for a short beam under bending matrix of the Q1/E4 element has two negative eigenvalues. However, the element tangent matrix should have only one negative eigenvalue for this deformation. At 39 per cent compression the element Q1/E4 yields a non-physical instability. The system tangent matrix has one negative eigenvalue at this deformation level. On the other hand the element B -QE4 keeps only one negative eigenvalue for the element tangent matrix, up to almost 100 per cent compression, and the system matrix has no negative eigenvalues. 5.6. Multi-element compression test for a Neo-Hookean material model One-half of a block with unit sides is discretized with a 10;20 element mesh. The material parameters for the Neo-Hookean model are chosen as in example 5.5. The element Q1/E4 yields a non-physical instability at 30)6 per cent. However, the first physical instability should appear at Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 635 Figure 5. (a) Stress p obtained with element QM6; (b) Stress p obtained with element QE2 and B -QE4 VV VV Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 636 R. PILTNER AND R. L. TAYLOR Figure 6. Finite element under compression Figure 7. (a) undeformed beam; (b) deformed beam and FE mesh 50)12 per cent compression. The element B -QE4 is able to detect this physical instability and the following instabilities: At 50)12 per cent compression the system matrix based on the B -QE4 elements has one negative eigenvalue. At 50)26 per cent compression the system matrix has two negative eigenvalues, and at 58)2 per cent compression the discretization with B -QE4 elements shows three negative eigenvalues. 5.7. Objectivity test for large deformation The beam shown in Figure 7(a) is fixed at the right end and subjected to a prescribed displacement as illustrated in Figure 7(b). After calculating the displacements and stresses the Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 637 Figure 8. Finite element discretization for example 5)8 Figure 9. Material 2. Cauchy stress vs. axial stretch structure is rotated and displacements and stresses are calculated for each rotation angle (h"10°, 20°, . . . , 180°). For the B -QE4 element the results are independent of the rotation angle. It should be noticed that the results are not frame independent if we use enhanced strains based on equation (76) instead of (70). 5.8. Material 2 subjected to a homogenous plane deformation One-half of a block of dimensions 1;1 is discretized with non-rectangular elements (Figure 8). The block is subjected to plane strain uni-axial stress parallel to the x -axis. Material model 2 is used with k"100. The analytical solution for p in this problem is given as a function of the stretch j (Figure 9): p "100(1!j\) (102) For every element in the mesh stresses can be calculated exactly. Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) 638 R. PILTNER AND R. L. TAYLOR In Figure 9, the solid line represents the analytical solution and the circles stand for stress values obtained numerically for a chosen finite element in the mesh. 6. CONCLUDING REMARKS For linear problems, the elements B -QE4 and QE2 give identical results and show an overall very good behaviour. The advantage of element B -QE4 is that it can be coded very similar to a displacement element and therefore runs faster than element QE2. The key for a fast linear and non-linear element implementation is to use orthogonal stress and strain functions in the mixed variational formulation. Due to the use of orthogonal functions matrix inversions for strain and stress parameters at the element level can be avoided in a finite element program. In the current implementation standard shape function derivatives N , N are exchanged by NM and NM V W V W functions derived in the paper. For non-linear problems the displacement like implementation of the mixed element B -QE4 has the obvious advantage that it is time saving compared to an implementation of the mixed variational formulation with matrix inversions for strain and stress parameters. REFERENCES 1. R. Piltner and R. L. Taylor, ‘A quadrilateral finite element with two enhanced strain modes’, Int. J. Numer. Meth. Engng., 38, 1783—1808 (1995). 2. T. J. R. Hughes, ‘Equivalence of finite elements for nearly incompressible elasticity’, J. Appl. Mech., 44, 181—193 (1977). 3. T. J. R. Hughes, ‘Generalizing of selective integration procedures to anisotropic and nonlinear media’, Int. J. Numer. Meth. Engng., 15, 1413—1418 (1980). 4. R. L. Taylor, P. J. Beresford and E. L. Wilson, ‘A non-conforming element for stress analysis’, Int. J. Numer. Meth. Engng., 10, 1211—1219 (1976). 5. J. C. Simo and M. S. Rifai, ‘A class of mixed assumed strain methods and the method of incompatible modes’, Int. J. Numer. Meth. Engng., 29, 1595—1638 (1990). 6. T. H. H. Pian and K. Sumihara, ‘Rational approach for assumed stress elements’, Int. J. Numer. Meth. Engng., 20, 1685—1695 (1984). 7. K.-Y. Yuan, Y.-S. Huang and T. H. H. Pian, ‘New strategy for assumed stresses for 4-node hybrid stress membrane element’, Int. J. Numer. Meth. Engng., 36, 1747—1763 (1993). 8. J. C. Simo and F. Armero, ‘Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes’, Int. J. Numer. Meth. Engng., 33, 1413—1449 (1992). 9. J. C. Simo, F. Armero and R. L. Taylor, ‘Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems’, Comput. Methods Appl. Mech. Engng., 110, 359—386 (1993). 10. U. Andelfinger and E. Ramm, ‘EAS-elements for two-dimensional, three-dimensional, plate and shell structures and their equivalence to HR-elements’, Int. J. Numer. Meth. Engng., 36, 1311—1337 (1993). 11. M. A. Crisfield, G. F. Moita, G. Jelenic and L. P. R. Lyons, ‘Enhanced lower-order element formulations for large strains’, in D. R. J. Owen and E. Onate (eds.), Computational Plasticity—Fundamentals and Applications, Pineridge Press, Swansea, 1995, pp. 293—320. 12. C. Freischläger and K. Schweizerhof, ‘On a systematic development of trilinear three-dimensional solid elements based on Simo’s enhanced strain formulation’, Int. J. Solids Struct., 33, 2993—3017 (1996). 13. S. Glaser and F. Armero, ‘Recent developments in the formulation of assumed enhanced strain finite elements for finite deformation problems’, Report No. ºCB/SEMM-95/13, University of California, Berkeley, 1995. 14. S. Glaser and F. Armero, ‘On the formulation of enhanced strain finite elements in finite deformations’, Engng. Comput., 1996, submitted. 15. J. Korelc and P. Wriggers, ‘Consistent gradient formulation for a stable enhanced strain method for large deformations’, Engng. Comput., 13, 103—123 (1996). 16. P. Wriggers and J. Korelc, ‘On enhanced strain methods for small and finite deformations of solids’, Comput. Mech., 18, 413—428 (1996). 17. J. C. Nagtegaal and D. D. Fox, ‘Using assumed enhanced strain elements for large compressive deformation’, Int. J. Solids Struct., 33, 3151—3159 (1996). 18. D. Roehl and E. Ramm, ‘Large elasto-plastic finite element analysis of solids and shells with the enhanced assumed strain concept’, Int. J. Solids Struct., 33, 3215—3237 (1996). Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999) SYSTEMATIC CONSTRUCTION OF B-BAR FUNCTIONS 639 19. E. A. de Souza Neto, D. Peric, G. C. Huang and D. R. J. Owen, ‘Remarks on the stability of enhanced strain elements in finite elasticity and elastoplasticity’, in D. R. J. Owen and E. Onate (eds.), Computational Plasticity—Fundamentals and Applications, Pineridge Press, Swansea, 1995, 20. P. Wriggers and S. Reese, ‘A note on enhanced strain methods for large deformations’, Comput. Methods Appl. Mech. Engng., 135, 201—209 (1996). 21. J. Jirousek and L. Guex, ‘The hybrid-Trefftz finite element model and its application to plate bending’, Int. J. Numer. Meth. Engng., 23, 651—693 (1986). 22. J. Jirousek, ‘Hybrid-Trefftz plate bending elements with p-method capabilities’, Int. J. Numer. Meth. Engng., 24, 1367—1393 (1987). 23. J. Jirousek, A. Venkatesh, A. P. Zielinski and H. Rabemanantosa, ‘Comparative study of p-extensions based on conventional assumed displacement and hybrid-Trefftz FE models’, Comput. Struct., 46, 261—278 (1993). 24. J. Jirousek, ‘Improvement of computational efficiency of the 9 DOF triangular hybrid-Trefftz plate bending element’, Int. J. Numer. Meth. Engng., 23, 2167—2168 (1986). 25. A. Venkatesh and J. Jirousek, ‘An improved 9 DOF Hybrid-Trefftz triangular element of plate bending’, Engng. Comput., 4, 207—222 (1987). 26. R. Piltner, ‘Recent developments in the Trefftz method for finite element and boundary element applications’, Adv. Engng. Software, 24, 107—115 (1995). 27. R. Piltner, ‘A quadrilateral hybrid-Trefftz plate bending element for the inclusion of warping based on a threedimensional plate formulation’, Int. J. Numer. Meth. Engng., 33, 387—408 (1992). 28. A. P. Zielinski and O. C. Zienkiewicz, ‘Generalized finite element analysis with T-complete solution functions’, Int. J. Numer. Meth. Engng., 21, 509—528 (1985). 29. M. F. Beatty, ‘Introduction to Nonlinear Elasticity’, in M. M. Carroll and M. A. Hayes (eds.), Nonlinear Effects in Fluids and Solids, Plenum Press, New York/London, 1996, pp. 13—112. 30. M. F. Beatty, ‘Topics in finite elasticity: hyperelasticity of rubber, elastomers, and biological tissues—with examples’, Appl. Mech. Rev., 40, 1699—1734 (1987). 31. J. K. Knowles and Eli Sternberg, ‘On the failure of ellipticity and the emergence of discontinuous deformation gradients in plane finite elastostatics’, J. Elasticity, 8(4), 329—379 (1978). 32. I. C. Taig, ‘Structural analysis by the matrix displacement method’, English Electric Aviation Report No. S017, 1961. 33. T. J. R. Hughes, ¹he Finite Element Method, Prentice-Hall, Englewood Cliffs, N.J., 1987. 34. K.-J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, N.J., 1976. Copyright 1999 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 44, 615—639 (1999)

1/--страниц