A simple isoparametric three-node shell finite element

Compurers d S~rucfures Vol. 44. No. 6, pp. 1263-1273. Printed in Great Britain.

0045s7949192 $5.00 + 0.00 ST>1992 Pergamon Press Ltd

1992

A SIMPLE ISOPARAMETRIC THREE-NODE SHELL FINITE ELEMENT P. BOISSE,J. L. DANIELand J. C. GELIN Laboratoire de M~anique

Appliquee, University de Franche-ComtL, CNRS, Route de Gray, 25000 Besancon, France 4 June 1991)

(I&e&d

Abstract-A

three-node isoparametric shell finite element including membrane and bending effects is proposed. The element is based on the degenerated solid approach and uses an assumed strain method to avoid shear locking. An intermediate convected covariant frame is used in order to construct the modified shear strain interpolation matrix. Validation tests show that shear locking is avoided and that a reduced integration procedure can be used without any loss of accuracy which is useful for the numerical efficiency.

NOTATION the nodal value of A quantity at node i subscripts taking the value I or 2 subscripts taking the value I,2 or 3 two-component vector [E,, &Jr three-component vector [E,, E,, &]’ denotes the value of the (.) quantity in the configuration obtained in the undeformed configuration ifa = 1 ifa =2 nodal vectors of u on the whole structure and on the element under consideration, respectively

1. INTRODUffION Thin-structure analysis using shell elements is an important task of finite element research because of its large industrial application. Much effort has been and is still directed toward obtaining efficient elements. Impressive advances have been obtained over the last few years. In particular, efficient quadrilateral four-node shell elements have been constructed using a mixed interpolation of the strain components [l-3]. Triangular elements are very attractive in practical cases because they allow automatic meshing for arbitrarily shaped components. Triangular shell elements have been built using the superposition of membrane and plate bending flat elements[4-61 or by using curved shallow shell elements [7,8]. For these elements the Kirchhoff assumption has been widely used. This theory requires the normal continuity across the interelement boundaries for the formulation of conforming elements. This has been shown to be impossible[28] in the case of triangular plate or shell elements. Disadvantages of such formulations lead to two families of elements. The first is based on discrete Kirchhoff assumptions (DKT) [9-i I] or extensions of this concept to take into account transverse shear CAS 4M-G

(DST) [12, 19. In that type of formation, high order shape functions are used in Mindlin kinematic and discrete conditions to eliminate some of the variables and so adapt these interpolation functions to the number of DOF of the element under consideration. This can be seen as the application of the ‘shear constraint’ [14]. The second family is the so-called degenerated solid (or continuum base) shell element family which is the point of departure of the proposed element. A three-node bending and membrane shell element is presented taking into account out-of-plane shear. The goal of this formulation is to build a very simpie and natural isoparametric element which uses the classical linear shape functions of the three-node triangle and with a low numerical cost. This element will be well adapted to model arbitrary shaped structures with complicated meshes. Using an intermediate convected covariant frame constructed with reference to the element sides that join a given node, the transverse shear strain components are evaluated at each node taking the mean value of the components along the sides of the node under consideration. The out-of-plane shear strain components are interpolated using the linear shape function. The formulation does not refer to the special assumption of Kirchhoff-Love and is applicable to thin and moderately thick shells for modelling arbitrary geometries. It does not need any numerical adjustment factor. The element does not contain spurious zero energy nodes. A reduced integration procedure can be used without any accuracy loss by using a single set of Gauss points at the normal at the element centre of gravity. This leads to a low cost element from a computational point of view. Numerical tests are presented that show the efficiency of the element. The solutions have been reached for any thickness of the shell which demonstrates that the proposed element avoids any shear locking.

1263

P. BOISE et al.

1264

Fig. I. Element geometry definition. 2. ELEMENT

KINEMATICS

R and y are interpolated

The proposed element is included in the so-called continuum-based shell element family introduced by Ahmad et al. [15-171. Many authors have developed elements based on such an approach; especially fournode [l, 21 or higher order [I 8, 191elements which have proved to be among the most effective.

R = i:

,=I

dX

a=$; g2=ayl.

8X

g3=i3

(1)

where 5 and q are the coordinates in the plane reference element and < is the position in the thickness with < E[-1; +I]. The position of a point M is defined as the sum of its position y(M) along the pseudonormal X and the position C(M) of the associated point of the mid-surface x(M) = R(H) + y(M).

Fig. 2. Geometry.

(3)

N’(<, q)X’

and y(M)=

2.1. Geometry (Figs 1 and 2) Pseudonormal vectors X’ (1X’I = I) are defined at each node as vectors that join the top and the bottom of the shell for the node under consideration. So the continuity of the geometrical description will be enforced along the edges of the elements. A natural convected coordinate system (5, q, [) is defined associated to the covariant frame (g, , g,, g3) with

as

i N’(Lrl)Y’(5)= ,=I

(4)

i N’(Ls+, ,=I

where h’ is the thickness at node i and N’ are the classical linear shape functions of the triangular element

N’(<,v)= l-5

-v,

N2(5,v)=Cr N3(t>v)=r?. (5)

Finally x(M) = i N’(<, q)X’ + < i N’(& a); I= I ,=I

X’.

(6)

2.2. Kinematic (Fig. 3) Displacement of the point M of the shell is the sum of the displacements of the point H on the midsurface and of the displacement uR given by the pseudonormal rotation

u(M) = I(H) + u,@‘)

(2)

Fig. 3. Nodal

kinematics

(7)

A

which is interpolated

1265

three-node shell finite element

by

u(M) = i N’(L s)[gi + uR(i)l, ,=1

(8)

The uncoupled membrane and bending mechanical behaviour allows us to separate the virtual work of internal membrane-bending loads on the one hand and the virtual work of internal shear loads on the other; therefore ( 14) can be written

where I& is the displacement given by the rotation of the point Mb initially on the pseudonormal Xb at node i0 and at the same coordinate c as M

T;*+7-;-T,*,l=0

(15)

with .a=;c(x~-x;).

(9) 7% =

Following the inextensibility constraint along the pseudonormal, an orthonormal tensor R’ (det R’= + 1) transforms Xb to x’

T: =

Pzp(~Wm&I sY

[& (rl )I’Cs[-%I dV.

sY

The stress component

ujj is assumed

(10) zero.

uk=$(R’-I)X;,

where I is the second-order identity tensor. If the rotation components are small enough in the loading step under consideration, and calling 8’ the rotation vector at node i (R’ - 1)X; = 0’ A X; + O[(@)*]. An orthogonal frame (V\,, V&, Vi,) each node with Vi,, = X;, so

(11)

is defined

dV

(16) (17)

to be equal to

The in-plane strain components are related to the nodal displacements using the in-plane strain interpolation matrix

Ed = B,,u: which leads to the in-plane

at

stiffness

C, = q%rn+,

(18) matrix (19)

with 0’ = e: v;, + ejV&,

(12)

n

i.e. finally In the same way, the transverse shear strains can be related to the nodal displacements

xg(-e;v;,+e;v;,). The kinematics of the element five DOF per mid-surface node. correspond to global coordinate two rotations correspond to Vi,

3. MODIFICATION

(13)

are described with The displacements directions and the and Vi0 directions.

OF THE STIFFNESS

MATRIX

Small displacements and linear behaviour are assumed. Minimization of potential energy leads to the displacement field as a solution of

tlt~/q = 0 on the part of the boundary with prescribed displacements; f is the body load and i the surface load.

L&l= JW

(21)

which leads to a shear stiffness matrix. However, it is well known that in their standard forms such Mindlin elements lock, i.e. they produce exclusively stiff solutions when the thickness of the element is small compared to the other dimensions. In the case of the low-order element, and especially in the case of the three-node element, solutions are very inaccurate for any thickness. A number of remedies have been proposed. Among these, reduced or selective integration [20,21] is the simplest but is proved to be ineffective in the case of the three-node triangular element (see test 2). A mixed interpolation of out-of-plane shear strain components have been introduced successfully for four-node elements [l, 2, 22-241. The proposed triangle is based on such a procedure. Transverse shear components are not classically calculated from the nodal displacements but are interpolated from their own nodal values which are fixed by assumption. This procedure leads to a modified shear strain interpolation matrix

1266

P. BorssEet ai.

The method of obtaining this out-of-plane modified interpolation matrix is detailed in Sec. 4.2. The modified shear stiffness matrix is calculated as (23)

“W

Fig. 4. Nodal covariant vectors related to the element sides. and the displa~ment

field is found as a solution of (24)

4.2. InterFo~ation matrix of the out-of-crane strain components

Because modifications have to be made to the transverse shear components, the ‘transverse’ notion forbids the use of the basic Cartesian frame and has to work in a coordinate system involving the tangent plane of the element. The proposed formulation uses the convected covariant frame (g, , g,, gJ) defined in (1).

If the transverse shear strain components are calculated in the same manner as the previous in-plane strain components the obtained element is far too stiff. A review of the literature shows that a reduced integration method fails to produce accurate results for the three-node element (see test 2). A direct interpolation of these tensorial components is made from the value of those components at the

(Rs + K,)u, = F.

4. STRAIN INTERPOLATION MATRIX IN THE

nodes

COVARIANT FRAME

The strain tensor E has five si~ificant E=

components

Et& 6 g’,

(25)

where (g’,g*, g-‘) are the contravariant associated to (g,, g,, gj), i.e. g,g’ = 6.j.

vectors

(26)

In-plane components E,, , Ez2, E,, are classically derived from nodal displacements while out-ofplane or shear components E,, and E, are directly interpolated from values at nodes. 4.1. rnter~olation components

matrix

where N’ are the linear shape functions defined in (5). 4.2.1. Calculation of nodal transverse shear strain components. In order to fix the value of the nodal strain components EL, the transverse shear strains along the element sides are assumed to be constant (see Fig. 4). At each node, the following vectors are defined

of the in-plane strain

The components are classically obtained (g, , g,, g3) fixed by the nodal displacements

from

where r{ is the convected coordinate along the side (i, i + 1) with r; = 0 at node i,

Err = kg,,

- g,g,h

r{ = 1 at node i f 1

(33)

(27) and where ri is the convected coordinate along the side (i, i - 1) with

Le. EM--j:

du 1 du ~glW+dfpgIo (

>

.

(28)

1 [ +r

3 i;,

(

i; i=l

dN’ fQQgf”z

aN’ _i u > dN’ gXO~+QQ~ I) ( c

x (-e’,v:,+B;V$J

ri=l

atnodei-1.

(34)

Using the covariant frame (f, , f2, f,) the strain tensor can be expressed as

Equations (1) and (13) give

4’5

ri = 0 at node i,

1.

E = E;fk @ f’, where (f’, f2, f3) are the contravariant to (f,,f,,f,), i.e. (29)

(35) vectors related

f’t; = 6,;.

The interpolation functions W([, II) (5) give the in-plane strain interpolation matrix as

The transverse shear strain components sides are given by

&,I = Bm&.

EZ = f(b

(30)

- &,&A.

(36) along the

(37)

1267

A three-node shell finite element

Fig. 5. Side related coordinates at node 1.

r2i-r\ rj=l-5-11

Fig. 6. Side related coordinates at node 2.

These components are assumed to be constant along each side and equal to the mean value that corresponds the value at the mid-point of the side that is

(i,

At point j,j E [4,5,6], the middle of the nodes i + I), au/al is written using (13) as

+h”‘(-e:+‘V;;’ The different mean quantities ( . ), can be calculated in the case of the three-node triangular element. At each node, ri and ri are related to the < and 9 coordinates and (f&, ri,,, f$) vectors are related to the (g’,, , g:,, g’,,) vectors. In the particular case of the triangular three-node element, the g,, and g,, vectors are constant over the mid-surface of the element but g,, is not, so we shall use g,, and g,, but use g& at the nodes (see Figs 5-7). The definition (1) of g, gives (g30)mat mid-points of the sides denoted 4-5-6 (see Fig. 8).

+ eg+‘v$‘)].

(39)

On the other hand

W di+l

_

pi

ii+‘_8’

=i+l= rl -rI and

atI

0ari

=$‘-I

-

ii’

2 m

f; =

-g2

f3=-g2+gt

Fig. 7. Side related coordinates at node 3.

rj=l-5-11

rit= 5

(41)

(42)

P. BOISSEei al.

1268

&,

h’X& h2X; )

= j(

&, = ‘4( h2X& h3X; )

& Fig. 8. Vector g,, at mid-points Equation

(32), (34), (36) and (37) lead to [E~]=C’u:.

(43)

The explicit form of the C’ matrix is given in the Appendix. 4.2.2. Nodal transverse shear strain components in (g,, g,, g3). The nodal shear strain components in base (fi, f’,, f;) are related to the components in

= $ h3X;- h’X;

f

of the sides.

form of the behaviour matrix which is consequently known in the orthogonal frame (e, , e,, eX). In order to calculate the stiffness matrix a base transformation is necessary which gives the components of E in (e,,e,,e,) E = E,R,g’” Q g” = &,,e,, Q eq which leads to

(gi*gLg;)

which gives the interpolation matrix in (e,,ez, Then the stiffness matrix can be computed.

e3).

(45) 5. NUMERICAL EXPERIMENTATION

Taking

into acount

leads to

(47) In the case of the three-node

triangle,

matrices

D’ are

Some classical validation tests for the shell elements [25] are presented. The responses to the patch test for triangular plate bending element are analysed in Sec. 5.1. Example 5.2 shows the efficiency of the element based on the proposed mixed interpolation of strain components in comparison with formulations using full integration and reduced integration. The other tests show the results obtained with element in 3D shell configurations. An interesting conclusion is the effectiveness of the element when a single set of Gauss points along the normal at the triangle centre of gravity is used. 5.1. Patch test This test is defined in [26], see Fig. 9. In the case of a constant bending moment (a), the irregular mesh gives the exact transverse displacement at nodes 2 and

-1

D3 =

[ -1

1 01

(48)

.x

At each Gauss point values at nodes

&]

is interpolated

from

its

_-__ 20

My=10

40

4 35

;

3

ux=uy=oat n 1.4 c

5

[E,,] = i N’[E;,] = ,c, NiD’Ciu:,= &us, (49) ,=I which defines a modified interpolation the transverse shear strain components.

matrix

for

It-fa 1

IS

14

uy=o at node I

Y

2

-Oat

1.4 My='0

(a) Bending moment

h= 1, E= 1000, v=.3 Fz=lCl

4.3. Base transform The classical shell assumption es1 = 0 along the normaf to the shell direction e3 leads to a particuiar

node

Fig. 9. Patch test.

Fz=2

1269

A three-node shell finite element

using exact and reduced integration of shear energy and in the case of the element built on the proposed mixed interpolation. It is clear that only the proposed formulation based on the mixed integration lead to a correct solution and even with few elements. In Fig. 12, the accuracy of the solution is shown in the case of different thicknesses (mesh with ten elements per side). The solution is always achieved and no locking is observed.

Fig. 10. Beam bending.

3 and the exact stress component crYrcorresponding to MY = 1 in the element. In the case where the mesh is subjected to shear forces (b), the exact displacements are reached at nodes 2 and 3 and the exact stress or; = 1 is obtained. In the case of a twisting moment (c), the displacement solution is reached at node 2. The solution in term of stress of the problem is constant and corresponds to M,V,= 1. A very slight variation of the stress (M,. E [0.991, 1.0041) is obtained which disappears if the mesh is symmetric.

5.3. Curved beam Figures 13-I 5 show a fast convergence of the finite element model to the exact solution. An error smaller than 0.5% is obtained (considering the maximum stress as well as the maximum displacement) in the case of a mesh with ten elements per side. The computations have been made using a single set of two integration points along the normal at the centre of gravity of each element.

5.2. Beam bending (Fig. 10) In this simple test, the transverse displacement is compared to the theoretical value obtained when the Kirchhoff assumption is applied. The results are given Fig. 11 in the case of a triangular Mindlin element

12

5.4. Scordelis-Lo

roof

This more complex example (see Figs 16 and 17) requires a fine mesh in order to be able to reach a

I

1

solution

0.8 T 0,6 -

exact integration selective integration proposed triangle 0

10 20 30 40 Number of elements per Fig. 11. Displacements

50 side

for different

60

formulations.

1,05 1 ,oo ‘I-

I I

I I

I

0,95 0,90 0,85

-

I

0,80 10

. . “““I

’ ,.‘--I

100

10000

1000

- - “‘“‘I 100000

Lie Fig. 12. Displacements

for different

thicknesses.

’ . -#T-r 100

000

P. BolssE et al.

1270

Fig. 13. Curved

1,Ol

beam.

,

0998 0,97 -

__f_

Mesh A

U

hkshCD Exact

0,96 -

0,95

, 20

I 10

0 Number

solution

of

elements

Fig. 14. Maximum

per

30 side

displacement.

0,99 0,96 0,97 0,96

__t_

MeshA

V

MeshCD

1

Exact solution

0

20

10 Number

of

elements

Fig. 15. Maximum

stress.

per

30 side

1271

A three-node shell finite element

correct solution. The case of the 16 x 16 mesh has been computed using the STIFF63 shell elements of ANS~S[27](~/~=0.959~ and DKT element (f/h= 0.958)[2]. The solution obtained with the proposed element is slightly better. The element proposed is inexpensive from the computer time used because only one set of two Gauss points is used per element. 6. CONCLUSION

A three-node shell element has been developed that includes both membrane and bending strain components. It is based on the degenerated solid approach mixed with an assumed strain method in order to avoid shear locking. The formulation introduces convected vectors related to the element sides and assumes that the transverse shear strain along the element side is constant. It is simple and natural and

does not need any numerical adjustment factor. Classical linear interpolation functions of the threenode element are used for the discretization of all the degrees of freedom. The solutions of classical tests have been obtaind for any thickness of the shells and shows that the proposed element avoids any shear locking. It does not contain any spurious zero-energy modes. The formulation does not refer to the Kirchhoff assumption and the element is applicable both to thin and moderately thick shells. The numerical experimentation shows that a reduced integration of the shear energy can be made without any accuracy loss. This leads to efficient element programming. This element will be especially well adapted to the analyses that require complicated meshing, for instance when adaptive refinement is used. The method can be extended to other elements. The construction of a six-node

ed? E=3000000

Mesh A

Mesh CD Fig. 16. Scordelis-Lo roof.

1272

P.

BOISE et al.

W

MeshCD Exact Solution

0

5

15

10

Number

of elements

per

20

side

Fig. 17. Displacement at point A.

triangular shell element based on such a method is now studied.

13. REFERENCES

I.

2. 3. 4.

5. 6.

7.

T. J. R. Hughes and T. E. Tezduyar, Finite element based upon Mindiin plate theory with particular reference to the four node bilinear isoparametric elements. /. Appl. Mech. 587-596 (1981). E. N. Dvorkin and K. J. Bath, A continuum mechanics based four node shell element for general non-linear analysis. Engng Comput. I, 77-78 (1984). R. H. MacNeai, Derivation of element stiffness matrices by assumed strain distributions. J. Nuclear Engfzg Design 70, 3-12 (1982). 0. C. Zienkiewicz, C. J. Parikh and I. P. King, Arch dam analysis by a linear finite element shell solution program. Proc. of Symposium of Arches and Dams, pp. 19-22. I.C.E., London (1968). R. W. Ciough and C. P. Johnson, A finite element aproximation for the analysis of thin shells. In/. J. So/ids Struct. 4, 43-60 (1968). J. H. Argyris, P. C. Dunne, G. A. Maiejannakis and E. Scheike, A simple triangular facet shell element with application to linear and nonlinear equilibrium and elastic stability problems. Camp. Meth. Appl. Mech. Engng 10, 371-403; 11, 97-83 (1977). C. A. Brebbia, S. Sabanathan, U. Tahbiidar and H. Tottenham, Statics and dynamics of hydrostatically loaded shells by finite element method. Proc. IASS Symposium

on

Hydramech~icaIi~

Loaded

14.

15.

16. 17. 18.

19.

20.

S~ell.~,

Hawaii (1971). 8. D. 3. Dawe, High order triangular finite element for shell analysis. Int. J. Solids Struct. II, 1097-l I iO(i975). 9. J. L. Batoz, K. J. Bathe and L. 0. Ho, A study of three node triangular plate bending elements. Int. J. Numer. Merh. Engng 15, 1771-1812 (1980). 10. J. L. Batoz, An explicit formulation for an efficient triangular plate bending element. 1111.J. Numer. Meth. Engng 18, 1077-1089 (1982).

ii. M. Bernadou, Some finite element approximations of thin shell problems. In Finite Element Methodsfor Plates and She11 Structures, (Edited by T. J. R. Hughes and E. Hinton) Vol. I, pp. 62-85. Pineridge Press, Swansea (1986). 12. J. L. Batoz and P. Lardeur, A discrete shear triangular nine d.o.f. element for the analysis of thick to very

21. 22.

thin plates. Inr. J. ‘Turner, Meth. Engng 28, 533-560 (1989). 0. C. Zienkiewicz, R. L. Taylor, P. Papadopoulos and E. Onate, Plate bending elements with discrete constraints: new triangular elements. Comput. Slruct. 35, 505-522 (1990). M. A. Crisfieid, The application of shear-constraint to the generation of plate elements. In Finite Element Methods for Plates and Shell Structures, (Edited by T. J. R. Hughes and E. Hinton) Vol. 1, pp. 153-174. Pineridge Press, Swansea (1986). S. Ahmad. B. M. Irons and 0. C. Zienkiewicz, Analysis of thick and thin shell structures by curved finite elements. In<. J. Numer. Meth. Engng 2, 419-451 (1971). G. Stanley. Continuum-based she11 elements. Ph.D. thesis, Standford University (1983). T. J. R. Hughes, The Finite Element .~erhod. Linear Static and Dynamic Finite Element Analysis. PrenticeHail ( 1987). K. J. Bathe and N. Dvorkin, A formulation of general shell elements-the use of mixed interpolation of tensorial components. Int. J. Numer. Meth. Engng 22, 697-722 ( 1986). G. Stanley, K. C. Park and T. J. R. Hughes, Continuum based resultant shell elements. In Finite Element Methods .for Plates and Shell Structures, (Edited by T. J. R. Hughes and E. Hinton) Vol. 1. pp. l-45. Pineridge Press, Swansea (i 986). 0. C. Zienkiewicz, R. L. Taylor and J. M. Too, Reduced integration techniques in general analysis of plates and shells. Int. J. Numer. Meth. Engng 3,275-290 (1971). T. J. R. Hughes, M. Cohen and M. Haroun, Reduced and selective integration techniques in finite elements analysis of plates. Nut. Engng Des. 46, 203-222 (1978). R. H. MacNeai, A simple quadrilateral shell element.

Comput. Struct. 8, 175-183 (1978). 23. P. Boisse, J. L. Daniel and J. C. Gelin, A new class of

three nodes and four nodes shell elements for the finite inelastic strain analysis. Application in sheet metal forming, Proc. of the European Conf. on New advances in Camp. Struct. Mech., pp. 529-539. Giens, France (1991). 24. K. J. Bathe and N. Dvorkin. Short communication. A four node plate bending element based on Mindiin Reissner plate theory and a mixed interpolation. Int. J. Numer. Meth. Engng 21, 367-383 (1985).

A three-node shell finite element

1273

30. P. Boisse, J. L. Daniet and J. C. Gehn, Sheet metal forming using three and four node shell elements. Proc. oftheInt. Conf FE-S~maiat~on of 3f) Sheet Metal Forming Process in Automotive industries, pp. 343-357, Zurich ( 1991). 31. J. L. Batoz and G. Dhatt, Plaques et coques par elements finis; analyse lineaire et non IineEaire. Cours IPSI, Paris (1986). 32. 0. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, 4th Edn. McGraw-Hill (1989). 33. K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall (1982). 34. J. C. Simo and D. D. Fox, On a stress resultant geometrically exact shell model. Part I: Formulation and optimal parametrization. Camp. Meth. A&. Mech. Engng 72, 267-304 (1989). 35. S. Timoshenko, Theory of Plates and She/Is. McGrawHill (1959).

25. R. H. MacNeal, The evolution of low order plate and shell elements in MCS Nastran. In kite Element Methods for Plates and Shell Structures Vol. 1, pp. 85121. Pineridge Press, Swansea (1986). 26. J. Robinson, Element evaluation. A set of assessment points and standard tests., Proc. F.E.M. in the Commercial Environment 1, 217-248 (1978). 27. C. J. Desalvo and R. W. Gorman, ANSYS engineering analysis system. User’s manual. Swanson Analysis Systems (1989). 28. B. Irons and K. J. Draper, Inadequacy of nodal connections in a stiffness solution of plate bending. AIAA Jnl 3, 961 (1965). 29. M. A. Crislield, A four noded thin plate bending element using shear constraints. A modified version of Lyon’s element. Camp. Meth. Appl. Me&. Engng 38, 93-120 (1983).

APPENDIX

The out-of-plane linear strain components in the convected frame constructed with reference to the element side related to the nodal disolacements [E,*] = Cry,. (AlI I

The explicit form of Cl matrix is given here. (g&),, (g;,), , (g;,,), are the three components of g,, vectors in the global Cartesian frame [E;: E;: E;; E;: E;: E;:]‘= [C, Cz C,] u;

-(gR), 0 C, =f

C,=! 2 L

-

-f&f, (t&)X 0

(&I1, (gt;,), 0 _ -(g:&

-(&). 0

(&), (&0)X (8%)Y (&I ), (g&), (&), -(g%, -(I&, --(I&,), 0 0 _ 0

(&)X -(I&x 0 0

--GA

-(g;,), 0

-fh’V:,g,o 0

:h ‘Vt,g,, 0

$‘V:,g,o $‘V:,g,

- $‘V;,g,, - $?‘Vl,g, $i’Vl,g*, 0

-fh’V;,g,., 0

cg:, )X -$‘V:,g,, cs:, )Y -(g:o), -(g:o)x $h2V;& - g,o) 0 0 0 0 0 0 ;hrV&g,, -(&I, -4&f, (J&),

(g:o),

0

0

k:o), --(l&J), (II;”)). 0

(i& ), 0

- ( g:, )Y - (&),

$%g,, - $‘V:,(gzo- g,o) 0 0 - jh2V$,g,,

- th*V&(g.?:zo - g,,)

~h~V~~(g~-g,,)

0

k!:“), WGAg,, -g,o)

-w,

(A21

ahrV;,g, - $‘V:,g,, 0 - Sh%(g,i, -go)

0

-

1

$w&20 - T&h) _ jh $‘rl&O ;h3V;,,gzo 0 $“V:,(g,

- g,,) _

’

(A3)