Geometrically nonlinear finite element analysis of beams, frames, arches and axisymmetric shells

Computers & Structures, Vol. 7, pp. 715-135.

Pergamon Press 1977. Printed in Great Britain

GEOMETRICALLY NONLINEAR FINITE ELEMENT ANALYSIS OF BEAMS, FRAMES, ARCHES AND AXISYMMETRIC SHELLS R. D. WooDt and 0. C. ZIENKIEWICZJ: Department of Civil Engineering, University College of Swansea, Wales (Received 20 August 1976) Abstract-The geometrically nonlinear analysis of elastic inplane oriented bodies, e.g. beams, frames and arches, is presented in a total Lagrangian co-ordinate system. By adopting a continuum approach, employing a paralinear isoparametric element, the formulation is applicable to structures consisting of straight or curved members. Displacements and rotations are unrestricted in magnitude. The nonlinear equilibrium equations are solved using the Newton-Rauhson method for which a number of examples are given. The derivations are extended to include axisymmetric structures. INTRODUCl’ION

Geometrically Nonlinear (GNL) analysis may be formulated in either a total Lagrangian or Eulerian coordinate system, the former in terms of the initial position; the latter in terms of the final deformed state. Computationally an Eulerian formulation is strictly an updated Lagrangian approach[l] where the initial position becomes the current equilibrium state prior to some incremental change. The total Lagrangian approach offers advantages since the initial configuration remains constant which simplifies formulation and computation. Furthermore, anisotropy presents no problem[2]. The GNL analysis of oriented bodies, i.e. beams, plates or shells, whether straight or curved, may be formulated using theories such as those given by Von Karman for plates and Marguerre for shallow shells[383. The strain displacement relations employed in these total Lagrangian methods are restricted to small rotations. Also, in a finite element context difficulties are encountered in representing rigid body motions of oriented bodies undergoing finite displacements[9, lo]. In a total Lagrangian formulation small rotation and rigid body restrictions may be overcome, in the case of straight beams, by using Budiansky’s nonlinear shell theory[ll], but for curved beams or shells finite displacement rigid body modes only admit small rotations[9, lo]. An updated Lagrangian formulation[l2-141 seeks to overcome the above difficulties by restricting displacements, from a current initial position, to satisfy the small rotation criteria and minimise the rigid body motion error. Alternatively, the advantages of a total Lagrangian formulation may be gained by considering the oriented body as a continuum [l, 2,8] which immediately removes any restrictions on the magnitude of displacements for straight and curved elements. The resulting formulations are uncomplicated and generally applicable but suffer the penalty of large numbers of degrees of freedom. The purpose of this paper is to present such a con-

tinuum approach for the GNL analysis of oriented two dimensional structures using a modified isoparametric element which omits midside nodes in the “thickness” direction[8], thereby reducing the number of degrees of freedom, These, straight, curved or varying thickness, elements which incorporate reduced integration and modified elastic constitutive relations permit GNL analysis for any magnitude of displacement. The solution to the governing nonlinear equations is achieved using the Newton-Raphson method combined with load or displacement incrementation. A modified Newton-Raphson method is also employed whereby the tangent stiffness matrix is inverted on the second iteration of each increment thereafter remaining constant. Solutions are presented for the column elastica, toggle frame, deep arches and an axisymmetric shallow shell. FINlTEDEFORMATION THEORY Some basic concepts are reviewed, these being essential for the formulation of the Lagrangian Equilibrium equations and their finite element solution. The reader should consult references [2,8,15-171 for complete details. A bar above an item (-) denotes reference to the deformed position of the body; the bar omitted signifies reference to the initial undeformed position. Co-ordinates are with respect to fixed Cartesian axes. A material particle p in an initial undeformed position is identified by the Lagrangian co-ordinates, x = 1x,Y, z1*

(1)

the same particle is identified in the final deformed position p’ by the Eulerian co-ordinates: z = [f y, 51’.

(2)

If the displacements u of the particle p in deforming to p’ are given as a function of x or f by, u(x) = u(Z) = [u, u, WI’

(3)

then the final co-ordinates are given by, tlecturer. SProfessorand Head.

%=x+o. 125

(4)

R. D. WOOD and 0. C. ZIENKIEWKZ

726

de = de0t A#de.

Now the relation between % and x enables the deformation Jacobian matrix J to be found as,

(17)

In the virtual work equation Almansi strain is associated with the Cauchy stress, 6, referring to stresses in the deformed body and is given as, cr = [&,,&, & f,,, F=,,FX,,]‘.

(5)

Displacement gradient vectors and matrices may be defined as functions of x or f as,

(18)

These stresses may be redefined in terms of the initial undeformed Lagrangian co-ordinates of p as the PiolaKirchhoff stress vector, u, which is given by eqn (18) with the bar removed. NONLINEAREQUILIBRIUMEQUATIONS

au

aw -$ z, (6)

The equilibrium equation for the deformed body is established from the virtual work equation given in terms of the Eulerian co-ordinate system, by equating virtual internal and external work as,[2,8,17] d&‘edri= I6

and

,5du=qdg+ duTpdA I A‘ Ia

(19)

(7) where the external work is due to du acting on (i) the surfac_etractions, extending over the deformed surface area A and, given by

where typically aW ’

(8)

a.2 1

ey 0

ez 0

8x

ex 00 ez 0 ey 0 8x 82 0ey1

.

(9)

p = MX> FY,&I’

(20)

and (ii) the body forces per unit mass, acting within the deformed volume 6, given by q = 14x,qy’8lT

(21)

Similarly, 5, 8 and A, can be obtained by replacing x by %. In the Eulerian co-ordinates, Z, the Almansi strain is given as,

p being the density of the deformed body. Note that dg, is the variation of the linear part of the Almansi strain. In t_hevirtual work eqn (19) the limits of integration, 6 and A, are unknown since they are a function of the g = t&7r,,63%Y,%x, %I’ (10) deformed co-ordinates f. A Langrangian equilibrium equation is obtained by transforming the limits of inwhich in terms of the linear COand nonlinear & dis- tegration to the initial undeformed known values of u placement gradients is, and A giving[2,8,17] c=$t&

~d~Tud”=~pdu7qdvt~du’podA

(20)

where #g L =-1,,2 .9.

where

(12)

Similarly, the Green’s strain l is given in terms of the undeformed Lagrangian co-ordinates as E=EotEL

(13)

=l;ie

(14)

where Q

= 2 O* If a virtual displacement du (= dii) is given by du = [du, do, dw]r

(15)

then the corresponding variation in the Almansi strain is[171, da=d$-;i,d# and the variation in Green’s strain is,1171

(16)

Consequently, in a Lagrangian description the internal work is a function of the variation of the linear and nonlinear parts of Green’s strain. Equation (22) is now approximated by representing the continuum by a finite element discretisation. FMTE ELEMENT FORMULATION In the following total Lagrangian formulation the equilibrium equation is considered for a single element, the complete structure equilibrium equation being assumed established under conditions of nodal equilibrium and compatibility by the displacement method[l8]. The displacement u within an element is given as a function of n discrete nodal displacements 6 as, u=NS

(23)

where S = [S,, &, . . S,lr and N is the shape function

721

Analysis of beams, frames, arches and axisymmetric shells array whose co-efficients are functions of the initial position x within the element. From eqn (23) the virtual displacements du are expressed in terms of the nodal virtual displacements ds as, du=Nd&

(24)

The displacement gradients 8 of eqn (7) are given in terms of the nodal displacements 6 by the linear relation, 6=GS

(25)

where the coefficients of gradient matrix G are Cartesian derivatives of the shape functions contained in N. Due to the linearity of eqn (25) the variation dt9 of 8 may be written, d@=Gd&

(26)

The linear Green’s strain vector l0 of eqn (13) is written in terms of S as, lo =

Boa

(27)

where the constant strain matrix B. consists of Cartesian derivatives of N. The variation de0 of l0 is

which is a vector of equivalent nodal forces due to the body forces q and tractions p” and may for convenience include concentrated nodal forces. Then the nonlinear equilibrium equations become: BTudv-R=O.

This is identical with that obtained in the infinitesmal displacement case with the crucial exception that B is a linear function of nodal displacements S. Equation (35) has the dual role of representing either the element, or in its assembled form, the total equilibrium equation. In the present context the Piola-Kirchhoff stress u is considered as a linear function of Green’s strain and is given, together with its variation by, a=De

(36)

da=Dde

(37)

where D is a symmetric matrix of constitutive coefficients. From eqns (29) and (36) it is evident that the stress, a, is a nonlinear function of S. In what follows the loading vector R is conservative. SOLUTION TO NONLINEAR EQUILIBRIUM

de0 = Bodg.

(28)

The components of Green’s strain vector E in terms of the nodal displacements 6 may now be written from eqns (9), (13), (14) (25) and (27) as, E= Bo+;BL 6 [ I

(35)

EQUATIONS

The basis of the solution algorithms for the assembled nonlinear eqns (35) is the Newton-Raphson method[l9] involving a series of solutions to linear incremental equilibrium equations. Assuming an initial estimate of the total displacements is known as Si for which @((so# 0 then the value of the function $ for an increment A$ in iJ is given, by the Taylor’s series expansion of & about $ ignoring third and succeeding terms, as,

where the linear strain matrix $(Si + A&) = +(Si) + KrASi +a . . B,(6) = A,G

and since A0 is a function of 8, l is a nonlinear function of s. From eqns (17), (23) and (28), the variation de is given by the equation de=Bds

(31)

where the strain matrix B(S) = B”+ A,G.

(32)

The total Lagrangian virtual work expression (22) is approximated by the finite element idealization by the substitution of eqns (24) and (31), giving d@

I”

B=a dv = ds=

I”

N*pq du + diV

NTpodA. IA

Since the virtual nodal displacements ds are arbitrary, the element nonlinear equilibrium equations are established from eqn (33) as, =lNTpqdv+lNTpodA.

where Kr is called the tangent stiffness matrix and is given by

where k,l=l,n. From eqn (38) the linearized approximation to the relation between the residual force vector I,%and the resulting increment in nodal displacements A$ necessary to achieve equilibrium is, A$ = - Kr-‘&

(34)

Let the right hand side of eqn (34) be represented by R

(40)

from which a new approximation to the total displacements is obtained as. $+I = Si + A&

(33)

lB%d.

(38)

(30)

(41)

In order to find a complete equilibrium path, R is applied as a series of incremental loads. Iterations continue within a load increment until + satisfies a given convergence criteria. Since the Newton-Raphson (N-R) process requires repeated calculation and invertion of K, a modified N-R method is also used such that Kr is calculated once only on the second iteration of each load increment. Unless stated otherwise, a complete N-R method is used.

128

R.

D. WCODand 0. C.

From eqn (39) an element tangent stiffness matrix Kr evaluated at, S = Si is K

(42)

T

ZIENKIEWICZ TANGENTSTIFFNESMATRIXFOR PARALINEARRLmfRNT

For the two dimensional case, the co-ordinates and displacements reduce to x, y and u, D respectively. Green’s strain vector becomes

using eqn (35) this becomes

l

Kr = lB’[$]

dv +l(LB)r

dv

(43)

= [e,, 4. hylT

corresponding to the two dimensional Piola-Kirchhoff vector

where

u = [ox,o;, g. i=1,6

and j=l,n

(49)

(50)

The formulation given below is easily extended to include axisymmetric structures (see Appendix I).

and L is a differential operation defined by, L,,=u~$;

I

l=l,n

and

k=l,6.

Substituting eqn (31) into (37) gives,

!!!

[ as

1

=DB

.

Recalling that the strain matrix B(S) contains a constant part B0 and a linear part BL, where BL = AeG, the last term of eqn (43) reduces to, jjLB)’ dv = 1 (LA,G)r dv.

(45)

Noting that a typical component of a term in LAB, say (a*u/axas,) is in fact a co-efficient of the gradient matrix G of eqn (26), Kr may be expressed as, Kr =K+K,

(46)

where K=

I”

B=DBdu;

K=

I”

GrSGdu

Shape functions The paralinear isoparametric element is intended for the analysis of oriented structures where the geometry is such that the thickness is small compared to other dimensions. The characteristics of the element are defined by the geometry and displacement interpolation functions which are linear in the thickness direction and parabolic in longitudinal directions (see Fig. 1). Consequently, the element allows for shear strain energy since normals to a mid-surface before deformation remain straight, but not necessarily normal to the midsurface after deformation. By changing the subscript 6 to 4 and using shape functions shown in Fig. 1, the derivation is applicable to the linear element. The geometry is interpolated as x = [N,, N2, . . . &lx

= NX’ (51)

y=NY’

where typically x’ = [X,, X2,. . .X,1’ gives the nodal x co-ordinates. For an isoparametric element, the displacements u and o are similarly defined in terms of nodal displacements U and V as,

(47)

u=NU

and

(52)

v=NV S=[;;

2;

and I3 is a 3 x 3 unit matrix.

zj

(48)

where typically, U = [U,, U,, . . U,]’ and for reference 8 of eqn (23) is S=[U,V,,U*V* )... T&V,]?

~~~~ gi = 0, Ti=+-l Linear

2

+ 0

Paralinear

Fig. 1. Two-dimensional linear and paralinear isoparametric elements-geometry and shape functions.

(53)

729

Analysis of beams, frames, arches and axisymmetricshells The deformed position 57,y of a particle originally at x, y is obtained from eqns (4), (51) and (52) as,

From eqn (32)

alv, o 1

Bi=Boi+ABGi

i=NX'+NU=Np y=NY’+NV=N+

(54)

where

where typically !@= [xl, x2, . . . .%,]‘. Consequently, for finite deformation the isoparametric element deforms into another isoparametric element defined by the same shape functions N, this allows easy calculation of the deformation Jacobian matrix, J.

Boi =

!

Co-ordinate and deformation matrices The relationship between the derivatives of N with respect to the curvilinear 5 and n and those with respect to x and y necessary in the formation of B and G matrices is fumish~ by the co-ordinate Jacobian matrix J, where,

&‘N;r dy

ax 0

i=1,6

2% ay

(63)

ax

64

a~ a~ au a~ ---ay ay ax ax

1

(65)

Gi = J, =

(62)

(55) Noting from eqn (4) that

If the curvilinear derivatives of the shape functions are written as.

allows the submatrix Bj to be concisely written as then the co-ordinate Jacobian J, may be written, (57) The Cartesian derivative vectors G, and G, of the shape functions N are found from eqns (55) and (56), using the chain rule of differentiation as,

B1 =

i

a2 aN. _,I. ax ax ’ a.f aNi -.ay ay t az aN,+g ahi ay ax axSay

ag aNi _.ax ax _.-

1 (67)

Anisotropic constitutive relations. The co-eilicients of the elasticity matrix D given in eqn (36) are considered for the particular case of the linear or paralinear ele158) ments simulating oriented body behaviour. A prismatic beam under pure bending exhibits a parabolic variation of displacement across the depth[21] but the linear shape functions in the n direction constrain this displacement where derivatives of 5, n with respect to x, y are ob- to be constant. Consequently under pure bending the tained from J,-‘. The deformation Jacobian of eqn (5) is strain in the n direction is zero and bending occurs under calculated from, conditions of plane strain and not plane stress as is usually assumed in flexure theory. For a general loading a “spurious” local transverse stress cri (see Fig. 1) is introduced, which in a predominately bending situation increases the bending strain energy resulting in overstiff behaviour. Consequently, the isotropic elasticity matrices norStrain and gradient matrices. The gradient matrix of mally used for plane stress and strain ~haviour have to be modified. This must be accomplished in the local eqn (25) rewritten in terms of nodal submatrices is, Cartesian axes, x’, y’, (see Fig. l), where x’ is concurrent . . . G$. G = [G,,C&r (60) with the base vector in the t direction. For an isotropic material the local Green’s strains E’ are related to the Similarly the strain matrix, of eqn (31) is, expressed as, local Piola-Kirchhoff stresses u’ by the elasticity matrix D as, e’ = D-l,,‘. B = D,, B2, . . . B$‘. (&s) (61) CAS Vol. I, No. b-D

R. D. WOOD and 0. C.

730

The effect of the spurious stress crk on the longitudinal strain l: is eliminated by zeroing the co-efficient in the D-’ matrix connecting these two quantities. Inversion of the modified matrix yields a symmetric anisotropic local elasticity matrix D which ‘relates the stress vector u’ to the strain vector E’ by, @’= riel.

(69)

All that remains is to transform the local elasticity matrix 8 into the global equivalents D by the standard relationship [ 181 D = TilTT (70) where, if s = sin a and c = cos a

T= [;;

fr,

cj;:].

(71)

For the plane stress case the modified elasticity matrix is [8]

1

ZIENKIEWICZ

The above integrals are evaluated using numerical integration. A 2 x 2 Gauss quadrature being used for both the paralinear and linear elements, implying a “reduced integration”[8,20] for the paralinear element. Above a length to thickness ratio of 25 the stiffness matrix of the paralinear element becomes ill-conditioned. This difficulty may be overcome by the introduction of a matrix re-conditioning scheme,[8] in which top (T)= 1) nodal displacements are expressed in terms of bottom (T = - 1) displacements and the difference between the two. The unknown displacements now become the bottom (q = - 1) values and the differences. Examples Example l-Column elastic. Here, the problem is to determine the shape of a column under vertical end loads higher than the elastic buckling load. The column is divided into five paralinear elements and is loaded eccentrically (see Fig. 2) to initiate sway deformation. Two cases are considered, a constant section and a linearly varying section over the lower two elements. The results for the former case show negligible discrepancy when compared to analytical solutions [22,23] which assume inextensibility and no shear deformation. Example 2-Williams toggle frame. Williams [24] solved this problem analytically and compared the results to experimental observations. The finite element, results using a modified N-R method, are shown in Fig. 3, where good agreement is obtained. Due to the presence of “snap through’ behaviour displacements rather than loads were prescribed at the apex. Example 3-Two-hinged deep arch. This problem has been extensively treated by Huddleston] and by Da Deppo and Schmit[26]. Both references show that the lowest buckling load occurs in an asymmetric mode characterised by a bifurcation in the central load deflection curve. The asymmetric deformation is usually initiated by introducing a small perturbing load[27,28]. Alternatively, in the vicinity of a bifurcation point a small prescribed displacement proportional to the first buckling mode can be introduced[lO]. Here, the asymmetric deformation is induced by imposing a small perturbation to the radius of the arch, see Fig. 4, as,

(72)

where v is Poisson’s ratio. For plain strain [8],

(73)

These elasticity matrices now give stress strain relationships corresponding to beam, plate or shell behaviour. For example, from eqn (72) a: = EE: which is the relationship used in simple beam theory, from eqn (73) o: = [E/(1 - ~‘)]a: which is the expression obtained when considering the bending of a strip of an infinitely long plate. Element tangent stiffness matrix. The submatrix K, of the element stiffness matrix Kr may be defined as B,TDB,*t . det J, . dt dq t

GTSGj .t.det J, . d[ dq

(74)

where (75) and t = dimension normal to x-y plane. The later term is eqn (74), is more economically expressed as (76)

For the asymmetric buckling the maximum radius is only l/1000 times the thickness greater than the original radius, for the symmetric buckling k = 0 and only half the arch need be analysed. Plots of central point load vs central deflection and horizontal reactions (see Figs. 5 and 6) show good agreement with Ref. (25). The asymmetric buckling load of 13EI/R2 agrees with that interpolated from numerical results presented by Da Deppot261, as does the finite element central deflection of 0.106R compared to 0.108R. The finite element symmetric buckling load of 15.3EI/R* compares well with the value of 152EI/R* presented by Huddleston[25]. A NewtonRaphson solution was used with prescribed central displacement. Example 4-Clamped-hinged deep arch. This problem of arch instability after large asymmetrical prebuckling deflections has been investigated by Da Deppo and Schmit [29] under the assumption of inextensibility and no shear deformation. Crown load deflection curves presented in Fig. 8 show good agreement with the analytical

Analysis

of beams, frames, arches and axisymmetric

shells

Eccentric Loading

-O-

Finite

Analytical

+

B=

Element

1.10

Thickness E= 12

‘2’*22)

=1

Per -0.247

Y 1O-3 (B=l)

Fig. 2. Column elastica loaded eccentrically.

H lb. 0

7OT

loo

290

LOO

300

5qo

L&:; I+

12.9L3”

E= 10.3 x lo6 Lb/in2

cl

v.0.0

X-Section 0.753” -4 5

Paralinear

Williams (24 1

0.2L3”

Elements

over

Half

Experimental Analytical Finite

Frame

+ -

Element

o

V in

Fig. 3. Williams toggle: geometry and load-deflection curves.

results even for the very large deformations occuring in the region of the buckling load. The finite element buckling load of 9.24EI/R2 gives a 3% error when compared to the analytical value of 8.97EI/R*. Example S-Spherical cap. This example studies the large deflection behaviour of an axisymmetric shell subjected to a central point load and ring loads. With the

central point load, this shell is used as a “standard” example by a number of investigations [l, 17,301 who compare with analytical results presented by Mescall[31]. A prescribed displacement solution was used, the modified N-R method being adequate for r/R = 0.0 and 0.25 but a complete (N-R) approach being necessary for r/R = 0.42. Good agreement is obtained for the point load

132

R. D. WOOD and 0. C. ZIENKIEWICZ

Detail

’

Radius Area

of

Hinge

= 100 in. =

1 in?

E = lOxlO6

lb/in?

10 Paralinear

Elements

+ 2

Linear

Elements

Fig. 4. Two-hinged deep arch geometry. 16 r

Ref (25) Axisymmetric i

1

, Finit:

0.00

,

EIe;T

,

Symmetric

0.00

0.16

,

Buckling,

0.2L

k = 0.001

,

Buckling

,

x t

,

0.32

0.40

Fig. 5. Two-hinged deep arch: central load-deflection. 16

I4

I-

12

10

-

Ref

(25)

Notation

6

/

,

I

1

I

I

2

L

6

6

10

I

12

,

1L

I

16

I

16

B2 El

Fig. 6. Two-hinged deep arch: central point load-horizontal

reaction.

as

,

20

Fig. 5

Analysis of beams, frames, arches and axisymmetric shells

733

0.06598 in

E= IO x 10’ Y eo.3 6

Paralinear

Lb/in’ Elements

Fig. 9. Spherical cap.

R=lOO 16

, ‘3i 215’.t =1 0

Parolineor

Elements

EI I lo6 + 1 Linear

Element

Fig. 7. Clamped-hinged arch: geometry.

-0-----Analytical

-O-

Finite

-----

Finite

Element I30

)

Element

Analytical

(29

1 Fig. 10. Spherical cap: load-deflection curves for various ring loads.

I

0.2

I

0.4

I

0.6

0.6

1

,

I.0

I.2

uv Ti’Ti

Fig. 8. Clamped-hinged arch: load deflection.

case. As the radius of the ring load increases the emergence of “snap through” behaviour is more evident accompanied by a distinctive buckling load. CONCLUSIONS

A two dimensional finite element solution to geometrically nonlinear problems is presented. By introducing a modified isoparametric element, suitable for oriented structures, deformations of any magnitude are allowed. This combines the advantages of a total Lagrangian description with the concise isoparametric formulation. Numerical solutions have been compared to analytical results for a number of different types of structure and

excellent correlation has been obtained. The method is general and is not restricted to any particular type of inplane or axisymmetric oriented structure. Also the formulation presented can readily be extended[8] to the three dimensional case allowing the GNL analysis of shells. Furthermore, by introducing a stringer element [32] the possibility exists for the GNL analysis of combined beam and cable structures and for structures consisting of laminated materials. REFERFMXS K. J. Bath and M. Ozdemir, Elastic-plastic large deformation static and dynamic analysis. Compt. Struct. 6(2), 81-92 (Apr. 1976). 0. C. Zienkiewicz and G. C. Nayak, A general approach to problems of plasticity and large deformation using isoparametric elements. Proc. Conj. Matrix Meth. Struct. Mech. Wright-Patterson AFB, Ohio, (Wl), (pub. 1973). C. Brebbia and J. Conner, Geometrically nonlinear finite element analysis. J. Engng Mech. Do. AXE %(EM2), 463483 (Apr. 1969). H. T. Y. Yang, Flexible plate finite element on elastic foundation. J. Struct. LG. ADCE %(STlO), 2083-2101(Oct. 1970).

134

R. D. WOOD and 0. C. ZIENKIEW~CZ

5. T. Y. Yang, A finite element procedure for large deflection analysis of plates with initial imperfections. AIAA J. 9(8), 1468-1473 (Aug. 1971). 6. P. G. Bergan and R. W. Clough, Large deflection analysis of plates and shallow shells using the finite element method. Int. J. Num. M&h. Engng 5, 543-556 (1973). 7. T. M. Roberts and D. G. Ashwell, The use of finite element mid-increment stiffness matrices in the post-buckling analysis of imperfect structures. Int. J. Solids Struct. 7.805-823 (1971). 8. R. D. Wood, The application of finite element methods to geometrically nonlinear finite element analysis. Ph.D. Thesis, University of Wales, Swansea C/Ph/20/73 (1973). 9. D. J. Dawe, Finite deflection analysis of shallow circles by discrete element method. Int. J. Num. Meth. Engng 3(4) 529-552 (1971). 10. T. Matsui and 0. Matsuoka, A new finite element scheme for instability analysis of thin shells. Int. .J. Num. Meth. Engng 10(l) 145-170 (1976). 11. M. Epstein and D. W. Murray, Large deformation in-plane analysis of elastic beams. Compt. Struct. 6(l), l-9 (1976). 12. D. W. Murray and E. L. Wilson, Finite element large deflection analysis of plates. J. Engng Mech. D&I. AXE 95(EMl), 143-165 (1%9). 13. T. Y. Yang, Matrix displacement solution to elastic problems of beams and frames. Int. J. Solids Struct. 9, 829-842 (1973). 14. P. Sharifi and E. P. Popov, Nonlinear buckling of sandwich arches. J. Engng Me& Div. ASCE 97(EM5) 1397-1412 (1971). 15. Y. V. Novozhilov, Theory of Elasticity. Pergamon Press, Oxford (1961). 16. Y. C. Fung, Foundation of Solid Mechanics. Prentice Hall, Englewood Cliff, New Jersey (1965). 17. G. C. Nayak, Plasticity and large deformation problems by the finite element method. Ph.D. Thesis, University of Wales, Swansea C/Ph/l5/71 (1971). 18. 0. C. Zienkiewicz, The finite Element in Engineering Science, 2nd Edn. McGraw-Hill, London (1971). 19. A. S. Householder, Principals of Numerical Analysis. McGraw-Hill, New York (1953). 20. J. M. Too, Two dimension, plate, shell and finite prism isoparametric elements and their applications. Ph.D. Thesis, University of Wales, Swansea C/Ph/14/71 (1971). 21. S. Timoshenko and J. N. Goodier, Theory of Elasticity, 2nd Edn, p. 253. McGraw-Hill, New York (1951). 22. R. Frisch-Fay, Flexible Bars. Butterworths, London (1962). 23. S. Timoshenko and J. M. Gere, Theory of Elastic Stability, 2nd Edn. McGraw-Hill, New York (l%l). 24. F. W. Williams, An approach to the nonlinear behaviour of the members of a rigid jointed plane framework with finite deflection. Quart. J. Mech. Appl. Maths. 17(4), 451-469 (1964). 25. J. V. Huddleston, Finite deflection and snap through of high circular arches. J. Appf. Mech. 35, 763-769 (Dec. 1968). 26. D. A. Da Deppo and R. Schmidt, Side sway buckling of deep circular arches under a concentrated load. J. Appf. Me&. 36, 325-327 (June 1969). 27. S. L. Lee, F. S. Manuel and E. C. Rossow, Large deflections and stability of elastic frames. J. Engng Mech. Div. ASCE 94(EM2), 521-547 (Apr. 1968). 28. J. L. Batoz, A. Chattopadhyay and G. Dhatt, Finite element large deflection analysis of shallow shells. Int. J. Num. Meth. Engng 10, 35-38 (1976). 29. D. A. Da Deppo and R. Schmidt, Instability of clamped-hinged circular arches subjected to a point load. Trans. ASME. 894896 (Dec. 1975). 30. W. E. Haisler, J. A. Stricklin and F. J. Stebbins, Development and evaluation of solution procedures for geometrically nonlinear analysis. AMA J. 10, 264-272 (1972). 31. J. F. Mescall, Large deflections of spherical shells under concentrated loads. J. Appf. Me&s. 32, 936-938 (1965). 32. B. Schrefler and R. D. Wood, Geometrically Nonlinear Analysis-a formulation for isoparametric curved or straight line elements in a Lagrangian co-ordinate system. Research report C/R/280/76. Department of Civil Engineering, University of Wales, Swansea (1976).

APPENDIXI Axisymmetric case. The formulations given previously for the component matrices of the tangent stiffness matrix are extended to the axisymmetric case as follows. The x and y co-ordinates become r and z and the corresponding displacements u and w. The circumferental Green’s strain is given by[8]

-11

1 r2 ee=2 [(I ;

But f = r + u which allows ee to be separated nonlinear parts as, u Ee=-fr Hence the axisymmetric

into linear and

1 Ii2 2r0

Green’s strain vector becomes E = k,, E,, yrz, EslT

corresponding

to the Piola-Kiichhoff

stress vector

u = [(T,,uz, 7,2, usIT The matrices of eqn (62) are rewritten as

(83)

0

au

ar

aw

au

0

0

(84)

azar

(85)

From eqn (4) and noting that

f= ,+;

(86)

the submatrix B, is rewritten as

ai aN,

ar aN

_.2

ar

Bi =

-.-

ar ’

ar

ar aiq -.a2 a2 ’ a? aNi

ar

az aru, -.a2

aP aN,

az

ar aN,

ai

aN

(87)

az’ar+ar’x, az’ar+ar’z 0

For the axisymmetric

case the modified elasticity matrix is (8)

(88)

Analysis of beams, frames, arches and axisymmetric shells

135

from which the meridian stress

cm

E d==g-qjd

where

which is that given by thin shell theory. The element tangent stiffness matrix is calculated, as before, from eqn (74) where S becomes

1 _I

u.12 T,zL 0 s = TJZ a,I* 0 0 CT” i 0 and

(89)

xii=I_‘Il_:[g

$I[:

::][:I

and Aij= c,~?,

t.det.I;d*dq