Extrapolated fields in the formulation of the assumed strain elements Part II: Three-dimensional problems

Computer methods In applied mechanics and engineering Comput.

ELSEVIER

Methods

Appl. Mech. Engrg.

154 (1998)

l-29

Extrapolated fields in the formulation of the assumed strain elements Part II: Three-dimensional problems Yung-I. Chen’, Henryk K. Stolarski* Depanment of Civil Engineering, University of Minnesota, Minneapolis, MN 55455, USA Received

16 June 1994; revised 8 August

1996

Abstract Assumed strain eight-node hexahedral elements with significantly extended range of applicability are presented. These elements are formulated using only standard translational displacements at each of their eight nodes and provide accurate solutions for a variety of benchmark problems such as spatial beams, plates, shells as well as general three-dimensional elasticity problems. The elements with such characteristics are particularly useful to model problems of complex geometry, where different regions of the domain might impose entirely different modeling requirements. The formulation starts with introduction of a parallelepiped domain associated with the given eight-node hexahedral element. Then, the assumed strain field is constructed for that associated parallelepiped domain. It is done by identification of various modes of its deformation and by proper modification of the strain field in the constant and linear bending modes. Strain and displacement extrapolation from the associated parallelepiped to the original hexahedral domain is subsequently used to establish the assumed strain field for the given element. Solutions to some popular benchmark problems demonstrate that the proposed assumed strain hexahedral elements exhibit remarkably high accuracy even when severely distorted and high aspect ratio meshes are used. Another advantage of the present elements is that locking for nearly incompressible materials is also mitigated. Unfortunately, the elements pass the patch test only when their shapes are parallelepipeds.

1. Introduction The eight-node hexahedral element is perhaps one of the simplest and most widely used three-dimensional elements. One reason behind its popularity is that-despite its simplicity-the element possesses basic features needed to represent the elementary modes of deformation such as rigid body, constant strain, constant bending, and some linear bending modes. Combination of these modes of deformation should theoretically allow the eight-node hexahedral element to be applicable not only to the general three-dimensional problems but also to dimensionally degenerate situations typical for structural members such as beams, plates, and shells, where the bending properties are important. Development of the elements with such extended range of applicability is important in the analysis of geometrically complex problems. In those problems, different zones of the domain may have different geometric proportions and resulting different response characteristics and modeling requirements. Use of a single element that would accurately model all those zones is highly advantageous.

* Corresponding author. ’ Presently at Provisional Engineering North Road, Taipei, Taiwan 105.

Office of High Speed Rail, Ministry

0045-7825/98/$19.00 0 1998 EIsevier Science S.A. All rights reserved PII SOO45-7825(97)00084-4

of Transportation

and Communications,

No. 240, Tun-Hwa

2

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

