- Email: [email protected]
A unified approach to mixed finite element methods: Application to in-plane problems
Computer Methods in Applied Mechanics and Engineering 98 (1992) 127-151 North-Holland CMA 233
A unified approach to mixed finite element methods: Application to in-plane problems Shmuel L. Weissman and Robert L. Taylor Department of Civil Engineering, University of California at Berkeley, Berkeley, CA, USA Received 31 July 1991 Revised manuscript received 11 February 1991 and 16 October 1991
A general method to treat internal constraints within the context of mixed finite element methods has been presented in previous work. The underlying idea is to constrain the assumed stress and strain fields to satisfy the homogeneous equilibrium equations in a weak sense and thus satisfy a priori the internal constraints. This method is now applied to generate four-node plane stress/strain elements. For these elements, it is proved that locking at the nearly incompressible limit (plane strain) is avoided at the element level. The proposed elements are shown to yield excellent results on a set of standard problems. Furthermore, excellent stresses are obtained at the element level.
1. Introduction
1.1. Motivation and objective Over the last two decades much research has been devoted to the design of plane strain elements that do not lock at the nearly incompressible limit. Numerous methods have been proposed (see below). Unfortunately, most of these methods are either based on parameter tuning or cannot be applied to other problems involving internal constraints (e.g., plate bending). In this work, the method presented by Weissman and Taylor [1] is applied to generate assumed stress and strain fields for four-node plane stress/strain elements. In particular, emphasis is put on proving that the proposed elements do not lock at the nearly incompressible limit. In totality, the proposed elements should exhibit the following properties: - r e s o l v e locking at the nearly incompressible limit (plane strain) at the element level; - r e c o v e r stresses that are variationally consistent with the formulation used to obtain the displacement field; -produce reliable stresses at the element level (without recourse to smoothing); -possess correct rank (no spurious energy modes); - h a v e performance independent of coordinate system or user input data.
1.2. Background Solutions of plane strain problems at the nearly incompressible limit have traditionally posed severe problems in the development of finite element models. This setback is the result of the failure of finite element models to impose the incompressibility constraint pointwise
(e.g., [21).
Much effort has been directed at producing finite element models that avoid the locking phenomenon at the nearly incompressible limit. Initial approaches utilized reduced integra0045-7825/92/$05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved
128
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
tion. However, this approach often results in unstable elements due to the appearance of spurious zero energy modes [3]. To overcome this difficulty, Hughes [4] used selective reduced integration (SRI). Malkus and Hughes [5] showed the SRI scheme to fall within the concept of mixed finite element methods. Hughes [6] refined the SRI scheme into a general method, known as the B-bar method. Simo et al. [7] showed that the B-bar method can be derived from the Hu-Washizu variational principle. Using a projection operator, Flanagan and Belytschko [8] showed that the usual finite element approximation can be rewritten in a form that leads to decoupling of the stiffness matrix, which can be written as the sum of the decoupled under-integrated stiffness matrix and a 'stabilization' stiffness matrix. Following this approach, Belytschko and Bachrach [9] presented formulations based upon the Hu-Washizu functional. The best of the elements presented is known as the quintessential bending incompressible (QBI). This element, however, is not frame invariant (i.e., it is dependent upon the coordinate system) unless a unique local coordinate system is defined and the element stiffness matrix and load vector, defined in the local coordinate system, are transformed to the global coordinate system. Wilson et al. [10] took yet another approach They observed that when the four-node element (isoparametric displacement formulation) is subjected to pure bending, it deforms in shear rather than in bending. As a result of this anomaly, the element is stiff in bending. In order to improve the element bending behavior, they introduced a set of incompatible displacements. It was soon realized, however, that the distorted element did not pass the constant strain patch test. Taylor et al. [11] modified the incompatible modes element so that the resulting element passed the patch test. In addition to the improved behavior in bending this element avoids locking at the nearly incompressible limit (this property was not observed, however, when the elements were first presented). Conditions for convergence when incompatible modes are used were presented by Strang and Fix [12]. Pian and Sumihara [13] used incompatible displacements to generate an assumed stress field for a four-node plane stress element, which was formulated via the Hellinger-Reissner variational principle. This element exhibits excellent characteristics in bending applications. In addition, its sensitivity to mesh distortion from a parallelogram shape appears to be among the smallest of four-node elements. Furthermore, when modified to account for plane strain constitutive relations, it avoids locking at the nearly incompressible limit. Initially, the formulation required a quadratic perturbation of the element shape. Recently, Pian and Wu [14] presented a new formulation that avoids this perturbation. A general formulation to generate incompatible element functions was presented by Wu et al. [15]. Punch et al. [16] presented a theoretical framework for the selection of assumed stress fields in the context of elements based upon the Hellinger-Reissner functional. Punch and Atluri [17] presented a family of two- and three-dimensional elements. The four-node element yielded comparable results to the Pian and Sumihara [13] element. The performance of these elements when applied to model nearly incompressible materials was demonstrated by Yang and Atluri [18] for the case of the Navier-Stokes equations. Simo and Rifai [19] presented four-node elements based upon the Hu-Washizu functional which introduced 'enhanced strains'. Results reported for these elements are almost identical to the results obtained by the element presented by Pian and Sumihara [13]. To avoid the limitation principle [20], the assumed stress and strain fields must not contain as a subset the corresponding fields derived from the displacement field. In addition, to avoid locking behavior in the presence of internal constraints, the assumed stress and strain fields must a priori satisfy these constraints pointwise. Weissman and Taylor [1] proposed constructing these assumed fields to possess the aforementioned properties by constraining, in-
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
129
dependently, the assumed stress and strain fields to a priori satisfy the homogeneous equilibrium equations in a weak sense.
1.3. Paper overview The stress and strain field generation procedure is reviewed in Section 2. The strong form for the plane stress and plane strain problem is summarized in Section 3. The finite element approximation is contained in Section 4. In Section 5, four-node quadrilateral elements formulated via either the Hellinger-Reissner or the Hu-Washizu functionals are proposed for both the plane stress and plane strain problems. It is proved that in the case of plane strain the trace of the assumed strain vanishes, pointwise, as the nearly incompressible limit is approached. Consequently, locking is avoided. Numerical examples demonstrating the excellent properties of the proposed elements are presented in Section 6. Concluding remarks are given in Section 7. Before continuing to Section 2, there is an important point to be made. The equivalence between the formulation based on the Hu-Washizu functional and that based on the Heilinger-Reissner functional, as evident from the numerical examples in Section 6, holds only in the case of linear elasticity. When the constitutive equations are nonlinear, the derivations of a consistent tangent stiffness rely crucially upon the use of the Hu-Washizu functional [7], while the Hellinger-Reissner functional leads to restrictive (cannot model perfect plasticity) and cumbersome algorithms [21]. Consequently, it is the Hu-Washizu functional that is of interest if the methodology presented is subsequently to be used in nonlinear problems.
2. Summary of assumed fields generation procedure Internal constraints pose a difficult problem for finite element methods, commonly manifested by locking. This difficulty was put in the form of a limitation principle by Fraeijs de Veubeke [20]. In order to avoid this limitation principle, the assumed stress and strain fields must satisfy a priori the equilibrium equations, and thus effectively satisfy the internal constraints. In the context of finite element methods this requirement can be relaxed, and the equilibrium equations can be satisfied in a weak sense. Taking advantage of this observation, Weissman and Taylor [1] presented a general method to generate assumed stress and strain fields. These fields are constructed to satisfy in a weak sense the homogeneous equilibrium equations, and consequently the internal constraints. This method is formulated within the framework provided by the Hu-Washizu variational principle, and is summarized below. Before proceeding to describe the proposed method it is useful to recall the steps involved in deriving finite elements via the Hu-Washizu functional. To this end let e, o" and U denote the assumed strain, stress and displacement fields, respectively; L the strain displacement operator; and D the elastic coefficients matrix. The Hu-Washizu functional, for the case of small deformations and linear elastic materials (ignoring initial stresses and strains) is given by HH(Cr, e,
U):= fa [½etDe+ (rt(LU- e)]da-//~xT(U),
(2.1)
where g2 is the domain, and //EXT is the external work. The weak, or variational, form may be obtained from the energy functional,//H, by making
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
130
it stationary. This result is obtained by taking the first total variation of H H and equating it to zero. To this end, let ~e, ~o" and 5U denote virtual strain, stress and displacement fields, respectively. The weak form is given by
DII n := ~
[~et(oe
-
o') +
~o't(LU - e) + (L ~U)'o" - ~Utb] d O -
fc
5U'o'" d F = O,
t
(2.2) where b is the body force vector, o-" are applied traction boundary conditions, and C, is the part of the boundary on which traction boundary conditions are specified. The finite element model is obtained by approximating the assumed fields. To this end, let the assumed fields be approximated over each element domain by - Stress field: o- = Ss, - Strain field: e = Ee, - Displacement field: U = Nd, where S, E and N are the shape functions for the approximating stress, strain and displacement fields, respectively; and s, e and d are the vectors of independent stress, strain, and displacement parameters, respectively. Substituting the assumed fields into the weak form leads to the following system of equations (these are the approximated Euler-Lagrange equations):
-A
0
0
G t
where
= 0
(2.3)
,
d
A " = f, StE d O ,
n "= f, EtDE dO , f G ' = I StB dO
.It
f
"= fg~ Ntb dO + fc 't Nttra d F
(2.4a)
'
(2.4b)
and B is defined by B : = LN. Eliminating the ~tress and strain coefficients yields
Kd=f ,
(2.5)
where K is the element stiffness matrix, given by X=
.
(2.6)
To satisfy stability conditions for every admissible displacement, the following inequalities must hold [22]: n~+n d~n,,
n,~>n d,
(2.7)
where n~, n~ and n~ are the number of strain, stress and displacement parameters, respectively (n d is equal to the number of nodal degrees-of-freedom). ~ ~The stability conditions presented in (2.7) are necessary but not sufficient. In addition some constraints must be placed on the admissible stress and strain shape functions S and E, respectively. In particular these functions should be chosen so that the LBB condition is satisfied (see e.g. [23, 24]).
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
131
A significant reduction in the computational effort can be obtained if A is invertible. Henceforth, A is assumed to be invertible. 2 As a result n~ = n~ and, consequently, the mixed patch test [22] provides the minimal number of independent stress and strain parameters. Furthermore, in view of the desire to minimize computations, this is the optimal number. As was stated above, the underlying idea behind the proposed method is to construct assumed stress and strain fields so that they a priori satisfy the homogeneous equilibrium equations in a weak sense. To this end, let the assumed strain field be given as the sum of two independent fields as follows: e := g': + e i
(2.8)
and let, for convenience, e ~ be given by i
e "=LU
i
(2.9)
where U i is a displacement field not contained in U h, which is the finite element approximation of the admissible displacement solution space. The initially assumed stress, if, and compatible strain, F c, fields are constrained so that the contribution to the energy functional, H H, from the mixed terms containing e i and either U or ff vanishes in a weak sense over each element's domain. As a result, the reduced stress, or, and compatible strain, e~ satisfy the homogeneous equilibrium equations in a weak sense. Furthermore, it follows from the variational structure that the incompatible strain, e ~, is zero pointwise (at the solution). The proposed method is summarized in Box 1. The condition U ~ ¢ U h (in Step 2) is introduced in order to prevent the appearance of spurious zero energy nodes (see [1] for discussion). The constraint 3(a) was introduced by Strang and Fix [12] as a requirement to pass the constant stress/strain patch test. Constraints 3(b) and 3(c) arise from the variational statement of the problem. Requirements similar to 3(b) were used by Pian and co-workers (e.g., [13, 14]). However, while the objective was the same (i.e., to constrain the assumed stress field to a priori satisfy the homogeneous equilibrium equations in a weak sense), the constraint equation itself was different, and motivated by different arguments. Constraint 3(c) was also presented in the context of the free Box 1. Proposed method 1. Select initial assumed (enriched) strain field, U, and stress field, (Pascal triangle). 2. Select the incompatible strain field, e i ( U i ~ U h). 3. Reduce the number of independent variables in U and ~ by introducing the following constraints: (a)
f ~ U .,i d ~ = 0 ,
(b) (c)
fl ei'D~ c dg~ = 0
--~ e" .
4. Substitute the reduced fields into the variational principle.
2
This assumption crucially depends on the selection of S and E ,zo that the LBB condition is satisfied.
S.L. Weissman, R.L. Taylor, A unified approach to m&ed FEMs
132
formulation (e.g. [25]). However, this condition was introduced by Bergan et al. not as a constraint on the assumed compatible strain field, e c, but rather as a constraint on the admissible higher order modes (i.e., incompatible displacements).
REMARK. Recall
that the Hellinger-Reissner functional can be obtained from the H u Washizu functional by enforcing the stress/strain relationship (i.e., tr = pointwise. If this relationship is substituted directly into the Hu-Washizu functional, e i drops out and, consequently, the proposed method cannot be used. However, if the formulation is started with a Hu-Washizu functional and the constraints enforced there, while taking notice of the special choice for the strain field, then the same reduced stress field is obtained and e i drops out only after obtaining the desired stress field.
De)
3. Strong form The strong, or classical form for both the plane stress and plane Strain pa-oblems is presented. Throughout this section Greek subscripts take the values 1 and 2. Repeated indices imply the usual summation convention. All quantities are referred to a fixed system of rectangular, Cartesian coordinates. A general point in this system is denoted by (x~, x 2, %). The domain of interest is a three-dimensional body embedded in an Euclidian three-space R 3. In the undeformed configuration, the domain O is of the following form:
{
I
J'~ "-" (X I , X2, X3) {~ R 3 X3 E
[-h(xn'x2) 2} 2 ' +h(xI'x2)],(Xn,X2)~ACR 2
where h is the thickness in the x 3 direction ~ and A is a closed region bounded by a simple closed curve C with an outgoing unit normal v. Let C U and C, be subregions of C such that C U U C, = C and C u N C, = t~, where C u is the part of C on which displacements are specified and C, is the part of C on which tractions are specified. The strong form is summarized in Box 2 where b,, are the body forces, o-~ are the specified traction boundary conditions, U~ are the specified displacement boundary conditions and C , , ~ is the elastic coefficients tensor, which for the case of plane strain is given by C,~
:= hI/.~(3,~t$~ + ~ )
+ aS,~tSvsl,
(3.1)
where a and /z are the Lam6 coefficients, which can be eliminated in favor of Young's Box 2. Plane stress/strain: strong form Given b , , o'i', and U~',; find U,, and tr~ such that
in A
Io',,~,t~ + b = 0 , { o',~ = C,,~,t e,~ ,
l~.,~-- ~(u,,,~ + on ct, on
C,
U~.,,);
u,, = u"" a
vt3tr,v = or, .
3 In the case of plane strain h - 1.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
133
modulus, E, and Poisson's ratio, v, using vE A= (1 + v ) ( 1 - 2v) '
E /z = 2(1 + v) "
(3.2)
The plane stress elastic coefficients are obtained by replacing A by J(, where ;( is given by 2/xA = vE ' ( = A +2/x ------5 1-v "
(3.3)
Note that in the case of plane stress, the strain terms in the x 3 direction can be obtained from the constitutive equations. Similarly, in the case of plane strain, the stress terms in the x 3 directions can also be obtained from the constitutive equations. Consequently, the in-plane stress and strain fields characterize the three-dimensional state for both problems. It follows from Box 2 that the strain displacement operator L is given by "
b
oxz L:=
0
0
(3.4)
Ox 2
0 b Ox2 Oxl-
s
0]
and the elastic coefficient matrix D is given by
D= 1-V
2
(3.5a)
v 1 0 0 0 ½(1 -- V)
for the case of plane stress, and by
E
D = (1 + v ) ( 1 - 2 v )
I, 0v vv1 -0 v
° 1
0 ~(1-2v)
(3.5b)
for the case of plane strain.
4. Finite element approximation The method summarized in Section 2 is used to formulate four-node quadrilateral plane stress and plane strain elements. The initially assumed stress and strain fields as well as the incompatible displacements used to generate the assumed incompatible strains are presented. In Section 4.1, the assumed displacement field is introduced. The assumed stress field is given in Section 4.2, and the assumed compatible strains in Section 4.3. The incompatible displacements used to generate the assumed incompatible strains are summarized in Section 4.4. 4.1. Assumed displacement field
The assumed displacement field is taken as a standard isoparametric field. The displacements are approximated by
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
134
(4.1)
U = Nld I ,
where N t is the shape function associated with node I, and d I is the nodal displacement vector• In the case of four-node quadrilateral elements, the shape functions are given by (4.2)
Nt( ~, 7) = -~(1 + ~:t~)(1 + n, 7/),
where ~: and ~ are the element natural coordinates, on the interval [ - 1 , 1], and ~1 and r h are the values of the natural coordinates at node I. 4.2. A s s u m e d stress fields
The assumed in-plane stress field is a complete linear field 4 expressed in the element natural coordinates as
t~*-
~~ tT,~ -
1
~ 0 0 0 0 0
I *l 0' 0, [s,;j 1
0
o 0 o
0
0 lrs~ 0 s*
l! ]
(4.3)
e " LZ
Since complete polynomials are used to express the in-plane stresses, (4.3) could be used for ~ directly. However, the reduction to satisfy the constraints introduced in Box 1 would require selection of different parameters in each element (i.e., there would be a dependence on element orientation of ~ and 77with respect to x~ and x2). This dependence may be avoided by using the transformation procedure described below. The following definitions are introduced: xs= ~txtt,
x t = ~rltxt,,
xh=
¼(serl)txt,,
ys = ¼~,xai,
y t = ¼rhx2,,
y h = ~( grl)tx2t.
Following Zienkiewicz and Taylor [31], the Jacobian of the coordinate transformation from the (~, 7/) space to the (x~, x,) space is given by
J = Jo +
+ J2n,
where Jo = xs y t -
xt y s ,
(4.4) (4.5a)
Ji = xs y h - x h ys ,
(4.5b)
J2 = x h y t -
(4.5c)
xt yh .
The stress field in the physical space is obtained by using the following transformation:
1 crq = ~ Fi,FkKtr lr ,
(4.6)
4 The stress resultant and strain fields are assumed as complete linear fields since this is the best assumption that may be used in conjunction with a bilinear displacement field.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
135
where both i and j take the values x, or x 2, and both I and J take the values ~ or v, and Fxl ~ = Ox~/O~, etc. However, to maintain the ability to model constant fields (i.e., pass the constant strain patch test), the transformation is performed based on the Jacobian, J, and 'deformation gradient', F, at the center of the element. At the center of the element, F is given by Fx, ~ = xs ,
Fx,, = xt ,
(4.7a)
Fx._e = ys ,
Fx., = y t .
(4.7b)
After redefining the independent coefficients, the assumed stress field is given by
=
I°112[!oo [ =
tr~2J
1
0
0
1 xs ysrt
xt2, x2,
yt2rt
2 ys yt ~
2 ys yt71
xt yt ~ xs ys ~
xt ytrt
A(
Art
yt z ~
ys2rt
2
ys z ~
sx,
2xxt ][iIs2 9
(4.8)
where A = xs yt + xt y s , 1
S I ~'- ~
1
S2 = ~
,
,
( y s 2 S~ + yl z S 4 + 2 ys y t s 7 ) ,
1
s3=~(xsyss~+xtyts 1
Si=~Sj
,
(XS 2 S 1 + x t 2 S 4 "~- 2 X S X t S 7 ) ,
,
,
4 +As 7),
,
,
with (i, j ) E {(4, 3), (5, 5), (6, 2), (7, 6), (8, 8), (9, 9)}. It is convenient to write or in the following form:
where S,,t2 ._ ( $ 6 , S7, $8, $9 ) .
With the above construction, the parameter set g2 may always be selected as the set to be eliminated in satisfying the constraint equation on the stress field (see Box 1). 4.3. A s s u m e d strain fields
The assumed strain field is formulated in the element natural coordinates, and then transformed into the physical domain. Two types of transformations can be used: - t h e same transformation as used for the corresponding stress fields; - t h e inverse to the transformation used for the corresponding stress fields. The motivation for the first approach is to have the same 'shape functions', or interpolation, for the strain fields as for the corresponding stress fields. The second approach is motivated by
s.L. Weissman.R.L. Taylor. A unified approach to mixed FEMs
136
the invariance of the complementary energy (trijeij) under coordinate transformation. Accordingly, the transformation is given by 5
eo = JF~' F~' e,,*.
(4.9)
A similar approach was taken by Simo and Rifai [19]. Following the path established for the assumed stress field, the strain field is assumed as a complete linear field in the element natural coordinates, and then transformed into the physical space• Using the transformation given by (4.6) yields e4 , El2
[!
xsErl ys2~
0 0
=
1 0
0 1 xsysrl
xt2~
xs2~ ys 2~
xtErl yt2~
2xsxt~ 2ysyt~
xtyt~ xsys~
xtyt,1
A~
yt 2~:
2xsxt'o][i']e2 2ysyt~ (4.10a) Arl
and using the inverse transformation, (4.9), yields
i ~'C
E
0 0 yt 2 71 1 0 xtErl 0 1 -2xtyt77
=
X
ys2 + xs z ~ -2xsys~
yt 2 ~ xt 2 t~ -2xtyt~
ys2 ,1 xs2rl -2xsys~l
Ill e., I-
-ys yt + -ys yt,ll -xsxt ~ -xsxtrl I" A~ A~q J (4.10b)
•
tj
It is convenient to write e in the following form:
[A]
+: = [a, I a..,l e,
where
-
e2
'
At
el = (el, e2, e3, e4, es) ,
~ = ( e 6 , e7, e~, e g ) .
With the above construction, the parameter set ~2 may always be selected as the set to be eliminated in satisfying the constraint equation on the strain field (see Box 1).
4.4• Assumed incompatible displacements Let the assumed incompatible displacements be given by •
"
i
Ult = N'tAt + NEA2
(4.11a)
Ui2 = NilA3 + Ni2A4 ,
(4.11b)
and where At, '~2,. • •, and Aa are the independent incompatible displacement parameters, and Ni,~ are the incompatible shape functions. Two options are selected for the incompatible shape functions Nit and Ni2. The first set was The transformation is based on the values at the center of the element.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
137
"presented by Wu et al. [15]. The functions are given by •
N'~ = se2 and
2J2 3,/o ~: + ~oo r/ 2J l
(4.12a)
2J I 2J2 N~ = "r/2+ ~ ~ - 3J,--~7/.
(4.12b)
The second set was presented by Taylor et al. [26]. In order to obtain a more compact form, the incompatible modes given below are a linear combination of the modes originally presented. The functions are given by N , - ( 1 - ~ J2 n)(1 - ~:2) + ~J,: ( 1 - T I
2)
(4.13a)
and N'2 = ( 1 - J' ~ : ) ( 1 - 2 ) + J , f~ 71(1- ~:2).
(4.13b)
Note that the second set of incompatible modes is zero (compatible) at the nodal points while the first set is not.
5. Proposed elements
The fields presented in Section 4 are used to generate four-node plane stress/strain elements. The proposed elements are cataloged in Table 1 by the type of formulation, the incompatible shape functions, and assumed strain field (Hu-Washizu). For all these elements, the assumed displacement and stress fields are given by (5.1) and (4.8), respectively. The remainder of this section consists of proofs (for the case of plane strain) showing that as the incompressible limit is approached, the assumed strain vanishes pointwise for the proposed elements. Consequently, the well-known locking behavior is avoided at the element level. The proof for elements formulated via the Hellinger-Reissner and Hu-Washizu variational principles are given in Propositions 5.1 and 5.2, respectively. P R O P O S I T I O N 5.1. The Hellinger-Reissner elements presented above can, in the case o f plane strain, model a state in which the trace of the strain goes to zero, pointwise, as v goes to
0.5-.
Table 1 Proposed plane stress/strain elements Element
FormuNtion
Incompatible shape function equations
PSS1 PSS2 PSS3 PSS4 PSS5 PSS6
Hellinger-Reissner Hellinger-Reissner Hu-Washizu Hu-Washizu Hu-Washizu Hu-Washizu
(4.12) (4.13) (4.12) (4.13) (4.12) (4.13)
Strain field representation -(4. l 0a) (4.10a) (4.10b) (4.10b)
138
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
PROOF. Recall that in the case of the Hellinger-Reissner variational principle the assumed strain field is obtained from the assumed independent stress field. Consequently, in the case of plane strain the strain field e is given by
e
=
D-lo"
(5.1)
where D -l is given by
o
-l = _
A
a
o 0] 2/z
0
(5.2)
,
allx
with a = 4/z(/~ + A). As Poisson's ratio, v, goes to 0.5-, A goes to infinity while/z remains bounded. It follows that, at the limit, the following relation is obtained: --
D 111=
--|
D22
--
- D l,,l = -D2t I .
-"
(5.3)
After imposing the constraint equations introduced in Section 3, the assumed stress field presented in (4.8) is reduced to the general form
0 a2r/+b2G 1 a3'r/ + bag
1
LO'12J
0
asG+b~,H/.~/, a6G +
b6nJLssj
(5.4)
where a t, a 2 , . . . , and at,, and bl, b 2 , . . . , and b 6 are constants to be determined by the procedure presented in Box 1. Substituting (5.2) and (5.4) into (5.1) and taking notice of (5.3) yields the desired result. [] P R O P O S I T I O N 5.2. The Hu-Washizu elements presented above can, in the case o f plane strain, model a state in which the trace of the strain field goes to zero, pointwise, as the incompressible limit is approached. PROOF. After constraining the assumed strain field introduced in (4.10), the reduced field is given in a general form by
Jell]
e =/ezz/= Lel2]
al rl + blG a 2rl + b2G asG + b 5 1 a 3rl + bag a6G +
[1
0
0
[0
1
0
0 0
b6rt.lL/~s
(5.5)
The structure imposed by the constraint equations is such that the following relations are obtained, as v goes to 0.5-, between the at, a 2 , . . . , and a 6 and b l, b 2 , . . . , b 6 coefficients: a~ + a 2 = 0 ,
a4 + a5 = 0
(5.6a)
b I + b2 =0,
b4 + b5 =0.
(5.6b)
and
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
139
A formal proof for (5.6) can be obtained. This proof, however, is very cumbersome and tedious, and does not yield a better understanding of the formulation. For these reasons, a numerical 'proof' will be given. To obtain a more compact form, an equivalent form to (5.6) is given. Namely, the values of (a~ + aE)r/+ (b I + bE)g and (a 4 + as)~ + (b4 + bs)r/are reported at the Gauss points. Since the functions are linear, zero value at all Gauss points implies zero value pointwise. First, a square element (see Fig. l(a)) is examined. The values for the first Gauss point are reported in Table 2. Values for other Gauss points may be obtained as follows: given that the value at the first point is o, then at the other points it is -sign(C) x v, where sign(~) is the sign of the ~: or r/coordinate at the given Gauss point (7/for the first column and ~ for the second column). Secondly, an element distorted from a parallelogram shape (see Fig. l(b)) is examined. In this case the values at the first and second Gauss points are reported. The value for the third Gauss point is minus the value at the second Gauss point, and the value at the fourth Gauss point is minus the value at the second Gauss point. Results are summarized in Tables 3 and 4. The results summarized in Tables 1, 2 and 3 'prove' (5.6). It follows from (4.12) and (4.13) that as v goes to 0.5- the trace of e is given by
(5.7)
trace(e) = e I + e : .
In order to show that the trace of the strain goes to zero pointwise recall from Section 2 that the stress strair~ weak relation is given by He
-
As
X2
(5.8)
= O .
X2
11
' (o.1)
(0.25,1.25)
(1,1)
(1.2S,0.5)
(0,0)
(~ ,o) X1
X1
('o) Fig. 1. (a) Square element; (b) element distorted from a parallelogram shape.
Table 2 Square element
0 0.3 0.49 0.499 0.4999
(a I + a2)"O + (b I + b2)~
(a 4 + as)~ + (b4 + b5)7/
-0.144338 -0.082479 -0.005661 -0.000577 -0.000058
-0.144338 -0.082479 -0.005661 -0.000577 -0.000058
140
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
Table 3 Distorted element (a, + a2)'q + (b~ + b2)~
v
Point
PSS3
PSS4
PSS5
PSS6
0
1 2
-0.138299 -0.142666
-0.136318 -0.145547
-0.239751 -0.234665
-0.166003 -0.224427
0.3
1 2
-0.074126 -0.077283
-0.072729 -0.079388
-0.2255061 -0.204323
-0.231639 -0.169421
0.49
1 2
-0.004723 -0.004982
-0.004612 -0.005156
-0.051700 -0.055967
-0.037454 -0.042924
0.499
1 2
-0.000479 -0.000505
-0.000467 -0.000523
-0.006410 -0.006890
-0.004384 -0.004980
0.4999
1 2
-0.000048 -0.000051
-0.00U047 -0.000052
-0.000650 -0.000700
-0.000446 -0.000506
Table 4 Distorted element (a 4 + a~)~j + (b4 + bs)~ v
Point
PSS3
PSS4
PSS5
PSS6
0
1 2
-0.099318 +0.102311
-0.096704 +0.103024
-0.386263 +0.388999
-0.414609 +0.420426
0.3
1 2
-0.053029 +0.055193
-0.051169 +0.055728
-0.692734 +0.703892
-1.013897 +1.045209
(I,49
1 2
-0.003365 +0,003542
-0.003215 +0.003588
+0.027529 -0.029822
+{I.021353 -0.024105
0.499
i 2
-0.000341 +0.000359
-0.000326 +0.000364
+0.002510 -0.002771
+0.001978 -0.002279
0.4999
1 2
-0.000034 +0.000036
-0.000033 +0.000036
+0.000249 -0.000276
+0.000197 -0.000227
To simplify notations and without loss of generality, ~ and 7/in the a s s u m e d strain field, eq. (5.5), may be replaced by f and ~, given by
Go,
h = n - rto,
w h e r e ~0 and r/0 are shifts introduced in order to m a k e H " block diagonal. In order to o b t a i n this result the following r e q u i r e m e n t s must be met: f,,t ~ d A = O '
fA'~dA=O"
T h e resulting Go and 7/o are given by
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
J, 3Jo'
141
7/° = 3J0"
The H matrix is now of the following form: H=
[o0a /~o] ,
(5.9)
where D is given by (3.5), A is the element area, a n d / t is the 'material' properties for the linear part of the strain field. Based on results in (5.6), /~i (~'= {e 4, es} ) goes to zero pointwise as the incompressible limit is approached. Hence, it remains only to consider the constant part, given by D A . Combining (5.8) and (5.9) yields 1 e, = a A
and
1 e2 = a A
{(A + 2/x)§,- Ag2}
(5.10a)
{(A + 2/~)g2 - Ag,},
(5.10b)
where g = A s . Substituting (5.10) into (5.7) yields trace(e) =
2/x(~, + s2) aA
(5.11)
"
As v goes to 0.5-, a goes to infinity while/~ is bounded; hence the desired result is obtained.
6. Numerical examples
The performance of the plane stress/strain (PSS) elements proposed in Section 5 is evaluated with several discriminating problems selected from the literature. The purpose of these evaluations is to test the proposed formulation's sensitivity to the specific choice of the incompatible shape functions as well as the overall performance of the proposed elements. First tackled are the constant strain patch tests in Section 6.1. Following these are the beam problems suggested by MacNeal and Harder [27] in Section 6.2, sensitivity to mesh distortion in Section 6.3, the Cook membrane [28] in Section 6.4, the thick walled cylinder also introduced by MacNeal and Harder [27] in Section 6.5, and a strip of finite width with a hole in Section 6.6. Convergence of the results obtained by the four-node elements presented in this paper are compared with other well-known four-node plane stress/strain elements. A listing of these elements, and the abbreviations used to identify them is given in Table 5. In all the examples given in this paper, identical results are obtained for elements Table 5 Reference elements Element
Description
2x2 P-S QBI
Isoparametric plane stress/strain Pian and Sumihara [13], plane stress/strain Belytschko and Bachrach [9], plane stress/strain
142
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
formulated via the Hellinger-Reissner variational principle and the corresponding (i.e., same incompatible shape functions) elements formulated via the Hu-Washizu variational principle. Identical results are also obtained for the two approaches taken to transform the assumed strain field from the element's natural space into the physical space. Convergence is examined in the energy norm (energy norm reported is twice the internal strain energy) since this is the natural convergence test in the context of finite element methods [12]. Convergence of displacement at some characteristic points is also examined as this is the common practice in the literature. All tables and figures show the displacement/ energy norm as a function of the number of elements (denoted nel) used in the corresponding mesh. 6.1. Patch test
A rectangular domain is modeled by a single element shown in Fig. 2(a) and the skewed mesh shown in Fig. 2(b). The mesh is subjected to a constant state of tension/compression. All elements presented pass this test. 6.2. Beam problems: in-plane shear
These problems were suggested by MacNeal and Harder [27] as standard problems to evaluate the performance of different elements. The meshes used contain only one row of six elements for both the straight beams, shown in Fig. 3, and the curved beam, shown in Fig. 4. Geometrical and material properties used are summarized in Table 6. The 'exact' tip displacement, w, given by beam theory [27] for the case of a straight beam is w = 0.1081 and for the case of a curved beam is w = 0.08734. The results, normalized with respect to the exact solution, are summarized in Table 7. Results obtained for all proposed elements are identical to the results obtained by the Pian
x2
x
(0,
,12)
n~
2 10.00,0.12)
(0.24,0.12)
~0.08.0.08)
(0.24,0.00)~ ' x
(0.00,0.00) (a)
•~0.16.0.0B)I
" (0.~,0.00)
I, ~'X (0.24,0.00)
(b) Fig. 2. Patch test: (a) one-element mesh; (b) skewed mesh.
Table 6 Material and geometrical properties for beam problems
Straight b e a m Curved b e a m
E 1.0E+7 1.0E+7
v 0.3 0.25
Thickness 0.1 0.1
Length/ inner radius 6.0 4.12
Width / outer radius 0.2 4.32
Arc m 90°
143
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
i /
i \
J
L,,o%
(z)
(b)
Fig. 3. Straight cantilever beam. (a) Regular shape elements; (b) trapezoid shape elements; (c) parallelogram shape elements.
Table 7 In-plane shear
LqDqLWJW.
Mesh
PSS/P-S
2x2
QBI
Straight (a) Straight (b) Straight (c)
0.9929 0.2208 0.7963
0.0933 0.0278 0.0348
0.9929 0.0777 0.0890
Curved
0.9782
0.0740
0.2324
ql, W . ' 4
Fig. 4. Curved beam mesh.
and Sumihara [13] element. The 2 x 2 element locks for all meshes. However, if two rows of elements are used, the 2 x 2 element yields reasonable results. All elements exhibit locking when trapezoid-shape elements are used (straight mesh (b)). This is in accordance with the theorem put forth by MacNeal [29] regarding the locking of tapered four-node membrane elements. The QBI element appears to be much more sensitive to mesh distortion than the PSS or P-S elements.
6.3. Beam bending: sensitivity to mesh distortion In th~s standard test, a beam is modeled by two elements. The beam is fixed at one end and subjected to a bending moment on the other, as shown in Fig. 5. The material properties used
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
144
10.0
:::::::::::::::::::::::::
:::::::::::::::::::::::::::
-F
_.__.._.~
Fig. 5. Beam bending problem, sensitivity to mesh distortion.
are E = 1.0 and p = 0.4999. A state of plane strain is assumed. The edge separating the two elements is gradually rotated, as shown in the figure, to skew the mesh. Results, normalized with the exact beam theory solution (w - 562.575), are shown in Fig. 6. Results for the 2 x 2 element are not reported since it exhibits severe locking. All proposed elements yield identical results to the Pian and Sumihara [13] element. The QBI element exhibits rapid deterioration as A is increased.
6.4. Cook's membrane problem The problem consists of a trapezoidal plate clamped on one end and loaded by a uniformly distributed in-plane bending load on the other end, as shown in Fig. 7. This problem has a considerable amount of shear deformation. It was suggested by Cook [28] as an excellent problem to test membrane elements using skewed meshes. No analytical solution is available for this problem. The best known solution is given by Cook as V ¢ - 0.29 (vertical displacement), O"A 0.236 (maximum principal stress) and o-n = - 0 . 2 0 1 (minimum principal stress). The material properties used are E = 1.0, v = 1.0/3.0, thickness = 1.0; plane stress assumption was used. The results normalized with respect to these best known answers are shown in Figs. 8 (Vc), 9 (trA) and 10 (o'~3). The results presented for one element are obtained by interpolating the nodal values; for the finer meshes the nodal values are reported. The results obtained for all proposed elements =
1.5
m , . ,m
PSS/P-S
=="='=
~1 EXACT
E
1.0
o ¢1 .m m O,, r~
N Q=m
o o
¢1 ¢1
E t.. O
0.5
0
i,=
0
4
QOOOQQ~IIUQ
B U N ii iI U iN II ii Ii D a B U l l
0.0
•
0
;
j QB I n i m I B u M
|
3
|
,
s
A
Fig. 6. Beam bending, sensitivity to mesh distortion: normalized tip displacement.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
c
u
m b
145
ti.o 16.0 row,
44.0
44.0 d
:~.:.::-" .,... I[:i:~:i:!?i:~:i...!:~:
ti!~iii~!i!~ii:i i ... i.. .
q
...
.... I[:i!i:!.i:i:!::i}:::ii!i. ..... lt!!~!i :i~!i!:!::i:i
48.0
I[ !! ili: -
Fig. 7. Cook's membrane problem. ¢¢
1.25
cU
E CU c~
1.00
m
D. r~
em
"O
S
0.75
m
..o~
¢1
em
I_ CU
o° #
0.50 ,
o°~ ~ o°°;#
cU N
emil m
,-~.m
PSS/P-S
~mm
2X2
0.25
EXACT
M
E I,, 0.00
1'
;
,'s
6'4
nel
Fig. 8. Cook's membrane problem; vertical disp'~acement at point C.
yield identical results to the Pian and Sumihara [13] element. Excellent results are obtained even for very coarse meshes; about 70% of the vertical displacement and 80% of the stresses are obtained using only one element. The 2 x 2 and QBI elements yield poor results for coarse meshes; less than 50% and 60% of the vertical displacement and stresses are obtained by the 2 x 2 and QBI elements, respectively, when a four-element mesh is used.
6.5. Thick-walled cylinder This problem was suggested by MacNeal and Harder [27] as an excellent problem to test the effect of nearly incompressible material. The mesh used is shown in Fig. 11. The material
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
146
1.25
E E
1.00
e.m
M m
E
f..
o.zs
"Or~ . ag . ~"
•
0.50
-
- - .
@mm
""
"-- "
2X2
IN IN II ID II IN,
EXACT 0.00
I
I
1
4
l's
I
s'4
25s
nel
Fig. 9. Cook's membrane problem; maximum principal stress at point A. 1.25
E E
1.oo,
.,.~1;;~ "~-
•E E ~~
.,*
0.7S'
L..
~o
•, .
oO'~;"
S 4p
. m.
J
m
.-
•E
~
0,25
,,, "-
~ ~ ~ ~
~
,
0.00
-
PSSIP-S 2X2
"''''"
(~1
~
EXACT
,
,
,
le
64
256
nel
Fig. 10. Cook's membrane problem; minimum principal stress at point B.
.0 °
P = I.@ 3:0
3.$
,I.2
S.2
6.7S
9.0
Fig, 11. Thick-wailed cylinder problem.
properties are E = 1.0 and u varies from 0.0 to 0.4999. The exact radial solution is given by
u = E ( R ~ - R2,) T where p is the internal pressure; R~ is the inner radius; R• is the outer radius; and r is the radius where the radial displacement, u, is to be computed.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
147
Table 8 Thick-walled cylinder v
PSS/P-S
2x 2
OB!
0 0.3 O.49 O.499 O.4999
11.9987 0.'.~9~4 0.9913 0.9910 0.9910
11.9964 11.99118 0. 8490 11.36116 O.0534
0.9985 11.9952 0.9911 O.99119 O.99118
Results, normalized with respect to the exact solution, are summarized in Table 8. Once again all proposed plane elements yield identical results with the Pian and Sumihara [13] element. Almost no deterioration with v is observed, only 0.77% difference between the two extreme values of v. The QBI element yields almost identical results. The 2 x 2 element exhibits the well known locking at the nearly incompressible limit. 6.6. Finite strip with a hole This problem is introduced in order to test the accuracy of the stress distribution inside the proposed plane elements. A strip of finite width with a hole, of width b, is subjected to uniform tensile stress in the axial direction, as shown in Fig. 12. The axial stress along the critical section, line A - B , is investigated. First, the relatively simple case of plane stress with free upper and lower edges is investigated in Section 6.6.1. Then, the more challenging case of plane strain, at the nearly incompressible limit, with the upper and lower edges fixed in the vertical direction, is tackled in Section 6.6.2. 6.6. I. Finite strip with a hole: plane stress The material properties used are E = 1.0 and v = 0.3. First, convergence in the energy norm is examined. No analytical solution for the energy is known. Therefore, a converged finite element solution of 19.37 obtained for the Pian and Sumihara [13] element using 3072 elements is used as the reference solution. Results normalized with respect to this reference solution are shown in Fig. 13. The axial stress distribution along the critical section for the meshes containing 48 and 192 elements is presented in Fig. 14. The results are compared with a reference solution given by Savin [30]. 5.0
L I"
,t
Symmetry or Free .~-
B
tu
b -- 10/3
~
P = 1.0
A
Symmetry Fig. 12. Strip with a hole problem.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
148
1.01
1.00
0.99 ¢¢
0.98
~L~ ~
"¢1
-'"
' elm m
0.97
-*I
E Z
:/
0.96
0.95
""
PSS~-S
- - .-. u
2X2
''"'"
0131
-""--
EXACT
/ I
,
3
; 1
4
'
,
192
768
, 3072
nei
Fig. 13. Strip with a hole; plane stress; convergence in the energy norm. 3.5 I= O
omm
e,m
z,.. r,~
2.5 ' ~
imm
- '
EXACT
"'
PSS/48
~'
r~ z.. rj:
1,5
Ilmm
M ,<
0,5
0:4
o'.s
o'.s
o.r
0.8
0:9
',.o
x2..Ib
Fig. 14. Strip with a hole; plane stress; axial stress distribution along the line A-B.
Almost identical results are obtained for the proposed elements and the Pian and Sumihara element (about 5.0E-3 difference in energy norm for the three-element mesh). A very small sensitivity to the incompatible shape functions used is observed (about 2 . 0 E - 3 difference in the energy norm for the thee-element mesh). Excellent results are obtained for both the energy norm, 95.82% of the reference solution with only three elements, and stress distribution inside the elements. The 2 x 2 and QBI elements yield comparable results for this problem. 6.6.2. Finite strip with a hole: plane strain The material properties used are E = 1.0 and v = 0.4999. First, convergence in the energy norm is examined. No analytical solution for the energy is available. Therefore, a converged finite element solution of 2.389 obtained for the Pian and Sumihara [13] element using 12288 elements is used as the reference solution. Results normalized with respect to this reference solution are shown in Fig. 15. The axial stress distribution along the critical section for the meshes containing 48 and 192 elements is presented in Fig. 16. No analytical solution for the
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
149
1.25
1.00 "
I...
0.75
¢: (1#
m
¢!,)
"
0.50 ~
N
em iml
,-,,
E I,,,
0.25
Z: ;' . . . . . . 3
000
m
PSS2,4,6
m ,=.m ,
0131
='''"='
2X2
T .... 12
PSS1,3,5
= "
•, . . . , ,
P S
,P °°t .,o °~'~
' " , " " = " " " " " ' = =l ' 1 ' b = = ! 48 192 768
I
3072
nel
Fig. 15. Strip with a hole; plane strain; convergence in the energy norm. 2.5 ¢:: O
I,m q[ml
:= ,,¢:
I ~
ellml
I. qlm) r4~
2.0"I
"'''''
PSS/48 PSS/192
•- = - "
Smoothed/768
" " ~
""
emll
r/j r/)
Im r/) alto elm
;I<
.<
,o
"o:4"
o:s" o:s" o'7 "o'8 " o : o ' , o
x21b
Fig. 16. Strip with a hole; plane strain; axial stress distribution along the line A-B.
stress distribution is available. Therefore, the smoothed stresses obtained for the 768-element mesh are used as a reference. The results obtained show the Pian and Sumihara element to exhibit a marginal superiority over the proposed elements for the coarse meshes, only 0.8% for a twelve-element mesh. The formulation shows a small sensitivity of the type of incompatible shape functions used, only 1.3% and 0.04% for three- and twelve-element meshes, respectively. Excellent results are obtained; with only three elements in the mesh, 69% and 68% of the reference solution is obtained for the PSS1,3,5 and PSS2,4,6, respectively. Excellent results are also obtained for the stress distribution inside the elements.
7. Concluding remarks The main objective of this paper is to apply the method proposed by Weissman and Taylor [1] to generate assumed fields for four-node quadrilateral linear elastic plane stress/strain
150
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
elements. O n e of the method's main goals is to avoid locking, which is commonly observed in finite element models of problems involving internal constraints. As an illustration of the treatment of internal constraints, the problem of plane strain at the nearly incompressible limit is chosen due to the difficulty it presents to finite element methods, and the large body of work which has been directed toward eliminating this obstacle. The application of the p r o p o s e d m e t h o d to nonlinear problems will be addressed in subsequent papers. Both elements formulated via the H e l l i n g e r - R e i s s n e r and H u - W a s h i z u variational principle are considered. For both types of formulations, it is proven that as the nearly incompressible limit is approached (plane strain), the trace of the strain vanishes pointwise (i.e., the internal constraint is satisfied pointwise). Consequently, locking at the nearly incompressible limit is avoided. The proposed elements are subjected to an extensive set of problems reported in the literature. Excellent resul~.~; are obtained for all examples. In fact the results are equivalent to those obtained by the best known elements (e.g., [13]). F u r t h e r m o r e , the variationaUy consistent stresses recovered at the element level yield excellent results. Identical results are obtained for elements formulated via the H e l l i n g e r - R e i s s n e r and H u - W a s h i z u functionals, provided the same incompatible functions are used. Also, identical results are obtained for both approaches taken to transform the assumed strain field from the element's natural space into the physical space. F u r t h e r m o r e , results show no sensitivity to the incompatible shape functions used.
References
Ill S.L. Weissman and R.L. Taylor, Treatment of internal constraints by mixed finite element methods: Unification of concepts, Internat. J. Numer. Methods Engrg. 33 (1992) 131-141.
121 T.J.R. Hughes, The Finite Element Method (Prentice Hall, Englewood Cliffs, NJ, 1987). [31 N. Bicanic and E. Hinton, Spurious modes in two-dimensional isoparametric elements, Internat. J. Numer. Methods Engrg. 14 (1979) 1545-1557.
141 T.J.R. Hughes, Equivalence of finite elements for nearly incompressible elasticity, J. Appl. Mech. 44 (1977) 181-183.
151 D.S. Malkus and T.J.R. Hughes, Mixed finite element methods-reduced and selective integration techniques: a unification of concepts, Comput. Methods Appl. Mech. Engrg. 15 (1978) 63-81.
[61 T.J.R. Hughes, Generalizing of selective integration procedures to anisotropic and nonlinear media, Internat. J. Numer. Methods Engrg. 15 (1980) 1413-1418.
171 J.C. Simo, R.L. Taylor and K.S. Pister, Variational and projection methods for the volume constraint in finite deformation elasto-plasticity, Comput. Methods Appl. Mech. Engrg. 51 (1985) 177-208.
181 D.P. Flanagan and T. Belytschko, A uniform strain hexahedron and quadrilateral with orthogonal hourglass control, Internat. J. Numer. Methods Engrg. 17 (1981) 679-706.
191 T. Belytschko and W.E. Bachrach, Efficient implementation of quadrilaterals with high coarse-mesh accuracy, Comput. Methods Appl. Mech. Engrg. 54 (1986) 279-301.
[1o1 E.L. Wilson, R.L. Taylor, W.P. Doherty and J. Ghaboussi, Incompatible displacement modes, in: S.J. [11] [12] [23] [14] [15]
Fenves, N. Perrone, A.R. Robinston and W.C. Schnobrich, eds., Numerical and Computer Models in Structural Mechanics (Academic Press, New York, 1973). R.L. Taylor, P.J. Beresford and E.L. Wilson, A nonconforming element for stress analysis, Internat. J. Numer. Methods Engrg. l0 (1976) 1211-1219. G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice Hail, Englewood Cliffs, NJ, 1973). T.H.H. Pian and K. Sumihara, Rational approach for assumed stress finite elements, Internat. J. Numer. Methods Engrg. 20 (1984) 1685-1695. T.H.H. Pian and C.-C, Wu, A rational approach for choosing stress terms for hybrid finite element formulations, Internat. J. Numer. Methods Engrg. 26 (1988) 2331-2343. C.-C. Wu, M.-G. Huang and T.H.H. Plan, Consistency condition and convergence criteria of incompatible
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
151
elements: General formulation of incompatible functions and its application, Comput. & Structures 27 (1987) 639-644. [161 R. Rubinstein, E.F. Punch and S.N. Atluri, An analysis of, and remedies for, kinematic modes in hybrid-stress finite elements: selection of stable, invariant stress fields, Comput. Methods Appl. Mech. Engrg. 38 (1983) 63-92. [17] E.F. Punch and S.N. Atluri, Development and testing of stable, invariant isoparametric curvilinear 2- and 3-D hybrid stress elements, Comput. Methods Appl. Mech. Engrg. 47 (1984) 331-356. [181 C.-T. Yang and S.N. Atluri, An assumed deviatoric stress-pressure-velocity mixed finite element for unsteady, convective, incompressible viscous flow: Part l--Theoretical Development, Internat. J. Numer. Methods Fluids 3 (1983) 377-398. 119] J.C. Sima a~ld M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Internat. J. Numer. Methods Engrg. 29 (1990) 1595-1638. [2o1 B. Fraeijs de Veubeke, Displacement and equilibrium models in the finite element methods, in: O.C. Zienkiewicz and G.S. Holister, eds., Stress Analysis (Wiley, London, 1965). [211 J.C. Simo, J.G. Kennedy and R.L. Taylor, Complementary mixed finite element formulations of elastoplasticity, Comput. Methods Appi. Mech. Engrg. 74 (1990) 177-206. [221 O.C. Zienkiewicz, S. Qu, R.L. Taylor and S. Nakazawa, The patch test for mixed formulations, lnternat. J. Numer. Methods Engrg. 23 (1986) 1873-1883. [231 I. Babuska, J.T. Oden and J.K. Lee, Mixed-hybrid finite element approximation of second-order elliptic boundary-value problems, Comput. Methods Appl. Mech. Engrg. 11 (1977) 175-206. [24] W.M. Xue and S.N. Atluri, Existence and stability and discrete BB and rank conditions, for general mixed-hybrid finite elements in elasticity, in: R. Spilker et al., eds., Hybrid and Mixed Finite Element Methods, AMD Vol. 73 (ASME, New York, 1985) 377-398. [25l P.G. Bergan and M.K. Nygard, Finite element with increased freedom in choosing shape functions, Internat. J. Numer. Methods Engrg. 20 (1984) 643-663. [26] R.L. Taylor, O.C. Zienkiewicz, J.C. Simo and A.H.C. Chan, The patch test for mixed formulations, Internat. J. Numer. Methods Engrg. 22 (1986) 32-62. [27] R.H. MacNeal and R.L. Harder, A proposed standard set of problems to test finite element accuracy, Finite Elements Anal. Design l (1985)3-20. [281 R.D. Cook, A plane hybrid elemen~ with rotational d.o.f, and adjustable stiffness, Internat. J. Numer. Methods Engrg. 24 (1987) 1499-1508. [29] R.H. MacNeal, A theorem regarding the locking of tapered four-noded membrane elements, lnternat. J. Numer. Methods Engrg. 24 (1987) 1793-1799. 1301 G.N, Savin, Stress concentration around holes, International Series Monographs in Aeronautics and Astronautics (Pergamon, Oxford, 1961). [311 O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, 4th Edition (MacGraw-Hill, London, 1989).
A unified approach to mixed finite element methods: Application to in-plane problems Shmuel L. Weissman and Robert L. Taylor Department of Civil Engineering, University of California at Berkeley, Berkeley, CA, USA Received 31 July 1991 Revised manuscript received 11 February 1991 and 16 October 1991
A general method to treat internal constraints within the context of mixed finite element methods has been presented in previous work. The underlying idea is to constrain the assumed stress and strain fields to satisfy the homogeneous equilibrium equations in a weak sense and thus satisfy a priori the internal constraints. This method is now applied to generate four-node plane stress/strain elements. For these elements, it is proved that locking at the nearly incompressible limit (plane strain) is avoided at the element level. The proposed elements are shown to yield excellent results on a set of standard problems. Furthermore, excellent stresses are obtained at the element level.
1. Introduction
1.1. Motivation and objective Over the last two decades much research has been devoted to the design of plane strain elements that do not lock at the nearly incompressible limit. Numerous methods have been proposed (see below). Unfortunately, most of these methods are either based on parameter tuning or cannot be applied to other problems involving internal constraints (e.g., plate bending). In this work, the method presented by Weissman and Taylor [1] is applied to generate assumed stress and strain fields for four-node plane stress/strain elements. In particular, emphasis is put on proving that the proposed elements do not lock at the nearly incompressible limit. In totality, the proposed elements should exhibit the following properties: - r e s o l v e locking at the nearly incompressible limit (plane strain) at the element level; - r e c o v e r stresses that are variationally consistent with the formulation used to obtain the displacement field; -produce reliable stresses at the element level (without recourse to smoothing); -possess correct rank (no spurious energy modes); - h a v e performance independent of coordinate system or user input data.
1.2. Background Solutions of plane strain problems at the nearly incompressible limit have traditionally posed severe problems in the development of finite element models. This setback is the result of the failure of finite element models to impose the incompressibility constraint pointwise
(e.g., [21).
Much effort has been directed at producing finite element models that avoid the locking phenomenon at the nearly incompressible limit. Initial approaches utilized reduced integra0045-7825/92/$05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved
128
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
tion. However, this approach often results in unstable elements due to the appearance of spurious zero energy modes [3]. To overcome this difficulty, Hughes [4] used selective reduced integration (SRI). Malkus and Hughes [5] showed the SRI scheme to fall within the concept of mixed finite element methods. Hughes [6] refined the SRI scheme into a general method, known as the B-bar method. Simo et al. [7] showed that the B-bar method can be derived from the Hu-Washizu variational principle. Using a projection operator, Flanagan and Belytschko [8] showed that the usual finite element approximation can be rewritten in a form that leads to decoupling of the stiffness matrix, which can be written as the sum of the decoupled under-integrated stiffness matrix and a 'stabilization' stiffness matrix. Following this approach, Belytschko and Bachrach [9] presented formulations based upon the Hu-Washizu functional. The best of the elements presented is known as the quintessential bending incompressible (QBI). This element, however, is not frame invariant (i.e., it is dependent upon the coordinate system) unless a unique local coordinate system is defined and the element stiffness matrix and load vector, defined in the local coordinate system, are transformed to the global coordinate system. Wilson et al. [10] took yet another approach They observed that when the four-node element (isoparametric displacement formulation) is subjected to pure bending, it deforms in shear rather than in bending. As a result of this anomaly, the element is stiff in bending. In order to improve the element bending behavior, they introduced a set of incompatible displacements. It was soon realized, however, that the distorted element did not pass the constant strain patch test. Taylor et al. [11] modified the incompatible modes element so that the resulting element passed the patch test. In addition to the improved behavior in bending this element avoids locking at the nearly incompressible limit (this property was not observed, however, when the elements were first presented). Conditions for convergence when incompatible modes are used were presented by Strang and Fix [12]. Pian and Sumihara [13] used incompatible displacements to generate an assumed stress field for a four-node plane stress element, which was formulated via the Hellinger-Reissner variational principle. This element exhibits excellent characteristics in bending applications. In addition, its sensitivity to mesh distortion from a parallelogram shape appears to be among the smallest of four-node elements. Furthermore, when modified to account for plane strain constitutive relations, it avoids locking at the nearly incompressible limit. Initially, the formulation required a quadratic perturbation of the element shape. Recently, Pian and Wu [14] presented a new formulation that avoids this perturbation. A general formulation to generate incompatible element functions was presented by Wu et al. [15]. Punch et al. [16] presented a theoretical framework for the selection of assumed stress fields in the context of elements based upon the Hellinger-Reissner functional. Punch and Atluri [17] presented a family of two- and three-dimensional elements. The four-node element yielded comparable results to the Pian and Sumihara [13] element. The performance of these elements when applied to model nearly incompressible materials was demonstrated by Yang and Atluri [18] for the case of the Navier-Stokes equations. Simo and Rifai [19] presented four-node elements based upon the Hu-Washizu functional which introduced 'enhanced strains'. Results reported for these elements are almost identical to the results obtained by the element presented by Pian and Sumihara [13]. To avoid the limitation principle [20], the assumed stress and strain fields must not contain as a subset the corresponding fields derived from the displacement field. In addition, to avoid locking behavior in the presence of internal constraints, the assumed stress and strain fields must a priori satisfy these constraints pointwise. Weissman and Taylor [1] proposed constructing these assumed fields to possess the aforementioned properties by constraining, in-
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
129
dependently, the assumed stress and strain fields to a priori satisfy the homogeneous equilibrium equations in a weak sense.
1.3. Paper overview The stress and strain field generation procedure is reviewed in Section 2. The strong form for the plane stress and plane strain problem is summarized in Section 3. The finite element approximation is contained in Section 4. In Section 5, four-node quadrilateral elements formulated via either the Hellinger-Reissner or the Hu-Washizu functionals are proposed for both the plane stress and plane strain problems. It is proved that in the case of plane strain the trace of the assumed strain vanishes, pointwise, as the nearly incompressible limit is approached. Consequently, locking is avoided. Numerical examples demonstrating the excellent properties of the proposed elements are presented in Section 6. Concluding remarks are given in Section 7. Before continuing to Section 2, there is an important point to be made. The equivalence between the formulation based on the Hu-Washizu functional and that based on the Heilinger-Reissner functional, as evident from the numerical examples in Section 6, holds only in the case of linear elasticity. When the constitutive equations are nonlinear, the derivations of a consistent tangent stiffness rely crucially upon the use of the Hu-Washizu functional [7], while the Hellinger-Reissner functional leads to restrictive (cannot model perfect plasticity) and cumbersome algorithms [21]. Consequently, it is the Hu-Washizu functional that is of interest if the methodology presented is subsequently to be used in nonlinear problems.
2. Summary of assumed fields generation procedure Internal constraints pose a difficult problem for finite element methods, commonly manifested by locking. This difficulty was put in the form of a limitation principle by Fraeijs de Veubeke [20]. In order to avoid this limitation principle, the assumed stress and strain fields must satisfy a priori the equilibrium equations, and thus effectively satisfy the internal constraints. In the context of finite element methods this requirement can be relaxed, and the equilibrium equations can be satisfied in a weak sense. Taking advantage of this observation, Weissman and Taylor [1] presented a general method to generate assumed stress and strain fields. These fields are constructed to satisfy in a weak sense the homogeneous equilibrium equations, and consequently the internal constraints. This method is formulated within the framework provided by the Hu-Washizu variational principle, and is summarized below. Before proceeding to describe the proposed method it is useful to recall the steps involved in deriving finite elements via the Hu-Washizu functional. To this end let e, o" and U denote the assumed strain, stress and displacement fields, respectively; L the strain displacement operator; and D the elastic coefficients matrix. The Hu-Washizu functional, for the case of small deformations and linear elastic materials (ignoring initial stresses and strains) is given by HH(Cr, e,
U):= fa [½etDe+ (rt(LU- e)]da-//~xT(U),
(2.1)
where g2 is the domain, and //EXT is the external work. The weak, or variational, form may be obtained from the energy functional,//H, by making
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
130
it stationary. This result is obtained by taking the first total variation of H H and equating it to zero. To this end, let ~e, ~o" and 5U denote virtual strain, stress and displacement fields, respectively. The weak form is given by
DII n := ~
[~et(oe
-
o') +
~o't(LU - e) + (L ~U)'o" - ~Utb] d O -
fc
5U'o'" d F = O,
t
(2.2) where b is the body force vector, o-" are applied traction boundary conditions, and C, is the part of the boundary on which traction boundary conditions are specified. The finite element model is obtained by approximating the assumed fields. To this end, let the assumed fields be approximated over each element domain by - Stress field: o- = Ss, - Strain field: e = Ee, - Displacement field: U = Nd, where S, E and N are the shape functions for the approximating stress, strain and displacement fields, respectively; and s, e and d are the vectors of independent stress, strain, and displacement parameters, respectively. Substituting the assumed fields into the weak form leads to the following system of equations (these are the approximated Euler-Lagrange equations):
-A
0
0
G t
where
= 0
(2.3)
,
d
A " = f, StE d O ,
n "= f, EtDE dO , f G ' = I StB dO
.It
f
"= fg~ Ntb dO + fc 't Nttra d F
(2.4a)
'
(2.4b)
and B is defined by B : = LN. Eliminating the ~tress and strain coefficients yields
Kd=f ,
(2.5)
where K is the element stiffness matrix, given by X=
.
(2.6)
To satisfy stability conditions for every admissible displacement, the following inequalities must hold [22]: n~+n d~n,,
n,~>n d,
(2.7)
where n~, n~ and n~ are the number of strain, stress and displacement parameters, respectively (n d is equal to the number of nodal degrees-of-freedom). ~ ~The stability conditions presented in (2.7) are necessary but not sufficient. In addition some constraints must be placed on the admissible stress and strain shape functions S and E, respectively. In particular these functions should be chosen so that the LBB condition is satisfied (see e.g. [23, 24]).
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
131
A significant reduction in the computational effort can be obtained if A is invertible. Henceforth, A is assumed to be invertible. 2 As a result n~ = n~ and, consequently, the mixed patch test [22] provides the minimal number of independent stress and strain parameters. Furthermore, in view of the desire to minimize computations, this is the optimal number. As was stated above, the underlying idea behind the proposed method is to construct assumed stress and strain fields so that they a priori satisfy the homogeneous equilibrium equations in a weak sense. To this end, let the assumed strain field be given as the sum of two independent fields as follows: e := g': + e i
(2.8)
and let, for convenience, e ~ be given by i
e "=LU
i
(2.9)
where U i is a displacement field not contained in U h, which is the finite element approximation of the admissible displacement solution space. The initially assumed stress, if, and compatible strain, F c, fields are constrained so that the contribution to the energy functional, H H, from the mixed terms containing e i and either U or ff vanishes in a weak sense over each element's domain. As a result, the reduced stress, or, and compatible strain, e~ satisfy the homogeneous equilibrium equations in a weak sense. Furthermore, it follows from the variational structure that the incompatible strain, e ~, is zero pointwise (at the solution). The proposed method is summarized in Box 1. The condition U ~ ¢ U h (in Step 2) is introduced in order to prevent the appearance of spurious zero energy nodes (see [1] for discussion). The constraint 3(a) was introduced by Strang and Fix [12] as a requirement to pass the constant stress/strain patch test. Constraints 3(b) and 3(c) arise from the variational statement of the problem. Requirements similar to 3(b) were used by Pian and co-workers (e.g., [13, 14]). However, while the objective was the same (i.e., to constrain the assumed stress field to a priori satisfy the homogeneous equilibrium equations in a weak sense), the constraint equation itself was different, and motivated by different arguments. Constraint 3(c) was also presented in the context of the free Box 1. Proposed method 1. Select initial assumed (enriched) strain field, U, and stress field, (Pascal triangle). 2. Select the incompatible strain field, e i ( U i ~ U h). 3. Reduce the number of independent variables in U and ~ by introducing the following constraints: (a)
f ~ U .,i d ~ = 0 ,
(b) (c)
fl ei'D~ c dg~ = 0
--~ e" .
4. Substitute the reduced fields into the variational principle.
2
This assumption crucially depends on the selection of S and E ,zo that the LBB condition is satisfied.
S.L. Weissman, R.L. Taylor, A unified approach to m&ed FEMs
132
formulation (e.g. [25]). However, this condition was introduced by Bergan et al. not as a constraint on the assumed compatible strain field, e c, but rather as a constraint on the admissible higher order modes (i.e., incompatible displacements).
REMARK. Recall
that the Hellinger-Reissner functional can be obtained from the H u Washizu functional by enforcing the stress/strain relationship (i.e., tr = pointwise. If this relationship is substituted directly into the Hu-Washizu functional, e i drops out and, consequently, the proposed method cannot be used. However, if the formulation is started with a Hu-Washizu functional and the constraints enforced there, while taking notice of the special choice for the strain field, then the same reduced stress field is obtained and e i drops out only after obtaining the desired stress field.
De)
3. Strong form The strong, or classical form for both the plane stress and plane Strain pa-oblems is presented. Throughout this section Greek subscripts take the values 1 and 2. Repeated indices imply the usual summation convention. All quantities are referred to a fixed system of rectangular, Cartesian coordinates. A general point in this system is denoted by (x~, x 2, %). The domain of interest is a three-dimensional body embedded in an Euclidian three-space R 3. In the undeformed configuration, the domain O is of the following form:
{
I
J'~ "-" (X I , X2, X3) {~ R 3 X3 E
[-h(xn'x2) 2} 2 ' +h(xI'x2)],(Xn,X2)~ACR 2
where h is the thickness in the x 3 direction ~ and A is a closed region bounded by a simple closed curve C with an outgoing unit normal v. Let C U and C, be subregions of C such that C U U C, = C and C u N C, = t~, where C u is the part of C on which displacements are specified and C, is the part of C on which tractions are specified. The strong form is summarized in Box 2 where b,, are the body forces, o-~ are the specified traction boundary conditions, U~ are the specified displacement boundary conditions and C , , ~ is the elastic coefficients tensor, which for the case of plane strain is given by C,~
:= hI/.~(3,~t$~ + ~ )
+ aS,~tSvsl,
(3.1)
where a and /z are the Lam6 coefficients, which can be eliminated in favor of Young's Box 2. Plane stress/strain: strong form Given b , , o'i', and U~',; find U,, and tr~ such that
in A
Io',,~,t~ + b = 0 , { o',~ = C,,~,t e,~ ,
l~.,~-- ~(u,,,~ + on ct, on
C,
U~.,,);
u,, = u"" a
vt3tr,v = or, .
3 In the case of plane strain h - 1.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
133
modulus, E, and Poisson's ratio, v, using vE A= (1 + v ) ( 1 - 2v) '
E /z = 2(1 + v) "
(3.2)
The plane stress elastic coefficients are obtained by replacing A by J(, where ;( is given by 2/xA = vE ' ( = A +2/x ------5 1-v "
(3.3)
Note that in the case of plane stress, the strain terms in the x 3 direction can be obtained from the constitutive equations. Similarly, in the case of plane strain, the stress terms in the x 3 directions can also be obtained from the constitutive equations. Consequently, the in-plane stress and strain fields characterize the three-dimensional state for both problems. It follows from Box 2 that the strain displacement operator L is given by "
b
oxz L:=
0
0
(3.4)
Ox 2
0 b Ox2 Oxl-
s
0]
and the elastic coefficient matrix D is given by
D= 1-V
2
(3.5a)
v 1 0 0 0 ½(1 -- V)
for the case of plane stress, and by
E
D = (1 + v ) ( 1 - 2 v )
I, 0v vv1 -0 v
° 1
0 ~(1-2v)
(3.5b)
for the case of plane strain.
4. Finite element approximation The method summarized in Section 2 is used to formulate four-node quadrilateral plane stress and plane strain elements. The initially assumed stress and strain fields as well as the incompatible displacements used to generate the assumed incompatible strains are presented. In Section 4.1, the assumed displacement field is introduced. The assumed stress field is given in Section 4.2, and the assumed compatible strains in Section 4.3. The incompatible displacements used to generate the assumed incompatible strains are summarized in Section 4.4. 4.1. Assumed displacement field
The assumed displacement field is taken as a standard isoparametric field. The displacements are approximated by
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
134
(4.1)
U = Nld I ,
where N t is the shape function associated with node I, and d I is the nodal displacement vector• In the case of four-node quadrilateral elements, the shape functions are given by (4.2)
Nt( ~, 7) = -~(1 + ~:t~)(1 + n, 7/),
where ~: and ~ are the element natural coordinates, on the interval [ - 1 , 1], and ~1 and r h are the values of the natural coordinates at node I. 4.2. A s s u m e d stress fields
The assumed in-plane stress field is a complete linear field 4 expressed in the element natural coordinates as
t~*-
~~ tT,~ -
1
~ 0 0 0 0 0
I *l 0' 0, [s,;j 1
0
o 0 o
0
0 lrs~ 0 s*
l! ]
(4.3)
e " LZ
Since complete polynomials are used to express the in-plane stresses, (4.3) could be used for ~ directly. However, the reduction to satisfy the constraints introduced in Box 1 would require selection of different parameters in each element (i.e., there would be a dependence on element orientation of ~ and 77with respect to x~ and x2). This dependence may be avoided by using the transformation procedure described below. The following definitions are introduced: xs= ~txtt,
x t = ~rltxt,,
xh=
¼(serl)txt,,
ys = ¼~,xai,
y t = ¼rhx2,,
y h = ~( grl)tx2t.
Following Zienkiewicz and Taylor [31], the Jacobian of the coordinate transformation from the (~, 7/) space to the (x~, x,) space is given by
J = Jo +
+ J2n,
where Jo = xs y t -
xt y s ,
(4.4) (4.5a)
Ji = xs y h - x h ys ,
(4.5b)
J2 = x h y t -
(4.5c)
xt yh .
The stress field in the physical space is obtained by using the following transformation:
1 crq = ~ Fi,FkKtr lr ,
(4.6)
4 The stress resultant and strain fields are assumed as complete linear fields since this is the best assumption that may be used in conjunction with a bilinear displacement field.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
135
where both i and j take the values x, or x 2, and both I and J take the values ~ or v, and Fxl ~ = Ox~/O~, etc. However, to maintain the ability to model constant fields (i.e., pass the constant strain patch test), the transformation is performed based on the Jacobian, J, and 'deformation gradient', F, at the center of the element. At the center of the element, F is given by Fx, ~ = xs ,
Fx,, = xt ,
(4.7a)
Fx._e = ys ,
Fx., = y t .
(4.7b)
After redefining the independent coefficients, the assumed stress field is given by
=
I°112[!oo [ =
tr~2J
1
0
0
1 xs ysrt
xt2, x2,
yt2rt
2 ys yt ~
2 ys yt71
xt yt ~ xs ys ~
xt ytrt
A(
Art
yt z ~
ys2rt
2
ys z ~
sx,
2xxt ][iIs2 9
(4.8)
where A = xs yt + xt y s , 1
S I ~'- ~
1
S2 = ~
,
,
( y s 2 S~ + yl z S 4 + 2 ys y t s 7 ) ,
1
s3=~(xsyss~+xtyts 1
Si=~Sj
,
(XS 2 S 1 + x t 2 S 4 "~- 2 X S X t S 7 ) ,
,
,
4 +As 7),
,
,
with (i, j ) E {(4, 3), (5, 5), (6, 2), (7, 6), (8, 8), (9, 9)}. It is convenient to write or in the following form:
where S,,t2 ._ ( $ 6 , S7, $8, $9 ) .
With the above construction, the parameter set g2 may always be selected as the set to be eliminated in satisfying the constraint equation on the stress field (see Box 1). 4.3. A s s u m e d strain fields
The assumed strain field is formulated in the element natural coordinates, and then transformed into the physical domain. Two types of transformations can be used: - t h e same transformation as used for the corresponding stress fields; - t h e inverse to the transformation used for the corresponding stress fields. The motivation for the first approach is to have the same 'shape functions', or interpolation, for the strain fields as for the corresponding stress fields. The second approach is motivated by
s.L. Weissman.R.L. Taylor. A unified approach to mixed FEMs
136
the invariance of the complementary energy (trijeij) under coordinate transformation. Accordingly, the transformation is given by 5
eo = JF~' F~' e,,*.
(4.9)
A similar approach was taken by Simo and Rifai [19]. Following the path established for the assumed stress field, the strain field is assumed as a complete linear field in the element natural coordinates, and then transformed into the physical space• Using the transformation given by (4.6) yields e4 , El2
[!
xsErl ys2~
0 0
=
1 0
0 1 xsysrl
xt2~
xs2~ ys 2~
xtErl yt2~
2xsxt~ 2ysyt~
xtyt~ xsys~
xtyt,1
A~
yt 2~:
2xsxt'o][i']e2 2ysyt~ (4.10a) Arl
and using the inverse transformation, (4.9), yields
i ~'C
E
0 0 yt 2 71 1 0 xtErl 0 1 -2xtyt77
=
X
ys2 + xs z ~ -2xsys~
yt 2 ~ xt 2 t~ -2xtyt~
ys2 ,1 xs2rl -2xsys~l
Ill e., I-
-ys yt + -ys yt,ll -xsxt ~ -xsxtrl I" A~ A~q J (4.10b)
•
tj
It is convenient to write e in the following form:
[A]
+: = [a, I a..,l e,
where
-
e2
'
At
el = (el, e2, e3, e4, es) ,
~ = ( e 6 , e7, e~, e g ) .
With the above construction, the parameter set ~2 may always be selected as the set to be eliminated in satisfying the constraint equation on the strain field (see Box 1).
4.4• Assumed incompatible displacements Let the assumed incompatible displacements be given by •
"
i
Ult = N'tAt + NEA2
(4.11a)
Ui2 = NilA3 + Ni2A4 ,
(4.11b)
and where At, '~2,. • •, and Aa are the independent incompatible displacement parameters, and Ni,~ are the incompatible shape functions. Two options are selected for the incompatible shape functions Nit and Ni2. The first set was The transformation is based on the values at the center of the element.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
137
"presented by Wu et al. [15]. The functions are given by •
N'~ = se2 and
2J2 3,/o ~: + ~oo r/ 2J l
(4.12a)
2J I 2J2 N~ = "r/2+ ~ ~ - 3J,--~7/.
(4.12b)
The second set was presented by Taylor et al. [26]. In order to obtain a more compact form, the incompatible modes given below are a linear combination of the modes originally presented. The functions are given by N , - ( 1 - ~ J2 n)(1 - ~:2) + ~J,: ( 1 - T I
2)
(4.13a)
and N'2 = ( 1 - J' ~ : ) ( 1 - 2 ) + J , f~ 71(1- ~:2).
(4.13b)
Note that the second set of incompatible modes is zero (compatible) at the nodal points while the first set is not.
5. Proposed elements
The fields presented in Section 4 are used to generate four-node plane stress/strain elements. The proposed elements are cataloged in Table 1 by the type of formulation, the incompatible shape functions, and assumed strain field (Hu-Washizu). For all these elements, the assumed displacement and stress fields are given by (5.1) and (4.8), respectively. The remainder of this section consists of proofs (for the case of plane strain) showing that as the incompressible limit is approached, the assumed strain vanishes pointwise for the proposed elements. Consequently, the well-known locking behavior is avoided at the element level. The proof for elements formulated via the Hellinger-Reissner and Hu-Washizu variational principles are given in Propositions 5.1 and 5.2, respectively. P R O P O S I T I O N 5.1. The Hellinger-Reissner elements presented above can, in the case o f plane strain, model a state in which the trace of the strain goes to zero, pointwise, as v goes to
0.5-.
Table 1 Proposed plane stress/strain elements Element
FormuNtion
Incompatible shape function equations
PSS1 PSS2 PSS3 PSS4 PSS5 PSS6
Hellinger-Reissner Hellinger-Reissner Hu-Washizu Hu-Washizu Hu-Washizu Hu-Washizu
(4.12) (4.13) (4.12) (4.13) (4.12) (4.13)
Strain field representation -(4. l 0a) (4.10a) (4.10b) (4.10b)
138
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
PROOF. Recall that in the case of the Hellinger-Reissner variational principle the assumed strain field is obtained from the assumed independent stress field. Consequently, in the case of plane strain the strain field e is given by
e
=
D-lo"
(5.1)
where D -l is given by
o
-l = _
A
a
o 0] 2/z
0
(5.2)
,
allx
with a = 4/z(/~ + A). As Poisson's ratio, v, goes to 0.5-, A goes to infinity while/z remains bounded. It follows that, at the limit, the following relation is obtained: --
D 111=
--|
D22
--
- D l,,l = -D2t I .
-"
(5.3)
After imposing the constraint equations introduced in Section 3, the assumed stress field presented in (4.8) is reduced to the general form
0 a2r/+b2G 1 a3'r/ + bag
1
LO'12J
0
asG+b~,H/.~/, a6G +
b6nJLssj
(5.4)
where a t, a 2 , . . . , and at,, and bl, b 2 , . . . , and b 6 are constants to be determined by the procedure presented in Box 1. Substituting (5.2) and (5.4) into (5.1) and taking notice of (5.3) yields the desired result. [] P R O P O S I T I O N 5.2. The Hu-Washizu elements presented above can, in the case o f plane strain, model a state in which the trace of the strain field goes to zero, pointwise, as the incompressible limit is approached. PROOF. After constraining the assumed strain field introduced in (4.10), the reduced field is given in a general form by
Jell]
e =/ezz/= Lel2]
al rl + blG a 2rl + b2G asG + b 5 1 a 3rl + bag a6G +
[1
0
0
[0
1
0
0 0
b6rt.lL/~s
(5.5)
The structure imposed by the constraint equations is such that the following relations are obtained, as v goes to 0.5-, between the at, a 2 , . . . , and a 6 and b l, b 2 , . . . , b 6 coefficients: a~ + a 2 = 0 ,
a4 + a5 = 0
(5.6a)
b I + b2 =0,
b4 + b5 =0.
(5.6b)
and
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
139
A formal proof for (5.6) can be obtained. This proof, however, is very cumbersome and tedious, and does not yield a better understanding of the formulation. For these reasons, a numerical 'proof' will be given. To obtain a more compact form, an equivalent form to (5.6) is given. Namely, the values of (a~ + aE)r/+ (b I + bE)g and (a 4 + as)~ + (b4 + bs)r/are reported at the Gauss points. Since the functions are linear, zero value at all Gauss points implies zero value pointwise. First, a square element (see Fig. l(a)) is examined. The values for the first Gauss point are reported in Table 2. Values for other Gauss points may be obtained as follows: given that the value at the first point is o, then at the other points it is -sign(C) x v, where sign(~) is the sign of the ~: or r/coordinate at the given Gauss point (7/for the first column and ~ for the second column). Secondly, an element distorted from a parallelogram shape (see Fig. l(b)) is examined. In this case the values at the first and second Gauss points are reported. The value for the third Gauss point is minus the value at the second Gauss point, and the value at the fourth Gauss point is minus the value at the second Gauss point. Results are summarized in Tables 3 and 4. The results summarized in Tables 1, 2 and 3 'prove' (5.6). It follows from (4.12) and (4.13) that as v goes to 0.5- the trace of e is given by
(5.7)
trace(e) = e I + e : .
In order to show that the trace of the strain goes to zero pointwise recall from Section 2 that the stress strair~ weak relation is given by He
-
As
X2
(5.8)
= O .
X2
11
' (o.1)
(0.25,1.25)
(1,1)
(1.2S,0.5)
(0,0)
(~ ,o) X1
X1
('o) Fig. 1. (a) Square element; (b) element distorted from a parallelogram shape.
Table 2 Square element
0 0.3 0.49 0.499 0.4999
(a I + a2)"O + (b I + b2)~
(a 4 + as)~ + (b4 + b5)7/
-0.144338 -0.082479 -0.005661 -0.000577 -0.000058
-0.144338 -0.082479 -0.005661 -0.000577 -0.000058
140
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
Table 3 Distorted element (a, + a2)'q + (b~ + b2)~
v
Point
PSS3
PSS4
PSS5
PSS6
0
1 2
-0.138299 -0.142666
-0.136318 -0.145547
-0.239751 -0.234665
-0.166003 -0.224427
0.3
1 2
-0.074126 -0.077283
-0.072729 -0.079388
-0.2255061 -0.204323
-0.231639 -0.169421
0.49
1 2
-0.004723 -0.004982
-0.004612 -0.005156
-0.051700 -0.055967
-0.037454 -0.042924
0.499
1 2
-0.000479 -0.000505
-0.000467 -0.000523
-0.006410 -0.006890
-0.004384 -0.004980
0.4999
1 2
-0.000048 -0.000051
-0.00U047 -0.000052
-0.000650 -0.000700
-0.000446 -0.000506
Table 4 Distorted element (a 4 + a~)~j + (b4 + bs)~ v
Point
PSS3
PSS4
PSS5
PSS6
0
1 2
-0.099318 +0.102311
-0.096704 +0.103024
-0.386263 +0.388999
-0.414609 +0.420426
0.3
1 2
-0.053029 +0.055193
-0.051169 +0.055728
-0.692734 +0.703892
-1.013897 +1.045209
(I,49
1 2
-0.003365 +0,003542
-0.003215 +0.003588
+0.027529 -0.029822
+{I.021353 -0.024105
0.499
i 2
-0.000341 +0.000359
-0.000326 +0.000364
+0.002510 -0.002771
+0.001978 -0.002279
0.4999
1 2
-0.000034 +0.000036
-0.000033 +0.000036
+0.000249 -0.000276
+0.000197 -0.000227
To simplify notations and without loss of generality, ~ and 7/in the a s s u m e d strain field, eq. (5.5), may be replaced by f and ~, given by
Go,
h = n - rto,
w h e r e ~0 and r/0 are shifts introduced in order to m a k e H " block diagonal. In order to o b t a i n this result the following r e q u i r e m e n t s must be met: f,,t ~ d A = O '
fA'~dA=O"
T h e resulting Go and 7/o are given by
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
J, 3Jo'
141
7/° = 3J0"
The H matrix is now of the following form: H=
[o0a /~o] ,
(5.9)
where D is given by (3.5), A is the element area, a n d / t is the 'material' properties for the linear part of the strain field. Based on results in (5.6), /~i (~'= {e 4, es} ) goes to zero pointwise as the incompressible limit is approached. Hence, it remains only to consider the constant part, given by D A . Combining (5.8) and (5.9) yields 1 e, = a A
and
1 e2 = a A
{(A + 2/x)§,- Ag2}
(5.10a)
{(A + 2/~)g2 - Ag,},
(5.10b)
where g = A s . Substituting (5.10) into (5.7) yields trace(e) =
2/x(~, + s2) aA
(5.11)
"
As v goes to 0.5-, a goes to infinity while/~ is bounded; hence the desired result is obtained.
6. Numerical examples
The performance of the plane stress/strain (PSS) elements proposed in Section 5 is evaluated with several discriminating problems selected from the literature. The purpose of these evaluations is to test the proposed formulation's sensitivity to the specific choice of the incompatible shape functions as well as the overall performance of the proposed elements. First tackled are the constant strain patch tests in Section 6.1. Following these are the beam problems suggested by MacNeal and Harder [27] in Section 6.2, sensitivity to mesh distortion in Section 6.3, the Cook membrane [28] in Section 6.4, the thick walled cylinder also introduced by MacNeal and Harder [27] in Section 6.5, and a strip of finite width with a hole in Section 6.6. Convergence of the results obtained by the four-node elements presented in this paper are compared with other well-known four-node plane stress/strain elements. A listing of these elements, and the abbreviations used to identify them is given in Table 5. In all the examples given in this paper, identical results are obtained for elements Table 5 Reference elements Element
Description
2x2 P-S QBI
Isoparametric plane stress/strain Pian and Sumihara [13], plane stress/strain Belytschko and Bachrach [9], plane stress/strain
142
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
formulated via the Hellinger-Reissner variational principle and the corresponding (i.e., same incompatible shape functions) elements formulated via the Hu-Washizu variational principle. Identical results are also obtained for the two approaches taken to transform the assumed strain field from the element's natural space into the physical space. Convergence is examined in the energy norm (energy norm reported is twice the internal strain energy) since this is the natural convergence test in the context of finite element methods [12]. Convergence of displacement at some characteristic points is also examined as this is the common practice in the literature. All tables and figures show the displacement/ energy norm as a function of the number of elements (denoted nel) used in the corresponding mesh. 6.1. Patch test
A rectangular domain is modeled by a single element shown in Fig. 2(a) and the skewed mesh shown in Fig. 2(b). The mesh is subjected to a constant state of tension/compression. All elements presented pass this test. 6.2. Beam problems: in-plane shear
These problems were suggested by MacNeal and Harder [27] as standard problems to evaluate the performance of different elements. The meshes used contain only one row of six elements for both the straight beams, shown in Fig. 3, and the curved beam, shown in Fig. 4. Geometrical and material properties used are summarized in Table 6. The 'exact' tip displacement, w, given by beam theory [27] for the case of a straight beam is w = 0.1081 and for the case of a curved beam is w = 0.08734. The results, normalized with respect to the exact solution, are summarized in Table 7. Results obtained for all proposed elements are identical to the results obtained by the Pian
x2
x
(0,
,12)
n~
2 10.00,0.12)
(0.24,0.12)
~0.08.0.08)
(0.24,0.00)~ ' x
(0.00,0.00) (a)
•~0.16.0.0B)I
" (0.~,0.00)
I, ~'X (0.24,0.00)
(b) Fig. 2. Patch test: (a) one-element mesh; (b) skewed mesh.
Table 6 Material and geometrical properties for beam problems
Straight b e a m Curved b e a m
E 1.0E+7 1.0E+7
v 0.3 0.25
Thickness 0.1 0.1
Length/ inner radius 6.0 4.12
Width / outer radius 0.2 4.32
Arc m 90°
143
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
i /
i \
J
L,,o%
(z)
(b)
Fig. 3. Straight cantilever beam. (a) Regular shape elements; (b) trapezoid shape elements; (c) parallelogram shape elements.
Table 7 In-plane shear
LqDqLWJW.
Mesh
PSS/P-S
2x2
QBI
Straight (a) Straight (b) Straight (c)
0.9929 0.2208 0.7963
0.0933 0.0278 0.0348
0.9929 0.0777 0.0890
Curved
0.9782
0.0740
0.2324
ql, W . ' 4
Fig. 4. Curved beam mesh.
and Sumihara [13] element. The 2 x 2 element locks for all meshes. However, if two rows of elements are used, the 2 x 2 element yields reasonable results. All elements exhibit locking when trapezoid-shape elements are used (straight mesh (b)). This is in accordance with the theorem put forth by MacNeal [29] regarding the locking of tapered four-node membrane elements. The QBI element appears to be much more sensitive to mesh distortion than the PSS or P-S elements.
6.3. Beam bending: sensitivity to mesh distortion In th~s standard test, a beam is modeled by two elements. The beam is fixed at one end and subjected to a bending moment on the other, as shown in Fig. 5. The material properties used
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
144
10.0
:::::::::::::::::::::::::
:::::::::::::::::::::::::::
-F
_.__.._.~
Fig. 5. Beam bending problem, sensitivity to mesh distortion.
are E = 1.0 and p = 0.4999. A state of plane strain is assumed. The edge separating the two elements is gradually rotated, as shown in the figure, to skew the mesh. Results, normalized with the exact beam theory solution (w - 562.575), are shown in Fig. 6. Results for the 2 x 2 element are not reported since it exhibits severe locking. All proposed elements yield identical results to the Pian and Sumihara [13] element. The QBI element exhibits rapid deterioration as A is increased.
6.4. Cook's membrane problem The problem consists of a trapezoidal plate clamped on one end and loaded by a uniformly distributed in-plane bending load on the other end, as shown in Fig. 7. This problem has a considerable amount of shear deformation. It was suggested by Cook [28] as an excellent problem to test membrane elements using skewed meshes. No analytical solution is available for this problem. The best known solution is given by Cook as V ¢ - 0.29 (vertical displacement), O"A 0.236 (maximum principal stress) and o-n = - 0 . 2 0 1 (minimum principal stress). The material properties used are E = 1.0, v = 1.0/3.0, thickness = 1.0; plane stress assumption was used. The results normalized with respect to these best known answers are shown in Figs. 8 (Vc), 9 (trA) and 10 (o'~3). The results presented for one element are obtained by interpolating the nodal values; for the finer meshes the nodal values are reported. The results obtained for all proposed elements =
1.5
m , . ,m
PSS/P-S
=="='=
~1 EXACT
E
1.0
o ¢1 .m m O,, r~
N Q=m
o o
¢1 ¢1
E t.. O
0.5
0
i,=
0
4
QOOOQQ~IIUQ
B U N ii iI U iN II ii Ii D a B U l l
0.0
•
0
;
j QB I n i m I B u M
|
3
|
,
s
A
Fig. 6. Beam bending, sensitivity to mesh distortion: normalized tip displacement.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
c
u
m b
145
ti.o 16.0 row,
44.0
44.0 d
:~.:.::-" .,... I[:i:~:i:!?i:~:i...!:~:
ti!~iii~!i!~ii:i i ... i.. .
q
...
.... I[:i!i:!.i:i:!::i}:::ii!i. ..... lt!!~!i :i~!i!:!::i:i
48.0
I[ !! ili: -
Fig. 7. Cook's membrane problem. ¢¢
1.25
cU
E CU c~
1.00
m
D. r~
em
"O
S
0.75
m
..o~
¢1
em
I_ CU
o° #
0.50 ,
o°~ ~ o°°;#
cU N
emil m
,-~.m
PSS/P-S
~mm
2X2
0.25
EXACT
M
E I,, 0.00
1'
;
,'s
6'4
nel
Fig. 8. Cook's membrane problem; vertical disp'~acement at point C.
yield identical results to the Pian and Sumihara [13] element. Excellent results are obtained even for very coarse meshes; about 70% of the vertical displacement and 80% of the stresses are obtained using only one element. The 2 x 2 and QBI elements yield poor results for coarse meshes; less than 50% and 60% of the vertical displacement and stresses are obtained by the 2 x 2 and QBI elements, respectively, when a four-element mesh is used.
6.5. Thick-walled cylinder This problem was suggested by MacNeal and Harder [27] as an excellent problem to test the effect of nearly incompressible material. The mesh used is shown in Fig. 11. The material
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
146
1.25
E E
1.00
e.m
M m
E
f..
o.zs
"Or~ . ag . ~"
•
0.50
-
- - .
@mm
""
"-- "
2X2
IN IN II ID II IN,
EXACT 0.00
I
I
1
4
l's
I
s'4
25s
nel
Fig. 9. Cook's membrane problem; maximum principal stress at point A. 1.25
E E
1.oo,
.,.~1;;~ "~-
•E E ~~
.,*
0.7S'
L..
~o
•, .
oO'~;"
S 4p
. m.
J
m
.-
•E
~
0,25
,,, "-
~ ~ ~ ~
~
,
0.00
-
PSSIP-S 2X2
"''''"
(~1
~
EXACT
,
,
,
le
64
256
nel
Fig. 10. Cook's membrane problem; minimum principal stress at point B.
.0 °
P = I.@ 3:0
3.$
,I.2
S.2
6.7S
9.0
Fig, 11. Thick-wailed cylinder problem.
properties are E = 1.0 and u varies from 0.0 to 0.4999. The exact radial solution is given by
u = E ( R ~ - R2,) T where p is the internal pressure; R~ is the inner radius; R• is the outer radius; and r is the radius where the radial displacement, u, is to be computed.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
147
Table 8 Thick-walled cylinder v
PSS/P-S
2x 2
OB!
0 0.3 O.49 O.499 O.4999
11.9987 0.'.~9~4 0.9913 0.9910 0.9910
11.9964 11.99118 0. 8490 11.36116 O.0534
0.9985 11.9952 0.9911 O.99119 O.99118
Results, normalized with respect to the exact solution, are summarized in Table 8. Once again all proposed plane elements yield identical results with the Pian and Sumihara [13] element. Almost no deterioration with v is observed, only 0.77% difference between the two extreme values of v. The QBI element yields almost identical results. The 2 x 2 element exhibits the well known locking at the nearly incompressible limit. 6.6. Finite strip with a hole This problem is introduced in order to test the accuracy of the stress distribution inside the proposed plane elements. A strip of finite width with a hole, of width b, is subjected to uniform tensile stress in the axial direction, as shown in Fig. 12. The axial stress along the critical section, line A - B , is investigated. First, the relatively simple case of plane stress with free upper and lower edges is investigated in Section 6.6.1. Then, the more challenging case of plane strain, at the nearly incompressible limit, with the upper and lower edges fixed in the vertical direction, is tackled in Section 6.6.2. 6.6. I. Finite strip with a hole: plane stress The material properties used are E = 1.0 and v = 0.3. First, convergence in the energy norm is examined. No analytical solution for the energy is known. Therefore, a converged finite element solution of 19.37 obtained for the Pian and Sumihara [13] element using 3072 elements is used as the reference solution. Results normalized with respect to this reference solution are shown in Fig. 13. The axial stress distribution along the critical section for the meshes containing 48 and 192 elements is presented in Fig. 14. The results are compared with a reference solution given by Savin [30]. 5.0
L I"
,t
Symmetry or Free .~-
B
tu
b -- 10/3
~
P = 1.0
A
Symmetry Fig. 12. Strip with a hole problem.
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
148
1.01
1.00
0.99 ¢¢
0.98
~L~ ~
"¢1
-'"
' elm m
0.97
-*I
E Z
:/
0.96
0.95
""
PSS~-S
- - .-. u
2X2
''"'"
0131
-""--
EXACT
/ I
,
3
; 1
4
'
,
192
768
, 3072
nei
Fig. 13. Strip with a hole; plane stress; convergence in the energy norm. 3.5 I= O
omm
e,m
z,.. r,~
2.5 ' ~
imm
- '
EXACT
"'
PSS/48
~'
r~ z.. rj:
1,5
Ilmm
M ,<
0,5
0:4
o'.s
o'.s
o.r
0.8
0:9
',.o
x2..Ib
Fig. 14. Strip with a hole; plane stress; axial stress distribution along the line A-B.
Almost identical results are obtained for the proposed elements and the Pian and Sumihara element (about 5.0E-3 difference in energy norm for the three-element mesh). A very small sensitivity to the incompatible shape functions used is observed (about 2 . 0 E - 3 difference in the energy norm for the thee-element mesh). Excellent results are obtained for both the energy norm, 95.82% of the reference solution with only three elements, and stress distribution inside the elements. The 2 x 2 and QBI elements yield comparable results for this problem. 6.6.2. Finite strip with a hole: plane strain The material properties used are E = 1.0 and v = 0.4999. First, convergence in the energy norm is examined. No analytical solution for the energy is available. Therefore, a converged finite element solution of 2.389 obtained for the Pian and Sumihara [13] element using 12288 elements is used as the reference solution. Results normalized with respect to this reference solution are shown in Fig. 15. The axial stress distribution along the critical section for the meshes containing 48 and 192 elements is presented in Fig. 16. No analytical solution for the
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
149
1.25
1.00 "
I...
0.75
¢: (1#
m
¢!,)
"
0.50 ~
N
em iml
,-,,
E I,,,
0.25
Z: ;' . . . . . . 3
000
m
PSS2,4,6
m ,=.m ,
0131
='''"='
2X2
T .... 12
PSS1,3,5
= "
•, . . . , ,
P S
,P °°t .,o °~'~
' " , " " = " " " " " ' = =l ' 1 ' b = = ! 48 192 768
I
3072
nel
Fig. 15. Strip with a hole; plane strain; convergence in the energy norm. 2.5 ¢:: O
I,m q[ml
:= ,,¢:
I ~
ellml
I. qlm) r4~
2.0"I
"'''''
PSS/48 PSS/192
•- = - "
Smoothed/768
" " ~
""
emll
r/j r/)
Im r/) alto elm
;I<
.<
,o
"o:4"
o:s" o:s" o'7 "o'8 " o : o ' , o
x21b
Fig. 16. Strip with a hole; plane strain; axial stress distribution along the line A-B.
stress distribution is available. Therefore, the smoothed stresses obtained for the 768-element mesh are used as a reference. The results obtained show the Pian and Sumihara element to exhibit a marginal superiority over the proposed elements for the coarse meshes, only 0.8% for a twelve-element mesh. The formulation shows a small sensitivity of the type of incompatible shape functions used, only 1.3% and 0.04% for three- and twelve-element meshes, respectively. Excellent results are obtained; with only three elements in the mesh, 69% and 68% of the reference solution is obtained for the PSS1,3,5 and PSS2,4,6, respectively. Excellent results are also obtained for the stress distribution inside the elements.
7. Concluding remarks The main objective of this paper is to apply the method proposed by Weissman and Taylor [1] to generate assumed fields for four-node quadrilateral linear elastic plane stress/strain
150
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
elements. O n e of the method's main goals is to avoid locking, which is commonly observed in finite element models of problems involving internal constraints. As an illustration of the treatment of internal constraints, the problem of plane strain at the nearly incompressible limit is chosen due to the difficulty it presents to finite element methods, and the large body of work which has been directed toward eliminating this obstacle. The application of the p r o p o s e d m e t h o d to nonlinear problems will be addressed in subsequent papers. Both elements formulated via the H e l l i n g e r - R e i s s n e r and H u - W a s h i z u variational principle are considered. For both types of formulations, it is proven that as the nearly incompressible limit is approached (plane strain), the trace of the strain vanishes pointwise (i.e., the internal constraint is satisfied pointwise). Consequently, locking at the nearly incompressible limit is avoided. The proposed elements are subjected to an extensive set of problems reported in the literature. Excellent resul~.~; are obtained for all examples. In fact the results are equivalent to those obtained by the best known elements (e.g., [13]). F u r t h e r m o r e , the variationaUy consistent stresses recovered at the element level yield excellent results. Identical results are obtained for elements formulated via the H e l l i n g e r - R e i s s n e r and H u - W a s h i z u functionals, provided the same incompatible functions are used. Also, identical results are obtained for both approaches taken to transform the assumed strain field from the element's natural space into the physical space. F u r t h e r m o r e , results show no sensitivity to the incompatible shape functions used.
References
Ill S.L. Weissman and R.L. Taylor, Treatment of internal constraints by mixed finite element methods: Unification of concepts, Internat. J. Numer. Methods Engrg. 33 (1992) 131-141.
121 T.J.R. Hughes, The Finite Element Method (Prentice Hall, Englewood Cliffs, NJ, 1987). [31 N. Bicanic and E. Hinton, Spurious modes in two-dimensional isoparametric elements, Internat. J. Numer. Methods Engrg. 14 (1979) 1545-1557.
141 T.J.R. Hughes, Equivalence of finite elements for nearly incompressible elasticity, J. Appl. Mech. 44 (1977) 181-183.
151 D.S. Malkus and T.J.R. Hughes, Mixed finite element methods-reduced and selective integration techniques: a unification of concepts, Comput. Methods Appl. Mech. Engrg. 15 (1978) 63-81.
[61 T.J.R. Hughes, Generalizing of selective integration procedures to anisotropic and nonlinear media, Internat. J. Numer. Methods Engrg. 15 (1980) 1413-1418.
171 J.C. Simo, R.L. Taylor and K.S. Pister, Variational and projection methods for the volume constraint in finite deformation elasto-plasticity, Comput. Methods Appl. Mech. Engrg. 51 (1985) 177-208.
181 D.P. Flanagan and T. Belytschko, A uniform strain hexahedron and quadrilateral with orthogonal hourglass control, Internat. J. Numer. Methods Engrg. 17 (1981) 679-706.
191 T. Belytschko and W.E. Bachrach, Efficient implementation of quadrilaterals with high coarse-mesh accuracy, Comput. Methods Appl. Mech. Engrg. 54 (1986) 279-301.
[1o1 E.L. Wilson, R.L. Taylor, W.P. Doherty and J. Ghaboussi, Incompatible displacement modes, in: S.J. [11] [12] [23] [14] [15]
Fenves, N. Perrone, A.R. Robinston and W.C. Schnobrich, eds., Numerical and Computer Models in Structural Mechanics (Academic Press, New York, 1973). R.L. Taylor, P.J. Beresford and E.L. Wilson, A nonconforming element for stress analysis, Internat. J. Numer. Methods Engrg. l0 (1976) 1211-1219. G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice Hail, Englewood Cliffs, NJ, 1973). T.H.H. Pian and K. Sumihara, Rational approach for assumed stress finite elements, Internat. J. Numer. Methods Engrg. 20 (1984) 1685-1695. T.H.H. Pian and C.-C, Wu, A rational approach for choosing stress terms for hybrid finite element formulations, Internat. J. Numer. Methods Engrg. 26 (1988) 2331-2343. C.-C. Wu, M.-G. Huang and T.H.H. Plan, Consistency condition and convergence criteria of incompatible
S.L. Weissman, R.L. Taylor, A unified approach to mixed FEMs
151
elements: General formulation of incompatible functions and its application, Comput. & Structures 27 (1987) 639-644. [161 R. Rubinstein, E.F. Punch and S.N. Atluri, An analysis of, and remedies for, kinematic modes in hybrid-stress finite elements: selection of stable, invariant stress fields, Comput. Methods Appl. Mech. Engrg. 38 (1983) 63-92. [17] E.F. Punch and S.N. Atluri, Development and testing of stable, invariant isoparametric curvilinear 2- and 3-D hybrid stress elements, Comput. Methods Appl. Mech. Engrg. 47 (1984) 331-356. [181 C.-T. Yang and S.N. Atluri, An assumed deviatoric stress-pressure-velocity mixed finite element for unsteady, convective, incompressible viscous flow: Part l--Theoretical Development, Internat. J. Numer. Methods Fluids 3 (1983) 377-398. 119] J.C. Sima a~ld M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Internat. J. Numer. Methods Engrg. 29 (1990) 1595-1638. [2o1 B. Fraeijs de Veubeke, Displacement and equilibrium models in the finite element methods, in: O.C. Zienkiewicz and G.S. Holister, eds., Stress Analysis (Wiley, London, 1965). [211 J.C. Simo, J.G. Kennedy and R.L. Taylor, Complementary mixed finite element formulations of elastoplasticity, Comput. Methods Appi. Mech. Engrg. 74 (1990) 177-206. [221 O.C. Zienkiewicz, S. Qu, R.L. Taylor and S. Nakazawa, The patch test for mixed formulations, lnternat. J. Numer. Methods Engrg. 23 (1986) 1873-1883. [231 I. Babuska, J.T. Oden and J.K. Lee, Mixed-hybrid finite element approximation of second-order elliptic boundary-value problems, Comput. Methods Appl. Mech. Engrg. 11 (1977) 175-206. [24] W.M. Xue and S.N. Atluri, Existence and stability and discrete BB and rank conditions, for general mixed-hybrid finite elements in elasticity, in: R. Spilker et al., eds., Hybrid and Mixed Finite Element Methods, AMD Vol. 73 (ASME, New York, 1985) 377-398. [25l P.G. Bergan and M.K. Nygard, Finite element with increased freedom in choosing shape functions, Internat. J. Numer. Methods Engrg. 20 (1984) 643-663. [26] R.L. Taylor, O.C. Zienkiewicz, J.C. Simo and A.H.C. Chan, The patch test for mixed formulations, Internat. J. Numer. Methods Engrg. 22 (1986) 32-62. [27] R.H. MacNeal and R.L. Harder, A proposed standard set of problems to test finite element accuracy, Finite Elements Anal. Design l (1985)3-20. [281 R.D. Cook, A plane hybrid elemen~ with rotational d.o.f, and adjustable stiffness, Internat. J. Numer. Methods Engrg. 24 (1987) 1499-1508. [29] R.H. MacNeal, A theorem regarding the locking of tapered four-noded membrane elements, lnternat. J. Numer. Methods Engrg. 24 (1987) 1793-1799. 1301 G.N, Savin, Stress concentration around holes, International Series Monographs in Aeronautics and Astronautics (Pergamon, Oxford, 1961). [311 O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, 4th Edition (MacGraw-Hill, London, 1989).