Unfortunately, as is well known, standard displacement formulation of the eight-node hexahedral element includes parasitic shear terms which render unacceptable results in structural analysis. In the last decade, a number of papers have been written to improve performance of the eight-node hexahedral in those special applications. Representative techniques include incompatible modes, drilling degrees of freedom, mixed method, and reduced integration method. In this paper, a new approach to the problem is discussed. It is an extension to the three-dimensional problems of a special version of the assumed strain method introduced for two-dimensional analysis in [l]. To provide a historical context for the development presented here, we start with the work of Wilson et al. [2] who suggested an incompatible eight-node hexahedral element to improve the bending behavior of the standard compatible isoparametric eight-node element. Nine additional internal displacement parameters were introduced to achieve this improvement. A modified version of the incompatible modes method which allows to satisfy the patch test was introduced by Taylor et al, [3]. Recently, the concept of incompatible modes as well as the field-consistent scheme was adopted by Chandar and Prathap [4] to formulate an eight-node displacement based hexahedral element. Both formulations are still sensitive to mesh distortion and lead to poor stresses prediction for certain problems. Very recently, various eight-node hexahedral elements with additional drilling degrees of freedom have been successfully formulated by Yunus et al. [5], Ibrahimbegovic and Wilson [6], Sze and Ghali [7], and others. A common feature of the above-mentioned elements is that they are all derived starting from a twenty-node straight-sided hexahedral element and employing the original Allman’s displacement interpolation scheme [8]. In addition, penalty stiffness matrices are often required to stabilize spurious energy modes which are intrinsic to the Allman’s displacement interpolation. Elements presented in [6] and [7] are formulated on the basis of a variational principle employing the drilling degrees of freedom as independent variables. Eight-node hexahedral elements employing additional drilling degrees of freedom at comer nodes are shown to exhibit good bending behavior. However, the computational time required to use this kind of elements in the analysis of three-dimensional problems will be increased significantly since the number of degrees of freedom per node is doubled, as compared to the elements with only three displacement degrees of freedom at each node. The origins of developing hybrid assumed stress elements can be traced back to the work of Pian [9]. Since the concept was introduced in [9], Pian and his colleagues [lo-141 as well as others have devoted a great deal of attention to formulate successful hybrid assumed stress elements. Using the Hellinger-Reissner variational principle, Pian and Chen [ 1 l] presented an assumed stress field involving eighteen independent stress parameters which are required in order to suppress all kinematic modes of the eight-node hexahedral element. A similar idea was also used by Loikkanen and Irons [ 151. Pian and Sumihara [13] established an assumed stress field with five independent stress parameters to formulate a four-node, assumed stress, two-dimensional element based on the functional of [lo]. This four-node variationally consistent element can be viewed as the representative in the field of developing the assumed stress elements. Taking a variational point of view, Pian and Tong [14] have examined the relationship between the hybrid stress formulations and incompatible mode approach. The results of their work have confirmed the existence of a close relationship between the two approaches, that was originally pointed out by Irons [16]. In addition, methods used in [ 131 are extended by Pian and Tong [ 141 to establish an independent stress field with eighteen stress parameters for formulating the eight-node hybrid hexahedral element. By using an alternative approach to determine independent stress field, Dong and Freitas [17] presented an eight-node hexahedral element based on the modified functional of [lo]. The independent stress field of [ 171 is expressed in terms of twenty-four stress parameters. This element exhibits good accuracy for both displacements and stresses. By taking the hybrid hexahedral element of [ 141 as parent element, Sze [ 181 devised recently what he calls ‘admissible formulation’ and the orthogonal stress/strain interpolants to derive the eight-node hexahedral element. Numerical results show that accuracy of Sze’s hybrid elements is close to their parent elements. By extending their previous work of [ 181 and by using the so-called shear scaling factors, Sze and Ghali 1191 presented a hybrid hexahedral element, which can be used in the analysis of beams, plates, shells and three-dimensional elasticity problems. In their work, different selections of shear scaling factors are required for different types of applications. For the analysis of thin plate/shell problems, remarkably accurate results have been obtained when shear scaling factors are selectively used. Their element is also capable of avoiding locking associated with nearly incompressible materials. By comparing the numerical results reported in [ 191 and in [ 141

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) 1-29

3

one can note, however, that, without use of the shear scaling factors, the performance of their element is identical to its parent element which is the eight-node hybrid-hexahedral element introduced by Pian and Tong [14]. On the other hand, for problems of complex geometry where different regions may be characterized by different type of response, selective use of the shear scaling factors is, in fact, equivalent to using different types of elements in different zones of the domain (e.g. shell elements, three-dimensional elements, etc.). A better situation would exist if a single three-dimensional element was designed so as to facilitate a smooth transition in the behavior of the model without any additional factors dependent on the geometric properties of the element. A new functional involving three independent variables of displacement, stress and strain was recently presented by Chen and Cheung [20]. Based on this functional, alternative approaches to formulate the eight-node hybrid hexahedral elements have been proposed by Chen and Cheung [21-231 and others. Elements derived from the functional of [20] have excellent bending properties and improved accuracy of stresses. However, they are still sensitive to the use of distorted meshes, From a computational efficiency point of view, reduced integration method has been used quite successfully. Its effectiveness is particularly noteworthy in application to three-dimensional solid elements. However, this is somewhat offset by the presence of the so-called hourglass modes. In the last decade, a number of developments have been made to deal with hourglass modes problems caused by the use of reduced integration scheme. A method termed y-projection technique for isolating and controlling the hourglass modes has been introduced by Flanagan and Belytschko [24] and specified further by Belytschko [25] soon thereafter. An alternative approach to control spurious modes, based on the Hu-Washizu variational principle, was used to formulate both under-integrated quadrilateral and hexahedral elements by Belytschko et al. [26]. A similar concept of y-projection method has been pursued by Liu et al. [27], Wang and Belytschko [28], Bachrach [29], and others. A new version of the assumed strain formulation has been presented by Stolarski and Chen [30] to improve performance of the four-node quadrilateral element. The basic idea of the formulation consists in identification of various modes of deformation of the element and then in proper modification of the strain field in the so-called in-plane bending modes. This element exhibits extremely accurate in-plane bending properties. However, the methodology used in its original formulation was not applicable to three-dimensional problems. Search for a formulation applicable to both two-dimensional and three-dimensional problems led to a modification presented in [ 11. The basic idea there consists in introduction of a parallelogram associated with a given general quadrilateral element, subsequent formulations of the assumed strain field in this parallelogram domain, and extrapolation of that field to the domain of the original quadrilateral element. It has been proved by means of various numerical results that element’s performance for both bending and membrane dominated problems is extremely accurate. An extension of the idea presented in [II to three-dimensional problems is discussed herein. As in the two-dimensional formulation, in the context of three-dimensional problems, a parallelepiped domain which is associated with a given eight-node hexahedral element is introduced first. The concepts of the mode decomposition method and strain modification in some of the resulting modes presented in [30] is then extended to establish the assumed strain field for the associated parallelepiped. The strain and displacement extrapolation analogical to that given in [l] for two-dimensional problems leads to the assumed strain field for the general eight-node hexahedral element. Just like the two-dimensional element presented in [l], the three-dimensional elements proposed in this work satisfy the patch test only- for special (parallelepiped) shapes. This was understood to be a drawback of the present formulation and various attempts had been made to correct it. Unfortunately, the resulting elements either failed to pass the patch test or yielded results which were inferior and less consistent in comparison with those obtained using the formulations described herein (which is in agreement with findings described for two-dimensional problems in 1331). Thus, rather than present a formulation that satisfies the patch test, but provides results similar to those obtained with many other existing formulations, it was decided to demonstrate numerically that accurate and convergent results can be obtained without satisfaction of the patch test. It was hoped that this presentation might motivate someone to improve the idea presented here further. An additional encouragement to do SO was the fact that in the solution of some problems (for example, shells or axisymmetric continuum problems) no general formulation of the patch test seems to exist, yet successful numerical solutions of those problems have been obtained. The concept of the patch test in those cases is assumed to be valid only in the limit, when the size of the elements tends to zero; this concept is in agreement with the original idea of the

4

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) 1-29

test. In that sense the patch test can also be satisfied for the formulations presented here if, for example, the mesh refinement is achieved by successive subdivisions of the original mesh along isoparametric lines; in the limit this leads to parallelopiped elements. Existence of the problems which allow to satisfy the patch test only in the limit has also been acknowledged, for example, in [31]. However, doubts have been expressed there (without any supporting arguments) as to the practical engineering use of the elements satisfying the patch test only in the limit. While for some elements this skeptical opinion may be verified, its general validity seems impossible to assert. In fact, the results of a comprehensive set of tests, used to evaluate the formulations described here, seem to contradict that opinion. This was another argument behind the decision to present elements developed in this work. The outline of this paper is as follows. In Section 2, a parallelepiped domain associated with a general eight-node hexahedral element is introduced. The modes of deformation of the associated eight-node parallelepiped, important for its subsequent formulation, is discussed in Section 3. Construction of the assumed strain field in the associated parallelepiped is presented in Section 4. In Section 5, the concept of the strain and displacement extrapolation is introduced to obtain the assumed strain field of the given hexahedral element. Section 6 describes test problems and numerical results obtained. Conclusions form the last, seventh section of the paper.

2. Eight-node hexahedral element and its associated parallelepiped For a general eight-node hexahedral element, an associated parallelepiped shown in Fig. 1 is defined by its nodal position vectors X1 (of node i) given by the following expressions Xr = Xc + [Ge G,, G&,

iip,&>’ = $e,

(14

.

In the above equation, X, is the centre of the associated parallelepiped,

Fig. 1. A general eight-node

hexahedral

element and its associated

parallelepiped.

5

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

X,=:+gx,=x;e,,

(lb)

I-1

the kth components

of vectors G,, Gi and Gf are given by the following

expression

{&, iii, lf}’ is the vector of coordinates of the nodal point i in the space of the isoparametric parameters {& <, [} of the associated parallelepiped whereas (5,. q,, 5,) T is the coordinates of the nodal point I in the space of isoparametric parameters { ,$,q, J} of the original hexahedral element, X, is the position vector at node I of the hexahedral element, and subscript or superscript k represents kth component of the vector relative to the fixed global Cartesian system ek = ek (kth coordinate). Throughout this work, indices repeated on two different levels imply summation over their range, the superscript ‘T’ designates the transpose of a matrix, and a ‘-’ on the top indicates that a quantity is related to the associated parallelepiped, unless otherwise stated.

3. Deformation modes of the eight-node parallelepiped For the associated eight-node parallelepiped, there are totally twenty four modes of deformation corresponding to twenty four degrees of freedom. They consist of six rigid body motion modes, six constant strain modes, and twelve additional higher-order modes corresponding to higher-order strain fields. The higher-order strain modes are the only ones in which strain fields have to be modified while seeking improvement in the performance of the eight-node hexahedral element. They can be further grouped into three linear shear strain modes, six constant bending modes, and three linear bending modes, see Figs. 2-4. It is understandable that, to improve bending behavior of the element, no strain modification is needed for the strain field corresponding to linear shear strain modes. However, identification of these higher-order shear strain modes is also important to design successful eight-node hexahedral elements. Proper identification of the higher-order bending modes along with the corresponding modification of the strain field in those modes is particularly crucial to formulate an effective hexahedral element with improved bending properties. Before proceeding to the mathematics of the mode decomposition process, we will briefly describe the deformation modes involved. 3.1. Constant strain modes The rigid body and constant of parameters

strain modes of the associated parallelepiped

are represented

by the following

set

(2) where Lib, i ;, and li i are, respectively, translational displacements of node 1 in e’, e2 and e3-direction, @i is the rigid body rotation of the element about global e’-direction at node f, and it represent a constant part of the strain field. 3.2. Linear shear strain modes The linear shear modes are represented

by the following

parameters

nLs= {S,, s*, SJT.

P

In the above equation, the subscripts indicate that the deformation mode the corresponding parameters are related to will cause no shear deformation on the pair of faces which will be intersected by the vector G with the same subscript (cf. Eq. (1~)). For instance, the linear shear mode described by Sp is designed in such a way that ,.A_.. AA_,. faces 1485 and 2376 (see Figs. 1 and 2), which will be intersected by vector G,, deform with no change in length along the edges and diagonals of those two faces (they just warp). Thus, the nodal displacement vectors have to be perpendicular to those faces. However, the remaining four faces will undergo shear deformation, their

Y.-I. Chen, H.K. SIolarski I Cornput. Methods Appl. Mech. Engrg. 154 (1998) l-29

Fig. 2. Deformed

configurations

for linear shear modes s,, s, and si (for a cubic shape of the hexahedral).

diagonals will change in length while the length of the edges will remain unchanged. defined similarly. 3.3. Constant

The modes S+ and St are

bending modes

For the deformation

modes related to the constant bending

BCB ={{PZB}=, {g;“}‘, where subscripts

{P;“,“,}’

$ 6 and l denote,

={{&e q, respectively,

strain field, the following

set of parameters

{&, p,>, {F@ ?$>)‘, the axes of rotation

is used (4)

with the positive

direction

defined by

Y.-l. Chen, H.K. St&u-ski

Fig. 3. Deformed

configurations

of constant

f Comput. Methods Appl. Mech. Engrg. 1.54 (1998) I-29

bending

modes (~2~ &,}, {p,, Bi} and { $Q, j$} (f or a cubic shape of the hexahedral).

vectors G& G+ and Gi. The parameters listed above and named using the same letter (6, @ or y) are considered to be a pair and, accordingly, are grouped together so that they can participate in the deformation process simultaneously (more detailed explanation will be provided in the sequel). The constant bending mode h’s, for instance, represents the rigid body rotation of the pyramid &%% about the vector Gi while keeping pyramid *Ah*Cl265 fixed. This deformed configuration is illustrated in Fig. 3 (compare Fig. 1). Similar explanation can be given for the other five constant bending modes.

Y.-l. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

Fig. 4. Deformed

3.4. Linear bending

configuration

of linear bending

modes t$, 8+ and Bi (for a cubic shape of the hexahedral).

modes

The remaining three higher-order deformation modes are used to describe the linear bending and are represented by the following three parameters

configurations

PALB= (d,, B,, Be,, . The linear bending

mode described by 8 with a specified subscript (8, ?j or d) will produce deformation

(5) causing

9

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

no shear deformation on those two faces which will be intersected by the vector G with the same subscript; however, the remaining four faces will undergo ‘in-plane bending’. This is illustrated in Fig. 4. Modification of the strain field in these three linear bending modes and in the six constant bending modes are necessary for the associated parallelepiped to adequately simulate various bending states in a three-dimensional deformation process. The twenty four deformation modes participating in the deformation process and described in the previous paragraphs will now be specified quantitatively. This leads to the following equations representing the relationship between the parameters representing those deformation modes and the vector of element nodal displacements a=aC

+dLS

+aCB

+aLB

=

ep’

+

fiL”“pL”

+

&BpCB

+

fiL”pLB (6)

In

the

above equation, d is the vector of total element nodal displacements

of the associated parallelepiped

given

by a = {;a, ti;, . .,li;

)...,

and the remaining

symbols

3.5. Quantitative

specijcation

i;

ZQ;}‘,

will be described

= {l;),l;;,1;;>

Ua,b)

in the sequel.

of the deformation

modes

Matrix c of Eq. (6) is a 24 X 12 matrix related to rigid body and constant strain modes. It contains constants dependent on the nodal coordinates and is partitioned into eight 3 X 12 matrices corresponding to eight nodal points of the associated parallelepiped as follows

e = [CT,e;,...,e;,...,e;,CilT

ci=V,,,,

A/i,/,

Aiiun =

- $)eakn

C$

Aif ,,en, Aii3,,en, Bii,e’, Bpi2e2, Bii3e3, Bii2e’, Bii3e2, Bii,e3], ,

Btj,, = (Xy - Xl;) 9

(8b-d)

where I, x 3 is a 3 X 3 unit matrix and the symbol eakn stands for the component of the permutation tensor. Matrix k”” m . Eq. (6) is a,2,4 X 3 tnL$rix corresponding to three linear shear modes. In order to define the of Eq. (6), unit vectors t,, with M = 1,2,3, are constructed at every relationship between vectors d and P node i of the associated parallelepiped g the directions parallel to vectors G,, G+ and G@ respectively; and is composed of eight 3 X 3 matrices. Those 3 X 3 matrices for node / designated as shown in Fig. 5. Matrix H of the associated parallelepiped are given by the following expressions k’”

= [ti2 X tjj, t, X ffi 7tf, X tiz] = [ti, t;, t;] .

(9)

Matrix fi cB of Eq. (6) is a 24 X 6 matrix which corresponds to the six constant bending modes, its eight parts corresponding to nodes i are given by the following formulas * CB Hi

A CB =

[[Hiirg

. CB

H/&ilt

(104

where (Xt - Xs)Gye,k,en 0,

,

for i listed in Table 1 for other nodes i

1’

with c$~ representing one of the rotations listed in Eq. (4) about direction indicating the point e with the position vector XC specified in Eq. (lb). Matrix fiLB m Eq. (6) is a 24 X 3 matrix corresponding 3 X 3 matrices which, for each node i of the associated

(lob) ci = z, Q or l, and subscript

‘6’

to three linear bending modes. It is composed of eight parallelepiped, are given as follows:

(11)

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) 1-29

10

Fig. 5. Designation

of the unit vectors tiM used in the definition

of the assumed

strain field and matrices 8””

and a””

Table 1 Nodes i for which &FL # 0 Constant

bending

mode 4:

n

where subscripts 8, +j and C are referred to linear bending identical as those defined in Eq. (9).

3.6. Projection

A n modes 13~ 0, and i$, respectively;

and vectors tr are

operators

Inversion of Eq. (6) provides quantitative information regarding overall deformation process of the associated parallelepiped

{{PC}‘, {p”“}‘, {p”“}‘, {p”“}‘}’

= [&, fiLS,fi”” &““I_‘ci=

participation

[[k”]‘,

of the particular

[ff”“]‘,

[ZF”]‘,

modes in the

[k”“]‘]‘ci

(12)

which leads to

pc=pc(j,

p

= p”“(j )

P~cB~{{~~B}‘,{~~B}~,{~~B}‘}‘=~cBa=[[~~B]T,[~~B]T,[~~B]T]T(i, PnLB= {‘g,, (j,, Jt}T = ff

LB(j= [{k k”}‘, {@L,B}T, {ii ;“,‘]‘d .

(13a,b) (13c) (13d)

Combining the above results with Eq. (6) allows one to describe the element nodal displacements in the constant strain, linear shear strain, constant bending and linear bending modes in terms of a linear combination of the total element nodal displacements

11

( 140) (14c) (14d) ..CB

where f” and f’“” are, respectively, called the constant strain and linear shear strain projection operators, r, , with * = &, b or 9, represents the constant bending projection operators for constant bending modes {c$ G,>‘, {p,, fi,}‘, or {Tf, Tt}T, respectively; pk”, with * = {, 6 or [, stands for the linear bending projection operator with respect to linear bending modes 8,, di or $, respectively. The projection operators satisfy the following conditions ~C+~LS+~&B+~~B+~:B+~;B+~~+~ZB=I, FOB

*

=

#

fi,

forA=Band

0,

any other combination,

(15a)

* =#,

(15b)

where Z is the identity operator. The projection operators in the above equation determine the decomposition the deformed configuration of the associated parallelepiped into particular modes of deformation.

4. Assumed

strain field for the associated

Discretized form

standard

strain-displacement

of

parallelepiped relationship

for the associated

parallelepiped

takes the following

E=tici,

(16)

where fi is the so-called

strain-displacement

Z=fi~c+fi&LS+~~cB+&B

operator.

Eq. (16) can be decomposed

into four parts as follows

= gc + &LS + &cB + sLB,

(17) . LS ACB and kLB are the unmodified strain fields associated with constant strain, linear shear strain, where 2, E , E constant bending, and linear bending modes, respectively. The strains identified in Eq. (17) result from the usual strain-displacement relationship and standard isoparametric formulation. While no modification is required for the constant strain and linear shear strain fields, improvements in the bending behavior of the associated parallelepiped crucially depends on how the strain fields in the constant and linear bending modes are modified compared to those of Eq. (17). This modification leads to an assumed strain field, and the corresponding strain-displacement operator for the associated parallelepiped. The technique of constructing the assumed strain formulation of the associated parallelepiped is outlined below. The concept of modification of higher-order strain fields ScB and iLB in Eq. (17) for the associated parallelepiped is similar to the concept adopted in [l] for the four-node parallelogram. To modify higher order strain fields 6 cB and S LB m . Eq. (17), three new unit vectors for each node i have to be constructed in addition M = 1, 2,3, present in Fig. 5. These three additional unit vectors tiM, with M = 4, 5 and 6, are constructed to ti‘+f, along the diagonals of the faces defined by the vectors (t f,, ti2), (tf,, ti3) and (frZ, ti3), respectively. At each node f, the strains A,, are evaluated along each of those six directions instead of the usual strain components. They equivalently and completely describe the strain state at the nodes and are easier to use in the following modification of the higher-order strain fields due to the constant bending and linear bending modes. 4.1. Constant strain The formulas &=&,8, “C

modes

describing M=

constant

strain field in six different

directions

of each node i are

1-6,

where A f,+,can be derived directly by considering member).

(18) elongation

of the line under consideration

(as if it was a truss

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

12

4.2.

Linear shear strain modes

Since no strain modification is needed for the strain field corresponding to the linear shear modes, and since in those modes each face of the parallelepiped undergoes homogeneous in-plane straining, the linear shear strains ik” can be assumed to be (19) 4.3. Constant bending modes In the three-dimensional problems, it turned out beneficial to group properly constant bending modes and modify the strain field in the entire group rather than in each mode individually. When the modification is performed individually and superposed together afterwards, they interfere negatively and result in a somewhat inferior assumed strain field, particularly for nearly incompressible materials. A similar situation occurs when the selection of the modes in the group is incorrect and does not properly reflect what happens in the actual bending process. For those reasons, the three groups identified in Eq. (4) by different symbols (&, b, T) have been selected here. Each group contains the deformation resulting from two constant bending modes which will occur simultaneously when Poisson ratio v # 0 (similarly to anticlastic curvature in plates). In each of these groups, the strain field has been modified independently of the modifications in the remaining two groups. The essence of those modifications is to create two-dimensional state of stresses in the plane perpendicular to the faces of the two rotating pyramids, associated with both modes of each constant bending group. Further clarifications of those modifications are provided in Remark 4.1; because this remark refers to particular terms of the following equations, it is placed at th,eccnd of this section. The modified constant bending defined as follows:

strains A!,,,, (in M-direction)

at each node i are now split into three parts and

(204 where subscripts * = &, /? or 9 are, respectively, used to indic:kte_ a quantity related to the three groups of ah, constant bending modes {&& c?+}~, {p,, fi,}’ or {T@ ytjT, and ArM is a modifying term used together with Ai,,, identified in Eq. (18) to more adequately simulate the strain state gB the constant bending modes and to improve the element’s performance for nearly incompressible materials. A f,,, is defined similarly for all constant bending modes, but the same modifications are related to different directions M depending on which of the constant bending modes is considered. To unify and shorten the presentation of those modifications, the six involved directions are c$llld M = 0, P, Q, R, S, T and their identification for particular modes is given in Table 2. With this notation, At,,,, for all constant bending modes can be defined by the following formulas

Table 2 Variables

0, P, Q, R, S and T for each bending

Direction

Constant

M

68

0 P

1 2

Q

3

R

4 5 6

s T

bending

mode Linear bending

mode

mode

iy,

4

8;

4

4

@E

“.+

i$

1

2 3

2 3

I

6 4 5

6 4 5

3 1 2 5 6 4

1 2 3 4 5 6

2 3 1 6 4 5

3

I

3 1 2 5 6 4

2 3 4 5 6

I 2 5 6 4

13

Y.-l. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) I-29

(20c-e)

where tf is a vector orthogonal to the appropriate face of the parallelepiped identified in Eq. (9), i = v/( 1 - v), and Y is the Poisson’s ratio.

and identical

with one of the vectors

REMARK 4.1. As explained in the text following Eq. (4), the constant bending modes represent rotations of the appropriate pyramids about the vectors G given in Eq. (1~). As a result, the principal strains within the faces of the parallelepiped undergoing compression or tension in those modes have directions specified by Eq. (20d). Both brackets of Eq. (2Ob), which correspond to the two modes considered as a group (C;, fi or p), are designed so as to create one-dimensional state of strains along this direction. The first term in each bracket defines strains in the direction of tfM accompanying the one-dimensional strain along the principal direction defined by Eq. (20d). The second term in each bracket accounts for the fact that, if only one of the two modes forming the group is present, the plane states of strains exists. Under these conditions, the corresponding terms representing strains perpendicular to the faces strained in the mode considered should be added if use of the general three-dimensional constitutive equation is to lead to one-dimensional state of stresses. Finally, since Eq. (20d) specifies a vector which, in general, is different from any of the six vectors tiM associated with the nodes of the parallelepiped, coefficient SzL of Eq. (20e) is introduced to calculate the strain in the principal direction defined by Eq. (20d) in terms of the strain along one of the vectors tfM, M = 1,2, 3, which is close to the principal strain direction. 4.4. Linear bending modes Since in the linear bending modes two adjacent faces of the parallelepiped undergo simultaneous ‘in-plane bending’, the modifications present in the following equations are designed so as to create one-dimensional state of strains along the intersection of those faces. This modification gives the same value of the strain along each corner line (intersection of adjacent faces) no matter which of the adjacent facyLbhis line is approached from. In the case of the linear bending modes, the modified linear bending strains hi are defined for each of the three modes separately and superposed as follows:

(214

# where subscripts

# = $, +j or [ stands for a quantity corresponding

to one of the three linear bending

modes 6,,

8, or t!$, respectively; /i FMis identified in Eq. ( 18), and R Fi is the modifying term introduced to improve linear bending properties of the element including the effects of the Poisson’s ratio. Assuming the same notation that has been used in the presentation of the constant bending modes and the corresponding identification of the directions given in Table 2, the modifications for all three linear bending modes are written in a unified fashion as follows:

,. LB

A

iM

,. LB /iiN

=

[(tlo . tlM)* -

AC lo

v(ifo. tpM>*lA

,

AC

=

[(tfo.tiN)* - ~(i;~~t,)*lA,,

M = f’, R >

R ;I = [(tfo . t,)* - z&g . t&)* - v(i& . t&*1& where definition

@lb)

N=Qvs> ,

of tl”, is given in Eq. (20d).

REMARKS 4.2. Modifications introduced in Eqs. (20b) and (21b) involve the Poisson’s ratio v; thus, already at this level, they incorporate constitutive equations. For low-order elements, like the one discussed here, this seems to be necessary to improve their performance in a wide range of applications, including in-plane bending and incompressibility. To stay focused on the underlying mechanisms of the problem, the constitutive equations used here are those of linear and isotropic elasticity. However, other material models can be incorporated, as

14

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) 1-29

long as strains can be assumed to be small. In particular, anisotropic elasticity can be easily introduced by using appropriate material constants instead of the Poisson’s ratio Y in Eqs. (20b) and (21b). Nonlinear material models can also be implemented within the framework presented here. To this end tangent material properties should be used within an incremental approach. However, all that goes beyond the scope of this work and may be investigated in the future. REMARK 4.3. As in the case of the constant bending modes, the first term in the brackets of Eq. (21b) associates strains in the direction of ti,,,, with the one-dimensional state of strains along the comer line. The second term accounts for the Poisson’s ratio effects, but this time a general three-dimensional case of deformation is considered since no plane strain case is identifiable. 4.5. Total modijied nodal strains The total modified nodal strains at each node i of the associated parallelepiped are obtained the unmodified constant strain part with the modified higher-order strain part. The result is 7 ..C A,=A,+A,+A,

..LS

,CB

by superposing

,LB

M= 1,2,3,4,%6,

+A,,

where subscript M represents the Mth vector t, emanating from node i Using Eqs. (14a,b,c,d), (15a), (18), (19), (20a) and (21a), the above equation is presented in the matrix form and expressed in terms of the vector of total nodal displacements of the associated parallelepiped

wheretheMthrowofhF,/iFB

and Af”

are identified

in Eq. (18), Eqs. (20a,b), and Eqs. (21a,b), respectively;

p:“, with * = 6, p, T, is specified in Eq. ( 14~) whereas p’,“, with # = $, $ l, is given in Eq. (14d). Note that in Eq. (23) all modifications are included at the nodal level. As can be seen from Eq. (23), it is accomplished by a technically simple (but, admittedly, difficult to explain) procedure of a$Jing additional six modifying matrices (corresponding to constant and linear bending modes) Att the matrix Ai. These modified matrices are in turn obtained by linear combination of parallelepiped. Eq. (23) provides modified strain state different directions of that node. They are transformed containing standard engineering strain components in $=s^;‘&=$;‘&&j$,

rows of the matrix Al at each node i at each node i by specifying modified back into modified usual strain vector the fixed global Cartesian system. The

of the associated strains along six & at each node i result is GW

where the entries of the matrix s^, are given as follows:

jI is the assumed strain operator at each node ! of the associated parallelepiped. With the above-modified strains defined at each node ! of the associated parallelepiped, an extrapolation procedure (similar to that in [l]) can be applied to obtain the assumed strain field for the original hexahedral element. This will be explained in the next section. Here, the following remarks are intended to provide some insight leading to the modifications specified in this section. REMARK 4.4. /i FMof Eq. (18) in a specified direction sfcthe unit vector tpM is related to only two nodal points defining that unit vector. Therefore, the components in A!,,., corresponding to the other six nodal points will be ,. LB 1 CB zero. Since matrices A p and A f of Eq. (23) are obtained from the matrix A,” (cf. Eqs. (20b) and (21b)), most of the entries of A;” and /it” are zero (117 entries and 114 entries out of 144 entries of these two matrices, respectively). Exploiting this, the computational cost associated with the evaluation of the strain-displacement matrix lf in Eq. (23) can be reduced

substantially.

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29 5.

Strain and displacement extrapolation

5.1. Strain extrapolation The assumed strain field iii of Eq. (24a) is evaluated at the nodal points of the associated parallelepiped and expressed in terms of its nodal displacements. Therefore, in order to obtain the assumed strain field of the original hexahedral element, strain and displacement extrapolations similar to that described in [l] have to be performed. As a result, the following expressions are obtained (25a) N:^ = $ (1 + &,,U

+ +h,)(l

(25b)

+ &t-A 7

where 2, are the extrapolated values of the assumed strain field at the nodal points of the given hexahedral element, {&, ej, &} are the coordinates of the nodal point j in the space of isoparametric parameters {‘k +j, 8 of .. 1 the parallelepiped while (t,, 4,. [,} are the coordinates of the nodal point I in the space of isoparametric parameters {& ;i, [} given by @,, $9 &>’ = ]GP &,, Gtl - ’W, - X,> . 5.2. Displacement

(25c)

extrapolation

The nodal strain vahtes ?I of Eq. (25a) are stiII expressed in terms of the noda displacements of the associated parallelepiped. Therefore, a displacement extrapolation matrix is required to express the nodal displacements of the associated parallelepiped in terms of those of the original hexahedral. This allows the assumed strain field of the hexahedral element to be determined in terms of its own nodal displacements. To this end, the following relationship between deformation modes and the vector of nodal displacements of the given hexahedral element is postulated d=dC+dLS+dCB+dLB

= [C, HLS, IF*, HLB]{{PC}T, {PLy

{PC”>‘. {PLB}T}T .

(26)

The symbols used in the above equation are explained in the sequel. At this point, we note that the above equation is similar to Eq. (6), but it is written for the original hexahedral element assuming that the parameters P specifying particular modes of deformation of the hexahedral (to be described in the sequel) are identical as those for the associated parallelepiped. In what follows, this will be made more precise. The symbols appearing in Eq. (26) are defined as follows. Vector d is arranged in the following way d={u;,u;,u;

,...,

u; ,...,

u;}T,

u;={u;,u;,u;},

Wa,b)

where u, is the displacement vector at node I. Vectors dC, dLs, dCB and dLB represent extrapolated nodal displacements of the original hexahedral resulting from the constant strain, linear shear, constant bending, and linear bending modes, respectively. REMARK 5.1. The idea of the strain and displacement extrapolation is, to some degree, arbitrary. This is particularly so for the constant and linear bending modes. The main guiding requirement here is that the extrapolated fields should converge to those specified for the associated parallelepiped when the element’s shape approaches that of a parallelepiped. It was expected that fulfillment of this requirement would lead to a convergent element since the parallelepiped elements formulated as proposed herein do satisfy the patch test, and the numerical examples support this expectation. This requires that, with the mesh refinement, shapes of all the elements in the problem converge to a parallelepiped. It is hoped that this type of convergence will be rigorously proven in the future. REMARK 5.2. Even though the idea of extrapolation allows for some arbitrariness, it was quite natural to extrapolate strains using the usual trilinear shape functions. Extrapolation of the displacement field turned out to be much more intricate, and use of the same shape functions for the displacement fields related to constant and

16

Y.-I. Chen. H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

linear bending modes did not work well. Another natural idea that was pursued here was to apply to the original hexahedral the kinematic decomposition similar to that used to define particular modes of deformation of the associated parallelepiped and relate the two decompositions to one another. For constant strain and linear shear modes, this concept is relatively clear. For the remaining modes, however, this is a little more vague. Subsequent remarks are intended to add some clarification to that. 52.1. Constant strain modes Matrix C of Eq. (26) has the following c = [CT, c;,

. . . , c;, . . . ) c;,

form

C;f]’ ,

(284

C, = [Z3x3, A,,,,,en7 AILZnenyA 1,3nen,B,,,e’, &e2, where all symbols are the same as in Eqs. (Ba,b,c,d). associated parallelepiped.

BI13e3,B,,2e’&e23 It is constructed

BI,,e31 ,

in exactly

(28b)

the same way as for the

5.2.2. Linear .shear strain modes Since the linear shear strain modes need no strain modification, the part of the matrix HLs of Eq. (26) related to node I of the given hexahedral element is extrapolated from the nodal matrices fi:” of Eq. (9) by using trilinear shape functions &s

= N;fiiS

,

where N; is identified

(29) in Eq. (25b).

5.2.3. Constant bending modes To construct matrix HCB, it has been assumed that pyramids corresponding to the original hexahedral and its associated parallelepiped (defined by the nodes having identical numbers) rotate by the same angles and about parallel vectors that pass through different points. For the parallelepiped, it is the point C defined in Eq. (lb) and, for the hexahedral, six different points ‘Ci are selected (one corresponding to each constant bending mode) and defined in what follows. With these assumptions, matrix HCB of Eq. (26) is defined similarly to its equivalent for the associated parallelepiped

CB H I$; =

for I CXf- X~,)G~e,,,,,e”,

listed in Table 1 (without a ‘-‘)

for other nodes I

0,

(3Ob)

with q& representing one of the rotations listed in Eq. (4) about direction a^= $, < or [. Furthermore, to present a unified description of all those six points, the nodes of the original hexahedral are relabeled as Z, J, K, L, 0, R Q and R (see Fig. 6). The assignment of specified numbers to the above labels for each of the constant bending modes is given in Table 3. C, is located on line 13 (shown in Fig. 6) SO that the following conditions are satisfied I&; - $1

I&, - %I

Ix** -x,*1

= Ix,* -x3*1

(314

’

{XCi - Xi}. {X2* - x, *) = 0 ,

{Xq - Xi> 1{X4* - x,*1 = 0,

(31W

with

rx,*,x,*,x,*,x,*I=~~X,+x~,X,+X,,X,+x~,Xr+Xnl~ (see Fig. 6) .

(314 (31d

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998)

Fig. 6. Reference

point C, in the reference

Table 3 Variables I, J, K, L, 0, P, Q and R for each constant Node number

Constant

I 4 8 5 2 3 7 6

bending

4 3 7 8 1 2 6 5

bending

l-29

17

plane formed by points I*, 2*, 3* and 4*

mode

mode

1 2 3 4 5 6 7 8

4 8 5 1 3 7 6 2

REMARK 5.3. As shown in Fig. 6, the point Ci used in the extrapolation of constant bending modes to the given hexahedral is located on the isoparametric line 13 of a surface formed by four nodal points given in Eq. (31d). Using this definition of the reference point C, reduces the concept proposed here in the context of threedimensional analysis to a similar and rather successful approach used in [l] for the (plane) four-node quadrilateral element. The deformation pattern of the given hexahedral element corresponding to a constant bending mode is then obtained by rotating pyramid C,OPQR about the axes defined by the vectors given in Eq. (lc) passing through the point C, by the same angles that have been identified in Eq. (4).

5.2.4. Linear bending modes Two different definitions of matrix HLB are presented and implemented here. In the first version, it is assumed that the face of the original hexahedral corresponding to the face of the associated parallelepiped (having the same numbers, that is) that warps in the linear bending modes also warps without length changes of its edges. This requires that the nodal points belonging to this face displace in the direction perpendicular to the adjacent edges. The magnitude of those displacements is found by requiring that projections of those displacements on the nodal displacement vectors obtained from the extrapolation using the trilinear shape functions is equal to those extrapolated displacements. It is clear that such a procedure leads to the nodal displacements of the hexahedral which are identical to those of the parallelepiped if the original hexahedral is a parallelepiped. The mathematical expression of the above description requires unit vectors t,,, with M = 1,2,3, to be constructed at

Y.-I. Chen, U.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) 1-29

Fig. 7. Designation

used in the definition

of the unit vectors t,,

of matrices HLs and HLB.

every node I of the original hexahedral element and designated as shown in Fig. 7. The part of matrix HLB of Eq. (26) related to node I of the given hexahedral element is then given by the following formulas

J ,. LB t I#

=

NIH_?,’

#

=

5, v,

t

3

Wb)

$ee,re tf, K = 1, 2 and 3 are the vectors defined by (t,2 X t,,), (t13 X t,,) and (t,, X t,,), respectively; vectors HI,, # = 8, 71 and 6 are identified in Eq. (11). The second way of constructing matrix Hk is designed so as to create in the original hexahedral a situation similar to that existing in the linear bending modes of the associated parallelepiped. In the case of the parallelepiped, the linear bending modes create strains that along consecutive parallel comer lines alternate in sign but have the same magnitude. Even though the corresponding lines of the hexahedral element may not be parallel, the nodal displacements at the ends of those lines (still perpendicular to the warping faces of the element, as in the first approach) are selected so as to give identical projections onto the direction of the line at each of those ends and to preserve magnitude of the strains along the four lines involved. To connect deformation of the hexahedral to deformation of its associated parallelepiped, the procedure applied in the first approach at all nodal points is applied here again, but only at one of those nodes (say node 1). Again, one can see that when the hexahedral approaches a parallelepiped shape, the procedure described above recovers the corresponding linear bending mode of the associated parallelepiped. The approach leads to the following results

HLB = [HL;B, H$

where

HL;B] =

(334

19

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 1.54 (1998) l-29

(33b)

(no sum on K) . 5.2.5. Displacement extrapolation matrix Inversion of Eq. (26) provides quantitative information regarding overall deformation process of the original hexahedral element

participation

of the particular

modes in the

((PCjT, {PLSIT, {PcB}T, {PLB}T}T = [C,HLS, HCB, HLB]_‘d = l[RC]T, [RLSIT, [RCB]T, [RLBITjTd.

(34)

Identification of the parameters obtained in the above equation with those given in Eq. (12) allows the vector of nodal displacements of the associated parallelepiped to be expressed in terms of the vector of nodal displacements of the original hexahedral element. The result is ci = [e,fi”“,

tiCB, fiLB][C,HLS,HCB,

where P is called the displacement established by relating deformation element. Substitution of Eq. (35) into Eq. expressed in terms of its own nodal Z, = [A’;ii$P]d

HLB]-‘d

= ‘Pd,

(35)

extrapolation matrix. As mentioned at the beginning of this section, it is modes of the associated parallelepiped to those of the given hexahedral (25a) leads to the assumed displacements

strain field for the given hexahedral

=&d ,

(36)

where i, is the modified assumed strain matrix at node I of the original modified assumed strain field of the hexahedral element is obtained from 2 = i N,E;, = i N&d /=I I=1

= i

N,[N$,V]d

hexabedral

where @. is the domain

element.

=iid ,

of the hexahedral ‘1 - v V

and D is the matrix of elastic constant V

A

V

V

0

0

IL o-

D=(I+v)(l-2zf)

Finally,

the

(37)

I=1

where i is the modified assumed strain operator of the given hexahedral element. With the g matrix defined, the following expression provides the modified element original eight-node hexahedral element

E

element

0 0 0 l-21, 2

0 0 0

0 0 0

0

0 0 l-2v 2

0

0

0

0

I-2v ___2

0

0

0

0

0

stiffness

matrix for the

given by

-

(38b)

-

REMARK 5.4. From the preceding development one can conclude that, if nodal displacements d describe rigid body motion, no strains are generated in the elements presented here. For such a vector d, Eq. (34) must yield pc # 0 andPLS =pCB =pLB = 0. More specifically, only the part ofpC that corresponds to the rigid body motion does not vanish, the constant strain part does (cf. Eq. (2)). Since Eq. (34) constitutes a part of Eq. (35), this implies that ci = lyd = cp” also represents rigid body motion. Under these conditions Eq. (6) must yield dLs = dCB = dLB = 0 and dC = d. As a result, all strains specified by Eqs. (IS), (19), (20a) and (21a) vanish and so must the extrapolated strain field defined in Eq. (25a). Based on similar arguments, one can show that if the deformation of the original element includes only constant strain field, this field remains unchanged by the modifications described here.

20

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) 1-29

REMARK 5.5. As shown in Eq. (37), the assumed strain field defined over the entire hexahedral domain is obtained from the assumed strains at its nodal points by using the trilinear shape functions. As a result, the matrix relating the assumed strain field to the vector of element nodal displacements is a matrix of polynomials with respect to the isoparametric coordinates. Since the Jacobian of isoparametric mapping is also a polynomial, the expression defining the element stiffness matrix can be integrated analytically. Furthermore, this expression should not be ‘underintegrated’ since all the necessary strain modifications have been introduced in the B matrix. Any underintegration may, in fact, lead to an inferior element.

6. Numerical simulations and discussions Performance of the present three-dimensional assumed strain hexahedral solid elements and the elements due to the standard eight-node isoparametric formulation, Taylor et al. [3], Pian and Tong [14] and Sze and Ghali [19] is evaluated based on the analysis of several benchmark problems frequently used in the literature. The elements used in this comparison are labeled as follows: the standard eight-node isoparametric hexahedral element. H8: the Taylor’s non-conforming eight-node hexahedral element [3]. HM9: PT18p: the Pian and Tong’s assumed stress hexahedral element [14]. the Sze and Ghali’s assumed stress hexahedral element [19]. SS18/3: HMODl: the eight-node assumed strain hexahedral element described in Section 5 with the use of Eqs. (32a,b). HMOD2: the eight-node assumed strain hexahedral element described in Section 5 with the use of Eqs. (33a,b). 6.1. Cantilever

beam under pure bending load

The performance of the present three-dimensional models under pure bending conditions and the effect of mesh distortion on the accuracy is studied by using the two-element cantilever beam as shown in Fig. 8. The cantilever beam of dimension 10 X 2 X 1 is subjected to an end moment. The interelement boundary is then gradually rotated with respect to the direction of bending (a distance e on the top and bottom surfaces) to skew the mesh. For the mesh layouts stated above, Poisson’s ratio Y = 0.25 and v = 0.4999 are considered for plane stress state. In addition, analysis of this problem under the condition of plane strain is also performed by suppressing all the nodal displacements along the vector of applied bending moment. The displacements at points A and B are normalized with respect to the theoretical solution, [32]. These normalized displacements are listed in Table 4 for the proposed HMODl and HMOD2 elements. It is interesting to note that the predictions (including all displacement components and stress components) of both HMODl and HMOD2 elements are identical to the theoretical solutions (based on beam theory) even when e = 4. These results provide a good indication that the selected bending modes with the corresponding modification of the strain operator used in this work can lead to an excellent bending behavior for the eight-node hexahedral element. Thus, it is expected that the element formulated in this way is also capable of yielding high accuracy for coarse meshes and extremely low sensitivity to mesh distortion for various bending dominated problems. It has been pointed out by MacNeal [33] that all four-node membrane elements, with only two translational degrees of freedom per node, exhibit distortion sensitivity if they satisfy the patch test. Although no numerical results due to HS, HM9, PTlSp and SSlSj? elements are reported here, it is expected that similar unsatisfactory performance also exists

Fig. 8. Cantilever

beam under pure bending

load.

21

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29 Table 4 Normalized Parameter e

0 1 2 3 4

displacements

at points A and B for the cantilever

beam under pure bending

Y = 0.25 (plane stress and plane strain)

load

v = 0.4999 (plane stress and plane strain)

UA

rJ.4

WA

UB

bB

WB

UA

u,

WA

UB

b,

WI<

1.m 1.00 I .OO 1.00 1.00

1.m 1.00 1.cO 1.00 1.00

1.m 1.oo 1.00 1.OO

1.oo 1.00 1.00 1.00 1.OO

1.00 1.00 1.00 1.oo 1.00

1.00 1.00 1.oo 1.oo I .oo

1.00 1.oo 1.OO 1.00 1.00

1.00 1.00 1.00 1.00 1.00

1.00 1.oo I .oo

1.oo

I .oo 1.OO I .OO

1.co

1.oo

1.00

1.oo

1.00

I .oo

1.00

1.oo

in these eight-node hexahedral satisfy the patch test.

6.2. MaeNeal’s

elongated

elements

which have only three translational

1.oo I .oo

I .oo 1.00 1.00 1.oo

degrees of freedom per node and

beam

This elongated cantilever beam of dimension 6 X 0.2 X 0.1 is presented by MacNeal and Harder [34]. The beam is modelled using three 1 X 6 meshes consisting of rectangular prism, trapezoidal prism, and parallelepiped elements, as shown in Fig. 9. This beam is subjected to either an in-plane or an out-of-plane end shear load. This problem is used to test the locking problem which is due to the use of high aspect ratio and irregular meshes. The material properties are the Young’s modulus E = 10’ and the Poisson’s ratio v = 0.3. The normalized results are listed in Table 5. For this problem, the proposed elements exhibit excellent performance for all loading cases and mesh layouts. Severe locking phenomenon plague H8, HM9, PTlSp and SSlSp elements. Furthermore, unsatisfactory performance is noted for a parallelepiped shape of the elements.

in-plane 10ad I

/

1

I 6

/

/

I

/

(a)

I

0

/

4S’Y Fig. 9. MacNeal’s

Table 5 Normalized Element model

H8 HM9 PTlSp ss1ap HMODI HMODZ

elongated

displacements In-plane

beam: (a) Rectangular

at the free end of MacNeal’s

prism elements;

elongated

out-of-plane load

(b) trapezoidal

03

prism elements;

(c) parallelepiped

elements

beam

load

Out-of-plane

load

Rectangular

Trapezoidal

Parallelepiped

Rectangular

Trapezoidal

Parallelepiped

0.093 0.978 0.98 1 0.98 1 0.98 I 0.98 I

0.026 0.047 0.047 0.047 0.981 0.98 1

0.032 0.624 0.625 0.625 0.981 0.980

0.025 0.973 0.98 1 0.98 1 0.981 0.981

0.010 0.030 0.03 I 0.03 1 0.972 0.979

0.014 0.528 0.587 0.587 0.739 0.73 1

Y.-I. Chen, H.K. Stolarski I Cornput. Methods Appl. Mech. Engrg. I54 (1998) i-29

Fig. 10. Curved cantilever

6.3. Curved cantilever

beam.

beam

This clamped curved cantilever beam loaded with either an unit in-plane or an unit out-of-plane load at the free end, as shown in Fig. 10, is also presented by MacNeal and Harder [34]. The beam is modelled using 1 X 6 equispaced mesh. The out-of-plane load is used to test the element’s performance under a combination of bending and torsion. The material properties are Young’s modulus E = lo7 and Poisson’s ratio v = 0.25. The geometric parameters are inner radius R, = 4.12, outer radius R, = 4.32, and thickness t = 0.1. The reference tip displacements in the direction of loading for the plane stress curved cantilever beam are 0.08734 and 0.5022 for the unit in-plane load and the unit out-of-plane load, respectively. The results, normalized with respect to the analytical solutions, are summarized in Table 6. It is shown that prediction of the proposed HMODl and HMOD2 elements is superior to HS, HM9, PTlSp and SSl8/3 elements. For the case of out-of-plane load, when torsion is involved, the proposed HMODl and HMOD2 elements can still yield very satisfactory accuracy. 6.4. Twisted cantilever

beam

The beam twisted by the angle of 90” from the fixed end to the free end is also suggested by MacNeal and Harder [34]. It is used to test the element performance due to the use of warped configuration. The twisted beam is modelled by using 2 X 12 mesh and loaded with either an unit in-plane or an unit out-of-plane load as shown in Fig. 11. The material properties are Young’s modulus E = 29 X lo6 and Poisson’s ratio v = 0.22. The geometric parameters are length L = 12, depth D = 1.1, and thickness t = 0.32. The analytical tip displacements in the direction of loading are 0.005424 and 0.001754 for the unit in-plane and the unit out-of-plane loads, respectively. The normalized results are reported in Table 6. The results show that predictions of HM9, PT18/?, Table 6 Normalized

displacements

in the direction

Element model

Curved cantilever

H8 HM9 PTl@? SSl8/3 HMODl HMOD2

0.073 0.876 0.877 0.877 1.013 1.012

In-plane

load

of loading for curved and twisted cantilever Twisted cantilever

beam Out-of-plane 0.231 0.819 0.846 0.846 0.943 0.943

load

In-plane 0.208 0.995 1.001 1.001 0.999 0.999

load

beams beam Out-of-plane 0.108 0.989 0.992 0.992 1.001 1.001

load

Y.-I. Chen, H.K. Stolarski

23

I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

unit out-of-plane load

unit Fig. 11. Twisted cantilever

SSl@, HMODl and HMODZ elements poor performance. 6.5. Pinched

circular

all exhibit excellent

ill-plane

load t

beam.

accuracy.

However,

the H8 element

shows very

ring

A circular ring with rectangular cross section t X t and radius R = 1, as shown in Fig. 12, is compressed by two unit loads acting along a diameter. This problem is analyzed under the condition of plane stress state. Owing to symmetry, only one-quarter of the circular ring is modelled by using 1 X 4, 1 X 8 and 1 X 16 equispaced meshes. The material properties are Young’s modulus E = 10’ and Poisson’s ratio v = 0.3. For this problem, t = l/60, 1 / 100 and l/500 are used to test the locking problem due to the use of high aspect ratio and coarse meshes. The displacements in the direction of loading are normalized with respect to the theoretical solutions, [35], and listed in Table 7. The proposed HMODl and HMOD2 elements exhibit excellent

IP-l Fig. 12. Pinched circular

Table 7 Normalized Element model

displacements

in the direction

of loading of the pinched

circular

ring.

ring

Rlt = 60

R/t = 100

R/t = 500

Mesh layout

Mesh layout

Mesh layout

1x4

1X8

1 X 16

1x4

1X8

1 X 16

1x4

1x8

1 X 16

H8 HM9 PTlS/? SSl8/3

0.004 0.177 0.178 0.853

0.018 0.908 0.909 0.909

0.068 0.99 1 0.992 0.992

0.002 0.007 0.026 0.852

0.073 0.822 0.989 0.967

0.073 0.822 0.989 0.992

0.000 0.003 0.003 0.852

0.000 0.187 0.187 0.967

O.OCKl 0.925 0.925 0.992

HMODl HMODZ

0.952 0.952

0.988 0.988

0.997 0.997

0.952 0.952

0.988 0.988

0.997 0.997

0.952 0.952

0.988 0.988

0.997 0.997

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

24

performance for all three thicknesses. It also shows that it is insensitive to the change of R/t ratio at all. Its predictions are superior to SSlSj? element. HM9 and PTlSp elements show a severe locking problem for large R/t ratio and coarse meshes. The H8 element locks severely for all cases. 6.6. Simply supported

square plate

A simply supported square plate of length L = IO and thickness t = 0.1 is loaded by an uniform distributed loading of intensity 1. Owing to symmetry, only one-quarter of the plate is modelled by using 2 X 2, 3 X 3 and 4 X 4 evenly divided meshes as shown in Fig. 13. The material properties are Young’s modulus E = lo7 and Poisson’s ratio Y = 0.3. The vertical displacements at the centre of the plate are normalized with respect to the theoretical deflection, [36], and listed in Table 8. Results obtained by PTUP, HMODl and HMODZ elements are almost identical, and they are slightly superior to the results due to the SSlSp element. Although, eventually, the prediction of the HM9 element tends to converge to solution with the refinement of the mesh, unsatisfactory performance still exists for the mesh of 4 X 4. 6.7. Clamped circular plate loaded by central point load Displacement field due to concentrated force applied to a plate allowing for shear deformation is singular. Thus, numerical calculations concerning such cases should be viewed with proper caution. The results reported here serve only the purpose of comparison with other results obtained under identical conditions.

Y

Fig. 13. Simply supported

Table 8 Normalized Element model

displacements

at the center of simply supported Simply supported

square plate under uniformly

square plate and clamped

circular

load.

plate

Clamped

square plate

circular

plate

Number of element

Mesh layout

H8 HM9 PTlS/3 SSlSp HMODl HMOD2

distributed

2x2

3x3

4x4

3

12

48

0.009 0.101 0.990 1.038 0.990 0.990

0.020 0.466 0.997 1.017 0.997 0.997

0.035 0.737 0.999 1.011 0.999 0.999

0.006 0.116 0.517 0.807 0.667 0.670

0.020 0.579 0.869 0.899 0.907 0.917

0.072 0.926 0.983 0.983 0.974 0.972

25

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 1.~4 (1998) l-29

N-3

N-12 Fig. 14. Clamped

circular

N=48

plate under central point load.

A clamped circular plate with radius R = 5 and thickness t = 0.1 is loaded by an unit central point load. Owing to symmetry, only one-quarter of the circular plate is analyzed by using three meshes of 3, 12 and 48 elements as shown in Fig. 14. The material properties are Young’s modulus E = 1000 and Poisson’s ratio v = 0.3. The central displacements in the direction of the unit point load are normalized with respect to the theoretical solution 5.4312, [36], and are listed in Table 8. Due to the use of shear scaling factors, the SSlSp element exhibits very good behavior for irregular and coarse meshes. It is mainly related to the fact that undesirable shear strain caused by the use of distorted meshes is scaled down by using the shear scaling factors. Numerical performance of PTlS& HMODl and HMOD2 elements are similar for this problem. 6.8. Hemispherical

shell under alternating

point load

A hemispherical shell with a 18” hole at the top of the shell is subjected to two pairs of orthogonal and opposite forces (one pair in inward diameter direction, the other pair in outward diameter direction). Owing to symmetry, only one-quarter of the shell is modelled by using 2 X 2, 4 X 4, 8 X 8 and 12 X 12 meshes as illustrated in Fig. 15. The material properties are Young’s modulus E = 6.825 X 10’ and Poisson’s ratio v = 0.3. The geometric parameters are radius R = 10 and thickness t = 0.04. The displacements in the direction of loading at point A are normalized with respect to 0.094 quoted in [34]. These normalized results are listed in Table 9. The proposed HMODl element produces excellent predictions for all meshes. 98.6% of the final solution is obtained for 4 X 4 mesh. Again, use of the shear scaling factors improves the SSlSp element’s performance for coarse meshes. The proposed HMOD2 element exhibits softer behavior when coarse meshes are used, and its accuracy is very close to that of the SSltQ3 and HMODl elements when refined meshes are used. All of the above three elements are capable of yielding converged solutions with the use of 8 X 8 mesh. Although, the predictions of the proposed PSlS/? element are inferior to that of SSlSp, HMODl and HMOD2 elements, eventually, it converges to the reference solution for 12 X 12 mesh. It seems fair to note that the meshes used here consist of flat elements and that for a small number of elements (2 X 2, 4 X 4) the model is geometrically significantly different from the sphere. As a result, a significantly different behavior should be expected even if the problem of piecewise flat surface is solved exactly (that is, without any discretization). 6.9. Scordelis-L.o

roof

The Scordelis-Lo problem is loaded by gravity loading of intensity 90 unit per midsurface area. This is a membrane dominated problem to test the ability of the element to model complex states of membrane response. Owing to symmetry, only one-quarter of the shell is modelled by using 2 X 2, 4 X 4, 6 X 6 and 8 X 8 meshes as indicated in Fig. 16. The material properties are E = 4.32 X 10’ and Poisson’s ratio Y = 0. The geometric parameters are: length of the shell L = 50, thickness t = 0.25, radius R = 25 and angle C$= 40”. For the displacement in the direction of gravity loading at point A on the midsurface, two different reference solutions 0.3024 and 0.3086 (both are quoted in [34]) are commonly used in the comparison. In this work, the computed displacements at the point A are normalized with respect to the value of 0.3086 to make it comparable with the results of others. These normalized results are listed in Table 9. For this problem, PTl8/3, SSU/?, and the proposed HMODl and HMOD2 elements yield similar performance. For the mesh 2 X 2, the elements PTlS@,

26

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) 1-29

Fig. 15. Hemispherical Table 9 Normalized Element model

displacements

at point A in the direction

Hemisphericat

shell under alternating

of loading for hemispherical

shell and Scordelis-Lo Scordelis-Lo

shell

roof

roof

Mesh layout

Mesh layout

HS HM9 PTl8/3 SSlS#! HMODl HMODZ

point loads.

2x2

4x4

8X8

12 x 12

2x2

4x4

6X6

8X8

0.000 0.000 0.000 0.721 1.106 I.232

0.001 0.010 0.041 1.050 0.986 1.073

0.003 0.163 0.742 1.007 1.008 1.008

0.006 0.493 0.957 0.998 0.998 0.998

0.025 0.139 1.331 1.459 1.421 1.421

0.062 0.542 1.028 1.061 1.037 1.037

0.092 0.852 1.010 1.019 1.001 1.001

0.121 0.947 1.002 1.007 0.987 0.987

SSlSp, HMODl and HMODZ exhibit excessively soft behavior. This is probably mainly because of inadequate representation of geometry by so few flat elements. However, the SSlS/? element exhibits softest behavior due to the use of shear scaling factors. When regular meshes are used, similar softer properties given by the SSlSj? element can also be observed in the case of simply supported square plate under uniformly distributed load. The HM9 element yields improved performance with the refinement of the mesh while it locks for the coarse meshes. 6.10. Thick-walled

cylinder

A thick-walled cylinder has been chosen to provide an additional test whether the element locks for nearly incompressible material. This thick walled cylinder is subjected to unit internal pressure and is considered to be at the plane strain state. The material properties are Young’s modulus E = 1000 and various Poisson’s ratios. The geometric parameters are inner radius R, = 3 and outer radius R, = 9. Owing to axial symmetry, only a

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

Fig. 16. Scordelis-Lo

21

roof.

q&!!&_p& SYM Fig. 17. Thick-walled Table 10 Normalized

displacements

Element model

H8 HM9 PTlSp SSlSfl HMODl HMOD2

Poisson’s

in the radial direction

cylinder.

at the inner surface of the thick cylinder

ratio

v = 0.25

v = 0.49

v = 0.499

Y = 0.4999

0.986 0.991 0.991 0.991 0.969 0.969

0.845 0.986 0.986 0.986 0.966 0.966

0.359 0.986 0.986 0.986 0.966 0.966

0.053 0.986 0.986 0.986 0.966 0.966

wedge of the thick-wall cylinder with unit thickness is analyzed for the mesh layout shown in Fig. 17. In addition to the symmetric boundary conditions along the radial directions of the wedge, displacements along the axial direction are fixed on the top and bottom surfaces of the wedge to model the plane strain state. The normalized radial displacements at the inner wall are computed and given in Table 10. All elements except the H8 element are insensitive to the value of Poisson’s ratio even very close to incompressibility for this standard plane strain problem.

7. Conclusions An assumed strain formulation aiming at improved bending properties of the eight-node hexahedral element is introduced. It can be used in the analysis of spatial beam, plate, shell as well as general three-dimensional

28

Y.-I. Chen, H.K. Stolarski I Comput. Methods Appl. Mech. Engrg. 154 (1998) l-29

elasticity problems. Numerical evaluation has shown its remarkable success in the improvement of bending behavior for realistic mesh configurations. In particular, in the case of a three-dimensional straight cantilever beam under pure bending load, the theoretical solution for displacements and stresses is obtained even when severely distorted meshes are used for nearly incompressible materials. In this paper, comparison has only been made with a few well-known formulations. However, many other formulations not reported here have been reviewed. These included the same hexahedral elements with drilling degrees of freedom (which implies about twice as many total degrees of freedom) and mixed formulations. In all cases, the comparison was favorable for the method described herein. A particularly noteworthy feature of the present three-dimensional elements is their high accuracy in shell problems. Since these problems are modeled here with only one layer of elements through thickness, the number of degrees of freedom per one node of the shell is six (three degrees of freedom at the nodes located on the top and bottom shell surfaces). In terms of total number of degrees of freedom, this makes the present elements comparable to shell elements. The present formulation does, however, include through-thickness deformation and does that entirely within the context of the three-dimensional analysis (cf. Simo et al. [37]). Comparisons of the results obtained with the present elements and those obtained by somewhat equivalent four-node shell elements (not reported here) also reveals a very competitive behavior of the element described herein. From the implementation point of view, the concept involved in defining the assumed strain field of the eight-node hexahedral element is quite simple even though the explanation of the details is, admittedly, quite involved. It consists of: (1) An introduction of a parallelepiped domain which is associated with the given hexahedral element. (2) The mode decomposition and modification of the strain field in the associated parallelepiped domain. This includes identification of various modes of deformation and then proper modification of the strain field in the bending related modes. At this step, the nodal values of the assumed strain field for the associated parallelepiped are determined. (3) Development of the strain and displacement extrapolation from the associated parallelepiped to the given hexahedral element. The nodal values of the strain field of the original hexahedral element are obtained from the nodal values of the strain field of the associated parallelepiped by using standard isoparametric functions. Displacement extrapolation assumes that parameters describing deformation modes of the original hexahedral element and its associated parallelepiped are identical. As was mentioned in Remark 5.5, the element stiffness matrix proposed in this work can be evaluated analytically. This nice technical feature of the formulation is offset somewhat by the fact that two 24 X 24 matrices have to be inverted during evaluation of the element stiffness matrix to obtain projection operators and the displacement extrapolation matrix. However, in the analysis of large problems, the time spent on the formulation of the element stiffness matrices is typically only a small fraction of the total time of computations. Therefore, the increased effort in the formulation of the element stiffness matrix is not necessarily a big drawback of this approach, particularly that a payoff in the form of consistently increased accuracy is achieved. Although the proposed eight-node hexahedral elements do not pass the patch test for a general hexahedral shape, various numerical results have demonstrated its excellent convergence characteristics as long as the initially general form of the hexahedral element approaches a parallelepiped shape with the refinement of the mesh. Practical mesh refinement processes are usually based on the subdivision of the existing elements along parametric lines and satisfy this condition.

References [l] H.K. Stolarski and Y.-I. Chen, Extrapolated fields in the formulation of the assumed strain elements. Part 1: Two-dimensional problems, Comput. Methods Appl. Mech. Engrg. 123 (1995) 247-262. [2] E.L. Wilson, R.L. Taylor, W.P. Doherty and J. Ghaboussi, Incompatible displacement models, in: S.J. Fenves et al., eds., Numerical and Computer Methods in Structural Mechanics (Academic Press, New York, 1973) 43-57. [3] R.L. Taylor, P.J. Beresford and E.L. Wilson, A nonconforming element for stress analysis, Int. J. Numer. Methods Engrg. 10 (1976) 121 l-1219. [4] S. Chandra and G. Prathap, A field-consistent formulation for the eight-noded solid finite element, Comput. Struct. 33 (1989) 345-355. [5] S.M. Yunus, T.P Pawlak and R.D. Cook, Solid elements with rotational degrees of freedom: Part I-Hexahedral elements, Int. J. Numer. Methods Engrg. 31 (1991) 573-592.

Y.-I. Chen, H.K. Stolarski

I Comput. Methods Appl. Mech. Engrg.

154 (1998) 1-29

29

[6] A. Ibrahimbegovic and E.L. Wilson, Thick shell and solid finite elements with independent rotation fields, hit. J. Numer. Methods Engrg. 31 (1991) 1393-1414. [7] K.Y. Sze and A. Ghali, A hybrid brick element with rotational degrees of freedom, Comput. Mech. 12 ( 1993) 147-163. [S] D.J. Allman, A compatible triangular element including vertex rotations for plane elasticity analysis, Comput. Struct. 19 ( 1984) I-8. [9] T.H.H. Pian, Derivation of element stiffness matrices by assumed stress distributions, AIAA J. 2 (1964) 1333-1336. [lo] T.H.H. Pian and D.P. Chen, Alternative ways for formulation of hybrid stress elements, Int. J. Numer. Methods Engrg. 18 (1982) 1697-1684. [II] T.H.H. Pian and D.P. Chen, On the suppression of zero energy deformation modes, bit. J. Numer. Methods Engrg. 19 (1983) 1741-1752. [ 121 T.H.H. Pian, D.P. Chen and D. Kang, A new formulation of hybrid/mixed finite element, Comput. Struct. 16 ( 1983) 81-87. [ 131 T.H.H. Pian and K. Sumihara, Rational approach for assumed stress finite elements, Int. J. Numer. Methods Engrg. 20 (1984) 1685-1695. [ 141 T.H.H. Pian and P. Tong, Relations between incompatible displacement model and hybrid stress model, Int. J. Numer. Methods Engrg. 22 (1986) 173-181. [15] M.J. Loikkanen and B.M. Irons, An S-node brick finite element, Int. J. Numer. Methods Engrg. 20 (1984) 5233528. [16] B.M. Irons, An assumed-stress version of the Wilson S-node isoparametric brick, Computer Report CNMEICRI56, Centre for Numerical Methods in Engineering, Dept. of Civil Engineering, Univ. of Wales, Swansea, UK, May 1972. [17] Y.F. Dong and J.A. Freitas, An efficient eight-node incompatible solid element with stress interpolation, Comput. Struct. 44 (1992) 773-781. [18] K.Y. Sze, Efficient formulation of robust hybrid elements using orthogonal stress/strain interpolants and admissible matrix formulation, Int. J. Numer. Methods Engrg. (1992) I-20. [19] K.Y. Sze and A. Ghali, Hybrid hexahedral element for solids, plates, shells and beams by selective scaling, Int. J. Numer. Methods Engrg. 36 (1993) 1519-1540. [20] W.-J. Chen and Y.K. Cheung, A new approach for the hybrid element method, Int. J. Numer. Methods Engrg. 24 (1987) 1697-1709. [21] Y.K. Cheung and W.-J. Chen, Isoparametric hybrid hexahedral elements for three dimensional stress analysis, Int. J. Numer. Methods Engrg. 26 (1988) 677-693. [22] W.-J. Chen and Y.K. Cheung, Improvement of three-dimensional hybrid hexahedral elements by using orthogonal approach, Engrg. Comput. 8 (1991) 257-271. [23] W.-J. Chen and Y.K. Cheung, Three-dimensional S-node and 20.node refined hybrid isoparametric elements, Int. J. Numer. Methods Engrg. 35 ( 1992) I87 1- 1889. [24] D.P. Flanagan and T. Belytschko, ‘A uniform strain hexahedron and quadrilateral with orthogonal hourglass control, Int. J. Numer. Methods Engrg. 17 (1981) 679-706. [25] T. Belytschko, Correction of article by D.P. Flanagan and T. Belytschko, Int. J. Numer. Methods Engrg. 19 (1983) 467-468. [26] T. Belytschko, J.S.-J. Ong, W.K. Liu and J.M. Kennedy, Hourglass control in linear and nonlinear problems, Comput. Methods Appl. Mech. Engrg. 43 (1984) 251-276. [27] W.K. Liu, J.S.-J. Ong and R.A. Uras, Finite element stabilization matrices-a unification approach, Comput. Methods Appl. Mech. Engrg. 53 (1985) 13-46. [28] X.-J. Wang and T. Belytschko, An efficient flexurally superconvergent hexahedral element, Engrg. Comput. 4 (1987) 281-288. [29] W.E. Bachrach, An efficient formulation of hexahedral elements with high accuracy for bending and incompressibility, Comput. Struct. 26 (1987) 453-467. [30] H.K. Stolarski and Y.-I. Chen, Assumed strain formulation for the four-node quadrilateral with improved in-plane bending behavior, Int. J. Numer. Methods Engrg. 38 (1995) 1287-1305. [31] R.L. Taylor, J.C. Simo, O.C. Zienkiewicz and C.H. Chan, The patch test-a condition for assessing FEM convergence, Int. J. Numer. Methods Engrg. 22 (1986) 39-62. [32] S.P. Timoshenko and J.M. Goodier, Theory of Elasticity, 3rd edition (McGraw-Hill Book Co., Inc., New York, 1970). [33] R.H. MacNeal, A theorem regarding the locking of tapered four-noded membrane elements, Int. J. Numer. Methods Engrg. 24 (1987) 1793-1799. 1341 R.H. MacNeaI and R.L. Harder, A proposed standard set of problems to test finite element accuracy, Finite Elements Anal. Des. 1 (1985) 3-20. [35] S.P. Timoshenko and J.M. Gere, Theory of Elastic Stability, 2nd edition (McGraw-Hill Book Co., Inc., New York, 1961). [36] S.P. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, 2nd edition (McGraw-Hill Book Co., Inc., New York, 1959). [37] J.C. Simo, MS. Rifai and D.D. Fox, On a stress resultant geometrically exact shell model. Part IV: Variable thickness shells with through-the-thickness stretching, Comput. Methods Appl. Mech. Engrg. 8 1 (1990)91- 126.