Efficient large deflection analysis of rectangular plates with general transverse form of displacement

PERGAMON

Computers and Structures 71 (1999) 333±352

Ecient large de¯ection analysis of rectangular plates with general transverse form of displacement G.H. Little * School of Civil Engineering, The University of Birmingham, Edgbaston, Birmingham B15 2TT, UK Received 25 May 1997; accepted 4 August 1998

Abstract Thin, elastic, rectangular plates are analysed. Each plate has symmetry about the transverse centre line, but the transverse form of displacement can take any form, depending on the degree of support (if any) along each longitudinal edge, and on the manner of loading. Classical plate theory is used (Love±Kirchho€ hypothesis), and large de¯ections compared to the thickness are accommodated. The lateral displacement is expressed as a series, each term of which involves the product of a sinusoidal term in the longitudinal direction, and a polynomial term in the transverse direction. The von Karman/Marguerre compatibility equation is satis®ed, and equilibrium found at each load stage by directly minimising the total potential energy. The method is presented, and is then used to analyse a wide variety of di€erent plate types. The method is shown to be very ecient, and accurate, by comparison with previous work. The results clearly show the subtle variations in transverse form of displacement for the various examples considered. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Rectangular plates; Elastic behaviour; Classical plate theory; Large de¯ections; Series solution; Buckling; Imperfections

1. Introduction The problem considered is the behaviour of elastic rectangular plates undergoing moderately large de¯ections, i.e. de¯ections that are of the same order as the plate thickness, but small compared to the in-plane dimensions. Attention is restricted to thin plates, for which classical plate theory is appropriate (Love± Kircho€ hypothesis). In previous work, the author developed an ecient method for analysing such plates [1]. That method represents the lateral displacement w as a double Fourier series; satisfaction of the von Karman [2]/Marguerre [3] compatibility equation enables a corresponding series for the (Airy) force function to be de®ned in terms of the coecients in the displacement series. Until this stage, the analysis has proceeded along the lines adopted by Levy [4].

* Corresponding author. Tel.: +44-121-414-5054; fax: +44121-414-3675; e-mail: [email protected].

Levy's method would now proceed to direct substitution into the von Karman/Marguerre equilibrium equation; on the other hand, the method in Ref. [1] achieves the same e€ect indirectly by direct minimization of the total potential energy of the system comprising the plate and the loading arrangement. The type of formulation proposed in Ref. [1] is well suited to the development of a general computer program, and the resulting program is presented in Ref. [5]. The type of plate problem covered by the method described above is limited to cases with symmetry about both plate axes, and also there must be zero lateral de¯ection along all four edges. Both pairs of opposite edges remain straight and parallel. A need arose to analyse outstands, i.e. rectangular plates supported along three edges, and free along the fourth. Clearly such a case is not allowed for by the method of Ref. [1]. This prompted the development of the method proposed in the present paper. In its basic features it is the same as the previous one; however, the series chosen to represent the lateral displacement,

0045-7949/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 1 7 5 - 8

334

G.H. Little / Computers and Structures 71 (1999) 333±352

whilst maintaining symmetry about the transverse axis, allows any form for the transverse pro®le. It is shown in Ref. [6] that the transverse edges still remain straight, but can now rotate in-plane. The in-plane loading applied to these edges in the general case comprises a longitudinal force on a prescribed line of action, together with a moment. Thus, for any plate, the in-plane loading can be uniform or can vary across the width. Full consideration is given to providing prescribed load or prescribed displacement. Rhodes et al. [11±13] also analysed rectangular plates in compression, with the loaded edges simply supported and remaining straight in-plane, and with various boundary conditions along the unloaded edges, and they also used an energy method together with the von Karman/Marguerre compatibility equation. However, their post-buckling solution was based on a single term for the de¯ected shape, this term being a good approximation to the critical buckling mode. Thus, the form of the de¯ected shape was not allowed to change during post-buckling, and, as they accept [11], this may lead to inaccuracy as de¯ections become large. The present method is a full multiterm solution, readily allowing for any number of terms in the series for w, and in the similar series for the initial imperfection w0. Thus, the present method fully allows for changes in the mode shape as loading progresses. Having said that, Rhodes et al. compared their results with experimental results [13], including their own, and those of others, and obtained reasonable agreement. Work on the outstand problem also includes that of Ronalds [8] and Ronalds and Chapman [14], who developed a simpli®ed but e€ective energy method, and compared its predictions with results they obtained using a ®nite element package. Louca and Harding [15] also studied outstands using ®nite elements, including individual outstands, and also an outstand attached to an equivalent width of plating. Cases covered by the proposed method include plates supported along all four edges, outstands (one free edge), wide columns (two free, opposite, edges, inplane compression) and deep beams (two free, opposite, edges, pure bending in-plane). Examples of all these cases are included here. Comparisons with previous work, and with results from a ®nite element model, are shown to validate the method. The method will also be shown to be very rapid in operation. Two ways in which the scope of the method can be increased will now be mentioned brie¯y. Only in-plane loading is dealt with here. Lateral loading can be added simply by including an appropriate potential energy term. In this respect, the ability of the method to deal with non-symmetry about the longitudinal axis would make it very suitable for the analysis of plates with, for example, hydrostatic pressure loading. For

the sake of clarity, the analysis is restricted here to isotropic plates. The intention is to include orthotropic behaviour in the near future, but the modi®cations required are more involved than for lateral loading. Of course, an enormous amount of work on plate behaviour has been done over the yearsÐonly a brief discussion has been included here. In recent years, modern methods of computer analysis have become dominant in plate analysisÐthe ®nite element method, e.g. as used in Refs. [14] and [15], or the ®nite strip method, e.g. Ref. [16]. The current method is proposed as an alternative to these more general methods of analysis, for problems of the type described above. It will be shown to be very much more ecient, both in time and storage requirements, than the ®nite element method for these problems. It is hoped that it may prove useful in the design oce, and also as a research tool, both as a direct source of research data, and to provide comparative results for ®nite element/®nite strip models.

2. Plate details 2.1. Geometry and loading A typical rectangular plate as considered in the analysis is illustrated in Fig. 1. Each plate ABCD is of length a, width d, and uniform thickness t. The x, y, and z axes are longitudinal, transverse, and lateral, respectively, and z is zero at mid-thickness. There is symmetry about the Oy axis, and the Ox axis passes through the two corners A and B. An in-plane compressive force P is applied in the x direction to the transverse edges, at points E and F in Fig. 1, the line of action being prescribed at y = h. An in-plane moment M about an axis through E and about an axis through F may also be applied: in some examples this moment is prescribed to be zero. The displacements corresponding to P and M are the shortening ea and the in-plane edge rotation y, respectively, Fig. 1(b). At any position, (x, y), the in-plane displacements at mid-thickness are u and v in the x and y directions, respectively, and w represents the total displacement from the ¯at state. As there may be an out-of-¯atness w0 before any load is applied, the lateral de¯ection is (wÿ w0). Plane stress is assumed, and so the stress and strain vectors at any point (x, y, z) are: sT ˆ fsx

sy

txy g

and

ET ˆ fEx

Ey

gxy g;

…1†

and the material is assumed to be elastic, isotropic and homogeneous.

G.H. Little / Computers and Structures 71 (1999) 333±352

335

yˆh u ˆ 3ea=2:

…4†

It is shown in Ref. [6] that this boundary condition is automatically satis®ed by the proposed method. In-plane, static 9 Zd > > Nx dy …a† > P ˆÿ > > > 0 = Zd : …5† M ˆÿ Nx …y ÿ h† dy …b† > > > > 0 > > ; Nxy ˆ 0 …c† Out-of-plane:

 kinematic w ˆ 0 : static Mx ˆ 0

…6†

2.2.2. Along the longitudinal edges, y = 0 and y = d The only boundary condition along these edges assumed the same for all examples, is the condition of zero in-plane loading. Thus, along y = 0 and y = d:  Ny ˆ 0 …7† Nxy ˆ 0

Fig. 1. Plate details: (a) plate dimensions, coordinate axes and loading; (b) displacement of loaded edge.

At any position (x, y), the stress resultants (per unit length or width as appropriate) are: Z ‡t=2 Z ‡t=2 Nˆ s dz and M ˆ zs s dz: …2† ÿt=2

ÿt=2

Except where stated otherwise, any vector contains three components: a longitudinal component (subscript x), a transverse component (subscript y), and an inplane shearing, or twisting, component (subscript xy). 2.2. Boundary conditions 2.2.1. Along the transverse edges, x =2 a/2 The transverse edges (AC and BD in Fig. 1) remain straight, but are allowed to rotate in-plane, the angle of rotation being y, corresponding to the applied inplane moment M. The shortening strain along the line of action of the load P is e. Either M or y may be prescribed, and likewise for either P or e. These edges have zero in-plane shear loading, and are simply supported against out-of-plane movement. Therefore, the boundary conditions along x =2 a/2 are stated thus: In-plane, kinematic: @u ˆ 3y; @y

…3†

which are static conditions. These edges therefore adopt a de¯ected shape in-plane (symmetrical about Oy) which is consistent with zero in-plane loading along the length. All the out-of-plane boundary conditions along these edges depend on the particular plate being analysed, and will be stated later, when particular examples are analysed. 3. Method of analysis 3.1. Basis The basis of the method of analysis is to represent w and w0 by kinematically admissible series, choose a compatible series to represent the force function F, ensure that it satis®es the kinematic in-plane boundary conditions, and ®nally directly minimize the total potential energy to obtain the de¯ected shape which best satis®es out-of-plane equilibrium and static boundary conditions. A load stage is de®ned by prescribing a convenient in-plane loading parameter, and the minimisation process takes place at each load stage. The basic equations for plate analysis are given in the Appendix [Eqs. (A.1)±(A.10)]. 3.2. The compatibility equation and in-plane boundary conditions Firstly de®ne Y as follows: Y ˆ y=d:

…8†

336

G.H. Little / Computers and Structures 71 (1999) 333±352

In the implementation of the method, the gi( y) functions of Ref. [6] are taken to be polynomials in Y, and each gi cos(mipx/a) term in Ref. [6] is now written as the sum of the separate polynomial terms, thus:   Nw X mi px wˆ wi Yni cos ; a iˆ1 mi ˆ 1; 3; 5 . . . ; mmax ; ni ˆ 0; 1; 2 . . . ; nmax

…9†

where wi are the displacement coecients constituting the variables of the problem, and Nw is the total number of terms in the series for w. (Note that, unlike Eq. (7) of Ref. [6], more than one term in the series may have the same value of mi.) The initial imperfection w0 is represented using the same series terms, thus:   Nw X mi px w0 ˆ ; …10† …w0 †i Yni cos a iˆ1 which enables a wide range of imperfection shapes to be studied. Now the di€erential operator in Eq. (A.9) is given by:   Nw X Nw 1 p 2X N …w† ÿ N …w0 † ˆ w^ ij Yni ‡nj ÿ2 2 ad iˆ1 iˆ1  …mi ÿ mj †px  …mij †1 cos a  …mi ‡ mj †px ‡ …mij †2 cos ; …11† a where

9 …mij †1 ˆ m2i nj …nj ÿ 1† ÿ mi ni mj nj > > = …mij †2 ˆ m2i nj …nj ÿ 1† ‡ mi ni mj nj and > > ; w^ ij ˆ wi wj ÿ …w0 †i …w0 †j

Eq. (11) is written thus:   max   1 p 2 mX 2ppx cp cos N …w† ÿ N …w0 † ˆ 2 ad a pˆ0 where X cp ˆ w^ ij mij Yni ‡nj ÿ2

qn X qˆ0

cpq Yq

In order to satisfy the compatibility Eq. (A.10), the Zp functions are related to cp as follows: 2 00 4 Z0000 p ÿ 2kp Zp ‡ kp Zp ˆ cp ;

p ˆ 0; 1; 2; . . . qn

…17†

where each prime denotes di€erentiation once by Y and kp ˆ

2pp : a

…18†

3.2.1. Consider Z0 For p = 0, and noting that the cp functions are polynomials in Y, integration of Eq. (17) is straightforward. Only Z000 is required for the stress resultants (Nx only). It has been shown in Ref. [6] that Z000 is given as follows:  "   X imax a 2 2 h 1 2 2 00 Z0 ˆ 2 ÿY yÿe ‡ m …g ÿ g20i †; p a d 2 i i iˆ1 …19† in which gi( y) is a general function for the transverse variation of w ( g0i is the corresponding function for the initial out-of-¯atness). Eq. (19) may be written as the polynomial: 2 "   2n max X a 2 h 00 Z0 ˆ 2 ÿY yÿe ‡ Z000q Yq : …20† p a d qˆ0

…12† 3.2.2. Consider Zp, p>0 Each Zp is found from Eq. (17) and consists of the complementary function (Zp)CF plus the particular integral (Zp)PI. Thus: …13†

Zp ˆ …Zp †CF ‡ …Zp †PI where

…14†

where the summation extends over all ij for which vmi2mjv = 2p, and mij is either (mij)1 or (mij)2 as appropriate. The expression for cp can be written cp ˆ

includes the sum of one or more w^ ij terms. Some cpq may be zero. The series for F is:  max   Et p 2 mX 2ppx Fˆ : …16† Zp cos 2 a a pˆ0

…15†

where qn = (2nmax ÿ2). Each non-zero quantity cpq

…Zp †CF ˆ …bp1 ‡ bp2 Y†e‡kp Y ‡ …bp3 ‡ bp4 Y†eÿkp Y and …Zp †PI ˆ

qn X qˆ0

Zpq Yq :

}

…21†

The quantities Zpq are determined so as to ensure that Zp = (Zp)PI is a solution of Eq. (17), and the quantities bp1, bp2, bp3, bp4 are then determined from the zero inplane loading conditions along Y = 0 and Y = 1,

G.H. Little / Computers and Structures 71 (1999) 333±352

Eq. (7). The algorithm for determining Zpq is in the following. 1: For q ˆ qn ; and q ˆ qn ÿ 1 : Zpq ˆ

cpq =k4p :

2: For q ˆ qn ÿ 2; and q ˆ qn ÿ 3 : Zpq ˆ ‰cpq ‡ 2k2p …q ‡ 1†…q ‡ 2†Zpq‡2 Š=k4p : 3: For q ˆ qn ÿ 4; to q ˆ 0 : Zpq ˆ ‰cpq ‡ …q ‡ 1†…q ‡ 2†…2k2p Zpq‡2 ÿ …q ‡ 3†…q ‡ 4†Zpq‡4 †Š=k4p :

}

Cp ˆ

2

1 61 ˆ6 40 0

along these edges; then Eq. (21) gives 3

1 …1 ÿ kp †eÿkp

7 7 7 5

Z000 ˆ

Zpqn

…26†

…27†

qX n ‡4 qˆ0

w f0q …Xi †Y f0q …Y†;

…28†

secondly Zp, in which the bp terms are the last four (see Table 1)

…23†

Zp ˆ

qX n ‡4 qˆ0

w fpq …Xi †Y fpq …Y†:

…29†

In the above equations, Xi represents a typical variableÐ the variables include all the wi together with e and y (the latter two only appear in Z0). The nature of the functions in the above equations depends on the value of q, and they are speci®ed in Table 1.

which is A p
B p ˆ ÿCC p :

…24†

Also: Table 1 Speci®cation of functions in Eqs. (28) and (29) q value

Z 00 ( p = 0) w f0q

0E qE qn q = qn + 1 q = qn + 2(=2nmax)

Z 000q Z 000q Z 000q

q = qn + 3

e

q = qn + 4

y

…25†

Note that A p ; D and F p are constant matrices and therefore only have to be evaluated once for a given plate. For a given p-value, F p relates the complementary function quantities directly to the particular integral quantities. The opportunity is now taken to express both Z000 and Zp in a common, complete, form as follows; ®rstly Z000 , in which the y and e terms are the last two (see Table 1)

(A.6) and (16) show that both Zp and Z0p must be zero

eÿkp

> > > .. > > . > > > ;

ˆ ÿF F p
E p where F p ˆ A ÿ1 p
D:

For zero Ny and Nxy along Y = 0 and Y = 1, Eqs.

6 e‡kp e‡kp eÿkp 6 6 4 kp 1 ÿkp kp e‡kp …1 ‡ kp †e‡kp ÿkp eÿkp 8 9 0 1 …Zp †PI …0† > bp1 > > > > < …Z † …1† > = B bp2 C p PI B C B Cˆÿ 0 > @ bp3 A > …Zp †PI …0† > > > > : ; …Zp †0PI …1† bp4

... ... ... ...

B p ˆ ÿA Aÿ1 Aÿ1 Ep p
C p ˆ ÿ…A p
D †E

the former as Eq. (22) is applied.

0

0 1 0 3

9 > > > > > > > > =

Thus:

be used for both cpq and Zpq, the latter over-writing

1

0 1 0 2

Zp0 Zp1 Zp2 Zp3

C p ˆ D
E p:

problem, it is worth noting that the same storage can

0

0 1 1 1

8 > > > 3> > 0 > > > < 7 17 0 5> > > qn > > > > > :

which is

…22†

method of analysis is very unlikely to have a storage

1

p PI

…Z †0 …0† > > > > : …Zp †PI 0 …1† ; p PI

Although a computer program based on the current

2

337

9 8 …Z † …0† > > > = < …Zp †PI …1† >

Zp (1E pE mmax)

Y f0q

w fpq

Y fpq

Yq Yq Yq   a 2 ÿ2 p     4 a 2 h ÿY a p d

Zpq bp1 bp2

Yq e ‡ kp Y Ye ‡ kp Y

bp3

e ÿ kp Y

bp4

Ye ÿ kp Y

338

G.H. Little / Computers and Structures 71 (1999) 333±352

As a result of the above analysis, the force function F, as given in Eq. (16), in conjunction with the lateral displacement w given by Eq. (9), exactly satis®es inplane equilibrium, together with all the in-plane boundary conditions, except those in Eq. (5a) and (b). 3.3. Energy formulation and out-of-plane equilibrium Out-of-plane equilibrium has not yet been satis®ed. Also, the static boundary conditions in Eq. (5a) and (b) are not yet satis®ed, and some problems are likely to have static out-of-plane boundary conditions which are not automatically satis®ed by the displacement functions. For the displacement functions chosen, in order to ®nd the de¯ected shape at each load stage which best satis®es all these requirements, the total potential energy will be minimised. The membrane strain energy f Im in the plate is: Z ‡a=2 Z d 1 fIm ˆ ‰N2 ‡ N2y ÿ 2nNx Ny 2Et ÿa=2 0 x ‡ 2…1 ‡

n†N2xy Š dx dy:

…30†

From Eqs. (A.6) and (16), the in-plane stress resultants N are:     max   1 Et p 2 1 2 mX 2ppx Nx ˆ 2 F;YY ˆ ; Z00p cos d 2 a d a pˆ0      max  Et p 2 2p 2 mX 2ppx Ny ˆ F;xx ˆ ÿ Zp p2 cos ; 2 a a a pˆ1     max   1 Et p 2 2p mX 2ppx : Nxy ˆ ÿ F;xY ˆ Z0p p sin d 2 a ad pˆ1 a

}

…31† Note that the summations are, of for Ny and Nxy, but from p = 0 for Substitute from Eq. (31) for N integrate with respect to x, which to vanish, to leave the following: fIm ˆ

m max X p4 Et ‡ …Ip1 ‡ Ip2 ‡ Ip3 ‡ Ip4 †Š ‰I 0 16a4 ad pˆ1

where I0 ˆ 2a2 Ip1

course, from p = 1 Nx . into Eq. (30), and causes many terms

Z

1

0

…Z000 †2 dY;

16p4 p4 ˆ a2

Z

1

0

Z2p dY;

Ip2 ˆ 8p2 …1 ‡ n†p2 Ip3 ˆ a2

Z 0

1

Z 0

1

…Z0p †2 dY;

…Z00p †2 dY

Ip4 ˆ 8p2 np2

Z

1 0

Z00p Zp dY:

}

…32†

The 8p 2nIp4 term and the 8p 2nIp2 term combine to produce a term which, after integration, and satisfaction of the zero in-plane loading boundary condition along Y = 0 and Y = 1, vanishes. This then leaves an expression for the membrane strain energy which is independent of n. Using Eq. (29), Ip1 can be expressed as: Ip1 ˆ

qX n ‡4 q n ‡4 X

kpqr w fpq w fpr

where each kpqr is the result of an integration with respect to Y of a product of two Yfpq functions. It is a straightforward matter to write a general subroutine to calculate these. Clearly similar equations can be written for Ip2±4, and then Ip1±4 can be summed to give: Ip1 ‡ Ip2 ‡ Ip3 ‡ Ip4 ˆ

qX n ‡4 q n ‡4 X

Kpqr w fpq w fpr

…35†

qˆ0 rˆ0

where Kpqr simply represents the sum of the kpqr corresponding to Ip1±4. The Kpqr are invariant for any particular analysis, so are calculated only once. After similar integrations for I0, the membrane strain energy can now be written ®nally as: fIm ˆ

n ‡4 q n ‡4 max q X X p4 Et mX Kpqr w fpq w fpr : 16a4 ad pˆ0 qˆ0 rˆ0

…36†

The bending strain energy f Ib in the plate is: Z Z D ‡a=2 d 2 fIb ˆ ‰w ‡ w 2;yy ‡ 2nw ;xx w ;yy 2 ÿa=2 0 ;xx ‡2…1 ÿ n†w 2;xy Š dx dy

…37†

where Dˆ

Et3 12…1 ÿ n2 †

…38†

and, from Eq. (9), and using the notation in Eq. (A.1)   2 X  imax p mi px 2 ni w ;xx ˆ ; ÿw i mi Y cos a a iˆ1  w ;yy ˆ  w ;xy ˆ

1 d2 p ad

where w i ˆ wi ÿ …w0 †i :

…33†

…34†

qˆ0 rˆ0

X imax iˆ1

X imax iˆ1

  mi px ; w i ni …ni ÿ 1†Yni ÿ2 cos a ÿw i mi ni Yni ÿ1 sin



 mi px ; a

}

…39†

Substitute from Eq. (39) into Eq. (37); a number of terms such as cos(mipx/a)cos(mjpx/a) are produced. On integration with respect to x, all such terms vanish except those for which mi = mj. The resulting ex-

G.H. Little / Computers and Structures 71 (1999) 333±352

pression for the bending strain energy is fIb ˆ

imax X imax 2X

Dp 2ad

Kij w i w j

…40†

iˆ1 jˆ1

where, for mi = mj, Kij = 0, and for mi$mj:     p2 m4i a2 ni …ni ÿ 1†nj …nj ÿ 1† Kij ˆ 2 ‡ 2 2a ni ‡ nj ‡ 1 2p ni ‡ nj ÿ 3  2    m nj …nj ÿ 1† m2i ni nj ÿn i ‡ …1 ÿ n† : …41† ni ‡ nj ÿ 1 ni ‡ nj ÿ 1 Precautions are necessary, for certain combinations of ni and nj, to avoid dividing by zero. The Kij, like the Kpqr, are invariant for any given plate problem, and so are calculated once only. The potential energy f E of the loads is: fE ˆ ÿPae ÿ 2My:

…42†

Any lateral loading system could readily be added to the above, provided it is symmetrical about Oy. The total potential energy f of the system comprising the plate and the loading arrangement is f ˆ fIm ‡ fIb ‡ fE

…43†

and the best equilibrium shape is found by minimising the total potential energy with respect to the variables wi, e, and y, thus: for i ˆ 1; 2; . . . ; N : f ˆ0 …a† ) w

also : f;e ˆ fm;e f;y ˆ fm;y ÿ 2M

;wi

ÿPa ˆ 0 ˆ0

…b† …c†

…44†

Eqs. (42)±(44) are speci®cally for prescribed load, i.e. a load stage is de®ned by prescribing the values of P and M, and the corresponding values of e and y arise out of the minimization process. Sometimes it is more convenient to prescribe the corresponding displacement, indeed it can be essential to do so if snap-through behaviour is possible. To prescribe e, remove the Pae term from fE, remove Eq. (44b), minimise the total potential energy with respect to wi and y, and then apply Eq. (44b) to ®nd P. To prescribe y, the procedure is similar. Having the ability to prescribe either load or displacement in each case, increases the versatility of the method. Thus, for example, for the load P to be applied on a ®xed line of action, prescribe M to be zero. To simulate loading between ¯at platens in a testing machine, say, prescribe y to be zero. Then the load P will still be at a ®xed line of action, but there will be an additional moment M which will vary as loading progresses, and which will mean that the total load is statically equivalent to a point load on a line of action that varies as loading progresses. Other combinations are, of course, possible. There is a simple, obvious, yet powerful, check that can now be made on the theory and the coding: analyse a typical plate using prescribed-e loading, note the

339

calculated values of P and use these to analyse the same plate with prescribed-P loading. The calculated e values should agree exactly with those of the ®rst analysis. This check has been satis®ed, and also the equivalent M±y check. Other simple checks involve prescribing di€erent sets of statically equivalent P, M, and d values, and seeing that the same de¯ected shape is predicted; similarly, di€erent prescribed sets of kinematically equivalent e, y, and d values, should give the same de¯ected shape. All these checks have been satis®ed. The basic sequence of calculations is: 1. De®ne the plate geometry and material properties. 2. De®ne the terms in Eq. (9). The terms used depend on the support conditions along Y = 0 and Y = 1, and on the degree of accuracy required. De®ne also the initial imperfection. 3. Calculate the invariant quantities, including the Kpqr for the membrane strain energy, and the Kij for the bending strain energy. 4. De®ne the load stage by prescribing M or y, and P or e. 5. Minimise the total potential energy f using a quasiNewton (or variable metric) algorithm. Such a method was used in Ref. [1], and it requires, for any set of values in the variables, the corresponding values of f and D f ; where D f is the vector of `slopes' of f with respect to each variable, i.e.:  @f T T Df ˆ ˆ ff;w1 ; f;w2 ; . . . f;e ; f;y g: …45† @X It is not necessary to calculate the matrix of second derivatives of f with respect to the variables (the sti€ness matrix)Ðthe method builds up an estimate of this automatically. Convergence to the minimum is usually very ecient. 6. Output the results for the load stage. 7. Repeat steps (4), (5) and (6) until a sucient number of load stages has been covered. Calculating the contributions of the bending strain energy and the potential energy of the loads towards f and D f follows in a straightforward manner from Eqs. (40) and (42). The contributions of the membrane strain energy are more dicult to evaluate, but the basic procedure is: For p = 0, 1, 2 . . ., mmax 1. from the given values of wi, e and y (the latter two only feature for p = 0), calculate cpq and cpq,Xi, all q, iÐEqs. (10), (11) and (13); 2. calculate Zpq and Zpq,Xi, all q, iÐEq. (22); 3. calculate bp1±4 and bp1±4,Xi, all iÐEqs. (23) and (25);

340

G.H. Little / Computers and Structures 71 (1999) 333±352

4. for the current p value, calculate the contributions towards f Im and D f Im ; and sum into the totalÐ Eq. (36).

3.4. Coupled polynomial terms for transverse form of displacement In the general case, the variables wi are independent of each other. For some problems, the wi are not independent, and must combine so that, for example, certain kinematic boundary conditions are satis®ed. In these cases, the Nw variables wi are replaced by the NW variables Wi, where NW < Nw, and the relationship between Wi and wi is ®xed as: wˆG
W

…46†

where G is a Nw  NW matrix with constant coecients. Thus, the total potential energy is now minimized with respect to W, e and y. Each time y and D f are required, w is calculated from W using Eq. (46), the calculations then take place as for the general case, and ®nally @f/@W is found from @f/@w using: @f @f ˆ GT : @W @w

…47†

Examples using this technique will be given in the next section. 4. Applications The method of analysis described above is now applied to a series of plate problems, di€ering according to the degree of out-of-plane support along the longitudinal edges, and to the manner of loading. Some of the results are useful in the sense of validating the method by comparison with previous work, and some serve to provide new insight into the deformation of plates having one or two longitudinal edges free. The out-of-plane support along an edge is denoted as S (simply supported), or C (clamped), or F (free). For example, a plate stated to be SSSF has edges x = 0 and x = a both simply supported (the ®rst two lettersÐall examples, of course, have those two edges simply supported), with edge Y = 0 simply supported, and Y = 1 free. In cases with an initial imperfection, it is shaped as the ®rst term in the series for w, which is always similar to the initial buckling mode. Except for the comparisons with Yamaki [7], all plates are analysed using prescribed shortening (or prescribed rotation y for the deep beam), with a standard normalised increment size of 0.1. For initially ¯at plates, and some small imperfection amplitudes, the increment size is reduced in the region of the buckling load. For initially ¯at plates, the displacement function

coecient w1 corresponding to the ®rst term in the series for w is given a small positive initial guess of 0.001t at each load stage until buckling occurs. Before this occurs, the energy minimization process results in this initial guess being returned to zero. Graphs are drawn by joining calculated points by straight linesÐ the resulting smoothness of the curves is an indication that sucient points have been calculated. Plots of transverse pro®les are always for the cross-section at mid-length, and a pro®le is given for each load stage. For comparisons with Yamaki, Poisson's ratio is taken as 1/3; for all other cases it is 0.3. The method is applied via a FORTRAN computer program running in double precision on a Digital AlphaStation 600 5/266 hardware platform. An indication of the eciency of the method is given by the fact that the total CPU time required for all the results given in Figs. 3±9 (20 plate analyses), was only 2 s. A direct comparison with running a commercial ®nite element package will be given in Section 4.4. The results are presented in normalised form, as follows: 9 e ke > > e ˆ where Ed ˆ ; > > > Ed …d=t†2 > > > > > P k > e  > P ˆ E; > where sd ˆ > 2 > sd td …d=t† > > > > > > y ky >  > y ˆ where yd ˆ ; = 2 yd …d=t† …48† 3 > > M k Etd > y > M ˆ ;> where Md ˆ > Md …d=t†2 6a > > > > > > w max >  > wmax ˆ ; > > t > > > > > f >  > f ˆ ; …sd Ed atd† and the values given to ke and ky depend on the problem. Note that wmax is the lateral displacement at the point (0, d/2) for comparisons with Yamaki, and is at the point (0, d) for outstands.

4.1. Simply-supported plates (SSSS) Application of the method to plates simply supported against out-of-plane movement along all four edges requires the use of the coupled polynomial technique explained in the previous section. Also, the case of uniform longitudinal compression is analysed, which enables the results to be compared directly with those of Yamaki [7] for the same problem. A suitable displacement function series for w is as follows:

G.H. Little / Computers and Structures 71 (1999) 333±352

wˆ "    #   m `max max X X 1 2` 1 2` mpx cos : Wm` ÿ Yÿ 2 2 a mˆ1;3 `ˆ1

malisation is done using ke = p 2, which makes P* equivalent to Yamaki's l. For each set of displacement function terms, w*max and f* are given in the table for each load stage, and Yamaki's values of w*max are also included. The present results show the total potential energy decreasing monotonically as the number of terms used increases, as would be expected. The results for the 9-term set agree closely with those for the 6term set, indicating that there is little point in using more terms. Agreement with Yamaki is suciently close. Exact agreement should not be obtained, for two reasons. Firstly, unlike the present method, Yamaki makes use of Coan's method, and does not exactly satisfy either in-plane equilibrium, or the zero load condition along the longitudinal edges. Secondly, Yamaki's initial imperfection is of doubly sinusoidal shape, not quite the same as the (1, 1) function, used here for the imperfection.

…49†

Functions of the type in square brackets satisfy the kinematic boundary conditions along Y = 0 and Y = 1, namely w = 0 along both those edges. The functions also provide symmetry about the centre line Y = 1/2, which is appropriate for the case of uniform compression. In order to employ the above series for w, ®rstly decide on the terms to be included, then multiply out the Y functions to obtain the individual polynomial terms as in Eq. (9), and ®nally determine the G array relating the wi to the Wm` : This array is included in the data. A square plate has been analysed with di€erent combinations of functions used to represent w. These combinations are: 1 term …1; 1†; 2 4 6 9

terms terms terms terms

…1; 1†; …1; 2†; …1; 1†; …1; 2†; …3; 1†; …3; 2†; …1; 1†; …1; 2†; …1; 3†; …3; 1†; …3; 2†; …3; 3†; …1; 1†; …1; 2†; …1; 3†; …3; 1†; …3; 2†; …3; 3†; …5; 1†; …5; 2†; …5; 3†

}

341

4.2. Clamped longitudinal edges (SSCC) Now consider uniform compression applied to a plate with its longitudinal edges clamped against outof-plane movement (Yamaki's case IIb, Ref. [7]). A suitable displacement function series for w is as follows:   m `max max X X mpx wˆ : …51† Wm` ‰Y` …1 ÿ Y†` Šcos a mˆ1;3 `ˆ2

…50†

where, for example, (3, 1) means the term with m = 3, ` ˆ 1: The rotation y is prescribed to be zero (symmetry); the load P* is prescribed, with line of action h/d = 0.5. Symmetry then requires M to be zero at equilibrium, and the calculated values of M are all extremely small. An initial imperfection of amplitude D0 = 0.1t, and shaped as the (1, 1) term, is included in each analysis. The results are given in Table 2. To enable direct comparison with the results of Yamaki (his case Ib), nor-

Functions of the type in square brackets satisfy the clamped edge conditions along Y = 0 and Y = 1, and provide symmetry about the centre line Y = 1/2, which is appropriate for the case of uniform compression. The above series for w is applied in a similar way to the simply supported case. For this case, a plate of aspect ratio a = 0.5 has been analysed with di€erent

Table 2 Simply supported plate, comparison with Yamaki [7] P*

1 term w*max

f*

2 terms w*max

f*

4 terms w*max

f*

6 terms w *max

f*

9 terms w*max

f*

Yamaki w*max

0.212 0.278 0.342 0.390 0.438 0.546 0.672 0.812 0.963 1.126 1.300 1.486

0.210 0.312 0.518 0.758 1.013 1.504 1.954 2.365 2.743 3.102 3.445 3.777

ÿ0.023 ÿ0.039 ÿ0.061 ÿ0.081 ÿ0.106 ÿ0.182 ÿ0.310 ÿ0.501 ÿ0.767 ÿ1.123 ÿ1.583 ÿ2.166

0.229 0.358 0.606 0.849 1.084 1.534 1.961 2.363 2.745 3.118 3.484 3.849

ÿ0.023 ÿ0.040 ÿ0.061 ÿ0.082 ÿ0.108 ÿ0.185 ÿ0.312 ÿ0.502 ÿ0.767 ÿ1.127 ÿ1.599 ÿ2.208

0.229 0.358 0.605 0.848 1.083 1.534 1.961 2.363 2.748 3.130 3.518 3.920

ÿ 0.023 ÿ 0.040 ÿ 0.061 ÿ 0.082 ÿ 0.108 ÿ 0.186 ÿ 0.316 ÿ 0.513 ÿ 0.794 ÿ 1.182 ÿ 1.705 ÿ 2.399

0.229 0.358 0.605 0.848 1.085 1.542 1.982 2.403 2.811 3.223 3.641 4.060

ÿ0.023 ÿ0.040 ÿ0.061 ÿ0.082 ÿ0.108 ÿ0.186 ÿ0.318 ÿ0.522 ÿ0.816 ÿ1.229 ÿ1.790 ÿ2.539

0.229 0.358 0.605 0.848 1.085 1.542 1.982 2.405 2.814 3.230 3.644 4.070

ÿ0.023 ÿ0.040 ÿ0.061 ÿ0.082 ÿ0.108 ÿ0.186 ÿ0.318 ÿ0.522 ÿ0.816 ÿ1.230 ÿ1.795 ÿ2.551

0.225 0.349 0.596 0.839 1.076 1.535 1.971 2.388 2.793 3.191 3.587 3.985

342

G.H. Little / Computers and Structures 71 (1999) 333±352

Table 3 Clamped plate, comparison with Yamaki [7] P*

1 term w*max

f*

2 terms w*max

f*

4 terms w*max

f*

6 terms w*max

f*

Yamaki w*max

0.386 0.524 0.675 0.799 0.934 1.269 1.721 2.320

0.207 0.313 0.517 0.714 0.912 1.309 1.721 2.154

ÿ0.076 ÿ0.142 ÿ0.242 ÿ0.353 ÿ0.506 ÿ1.046 ÿ2.149 ÿ4.281

0.206 0.312 0.522 0.728 0.937 1.368 1.823 2.298

ÿ0.076 ÿ0.142 ÿ0.243 ÿ0.355 ÿ0.514 ÿ1.099 ÿ2.386 ÿ5.094

0.206 0.312 0.522 0.727 0.936 1.366 1.818 2.287

ÿ0.076 ÿ0.142 ÿ0.243 ÿ0.355 ÿ0.514 ÿ1.099 ÿ2.387 ÿ5.098

0.206 0.312 0.522 0.727 0.936 1.368 1.824 2.309

ÿ 0.076 ÿ 0.142 ÿ 0.243 ÿ 0.355 ÿ 0.514 ÿ 1.099 ÿ 2.388 ÿ 5.110

0.202 0.304 0.509 0.715 0.923 1.352 1.806 2.296

combinations of functions used to represent w. These …m; `† combinations are: 1 term …1; 2†; 2 terms …1; 2†; …1; 3†; 4 terms …1; 2†; …1; 3†; …3; 2†; …3; 3†; 6 terms …1; 2†; …1; 3†; …3; 2†; …3; 3†; …5; 2†; …5; 3†:

}

4.3. Outstands (SSSF) For the outstand problem, with Y = 0 simply supported and Y = 1 free, it is not necessary to use coupled polynomial terms; the basic series for w in Eq. (9) is used, with the following sets of (m, n) terms:

…52†

1 term …1; 1†; 2 terms …1; 1†; …1; 2†; 4 terms …1; 1†; …1; 2†; …3; 1†; …3; 2†;

An initial imperfection of amplitude D0 = 0.1t, and shaped as the (1, 2) term, is included in each analysis. Other details of the analysis are as for simply supported plates. The results are presented in Table 3 in a similar manner to Table 2, for simply supported plates, and similar comments apply. In particular, there is close agreement with Yamaki.

8 terms …1; 1†; …1; 2†; …1; 3†; …1; 4†; …3; 1†; …3; 2†; …3; 3†; …3; 4†:

}

…53†

Results for an outstand of aspect ratio a = 10 are given in Table 4. For these results

Table 4 Simply supported outstand P*

1 term w*max

0.100 0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.900 1.000 1.100 1.200 1.300 1.400 1.500 1.600 1.700 1.800 1.900 2.000

0.121 0.154 0.211 0.335 0.754 2.139 3.411 4.380 5.183 5.882 6.509 7.082 7.612 8.108 8.576 9.020 9.444 9.849 10.238 10.613

f* 0.005 0.020 0.045 0.080 0.125 0.177 0.233 0.289 0.347 0.406 0.465 0.525 0.586 0.647 0.710 0.773 0.836 0.901 0.966 1.032

2 terms w*max 0.121 0.154 0.211 0.335 0.754 2.139 3.411 4.379 5.182 5.880 6.506 7.079 7.608 8.104 8.572 9.015 9.438 9.842 10.231 10.605

f* 0.005 0.020 0.045 0.080 0.125 0.177 0.233 0.289 0.347 0.405 0.465 0.525 0.586 0.647 0.710 0.773 0.836 0.901 0.966 1.032

4 terms w*max 0.121 0.154 0.211 0.335 0.754 2.147 3.442 4.442 5.280 6.015 6.678 7.286 7.851 8.380 8.880 9.355 9.808 10.242 10.659 11.060

f* 0.005 0.020 0.045 0.080 0.125 0.177 0.233 0.289 0.347 0.405 0.464 0.524 0.585 0.646 0.707 0.770 0.833 0.896 0.960 1.025

8 terms w*max 0.121 0.154 0.212 0.336 0.757 2.152 3.446 4.446 5.283 6.018 6.681 7.289 7.854 8.383 8.883 9.358 9.811 10.245 10.662 11.063

f* 0.005 0.020 0.045 0.080 0.125 0.177 0.233 0.289 0.347 0.405 0.464 0.524 0.584 0.645 0.707 0.769 0.832 0.895 0.959 1.024

G.H. Little / Computers and Structures 71 (1999) 333±352

ke = 202  355/205,000, which enables a comparison to be made with those of Ronalds [8]. An initial imperfection of amplitude D0 = 0.1t, and shaped as the (1, 1) term, is included in each analysis. Longitudinal compression is applied by having h/d = 0.5, prescribing the shortening, and prescribing M = 0. The loaded edges are, therefore, free to rotate in-plane. Table 4 shows that the 8-term set agrees closely with the 4-term set, indicating convergence. An interesting point is the good level of accuracy provided by the 1-term results. The 8-term results are compared with those of Ronalds (scaled from Fig. 5(a) of Ref. [8]) in Fig. 2. Ronalds' analysis is of a simpli®ed nature, and yet her results agree reasonably well with the present 8-term solution for this case. Using the 8-term set, some additional results for outstands are given in Figs. 3±5. For these results ke is taken as follows:   1 p2 ke ˆ …54† 0:5…1 ÿ n† ‡ …1 ÿ n2 † 12a2 which ensures that sd is equal to the buckling stress corresponding to the ®rst term in the displacement function setÐEq. (53)Ðtogether with uniform compressive loading. Thus sd represents the buckling stress based on an assumed mode which has a straight line in the y-direction, and will be an approximation to the actual buckling stress. Timoshenko and Gere [10] quote this formula as being suciently accurate for long plates. The load-shortening curves in Fig. 3 are for a = 10, and include initial imperfections of amplitude D0 = 0,

343

0.05t, 0.1t and 1.0t. It can be seen from the curve for the initially ¯at plate that, for this long outstand, the approximate buckling stress is very accurate. This supports, for this plate geometry, the commonly made assumption that transverse pro®les remain straight. Further support for this can be seen in Fig. 3(b), in which the transverse pro®les at mid-length for the 0.1t imperfection are drawn for each load stage. The initially ¯at plate shows a small positive post-buckled sti€ness of approximately 6% of the full elastic (prebuckled) value. Fig. 4 shows results calculated on the same basis for outstands having a = 1. The load-shortening curve for the initially ¯at plate now indicates buckling to occur at P* 1 0.98. Thus, at this lower aspect ratio, buckling occurs with a curved transverse pro®le at a stress slightly less than that obtained by assuming the transverse pro®le to remain straight. Transverse pro®les for the 0.1t imperfection are given in Fig. 4(b) at each load stage, and they are noticeably curved. Of course, in order to satisfy the static boundary condition of zero transverse moment along each longitudinal edge, they must show negative curvature (transverse curvature = ÿ w,yy) at Y = 1, and zero curvature at Y = 0, as exhibited by the pro®les in Fig. 4(b). To clarify the transverse pro®les further, imagine, at each load stage, a straight line drawn from w* = 0 at Y = 0, to w*max at Y = 1. Deviations wdev(=wÿ Ywmax) of the pro®le from this straight line are plotted in Fig. 5, for each load stage. The zone of negative curvature near Y = 1 can clearly be seen. For most of the plate width from the supported

Fig. 2. Simply supported outstands, a = 10; D0 = 0.1t: comparison with Ronalds [8].

344

G.H. Little / Computers and Structures 71 (1999) 333±352

Fig. 3. Simply supported outstands, a = 10. (a) load-shortening curves; (b) transverse pro®les at mid-length, D0 = 0.1t. 20 pro®les are shown, drawn for e* = 0.1, 0.2, 0.3, . . . 2.0.

longitudinal edge, the curvature is positive. It should be noted that, for any aspect ratio, the pro®les have similar features, but for larger aspect ratios, deviations from the straight line are smaller (e.g. for the

a = 10 case, the deviations are an order of magnitude smaller). The post-buckled sti€ness of the a = 1 case is also positive, but is even smaller than that of the longer

G.H. Little / Computers and Structures 71 (1999) 333±352

345

Fig. 4. Simply supported outstands, a = 1: (a) load-shortening curves; (b) transverse pro®les at mid-length, D0 = 0.1t. 20 pro®les are shown, drawn for e* = 0.1, 0.2, 0.3, . . . 2.0.

plate, being only about 3% of the full elastic value. Looking now at the e€ect of imperfections on plate behaviour at the two aspect ratios, the imperfection± sensitivity is considerably greater for the shorter plate.

4.4. Wide columns (SSFF with uniform in-plane compression) For the wide column problem, both edges Y = 0 and Y = 1 are free; the basic series for w in Eq. (9) is used,

346

G.H. Little / Computers and Structures 71 (1999) 333±352

Fig. 5. Simply supported outstands, a = 1: transverse pro®les at mid-lengthÐdeviation from straight line joining edges; D0 = 0.1t. 20 pro®les are shown, drawn from e* = 0.1, 0.2, 0.3, . . . 2.0.

with the following set of (m, n) terms: 10 terms …1; 0†; …1; 1†; …1; 2†; …1; 3†; …1; 4†; …3; 0†; …3; 1†; …3; 2†; …3; 3†; …3; 4†

…55†

the two terms having n = 0 enabling both longitudinal edges to be free. For this case, we take ke as follows: ke ˆ p2 =12a2 ;

…56†

which makes sd equal to the Euler buckling stress. The moment M is prescribed to be zero (symmetry); the shortening e* is prescribed, with line of action h/d = 0.5. Results have been calculated for two aspect ratios, a = 1 and a = 10; load-shortening curves are given in Fig. 6(a) and 7(a) for various initial imperfection sizes D0, and transverse pro®les at mid-length in Fig. 6(b) and 7(b). For the a = 10 case, the behaviour is indistinguishable from that obtained using an ordinary elastic largede¯ection analysis of the imperfect column [9], as would be expected. In particular, the initially perfect case shows that buckling occurs accurately at the Euler buckling load, and the postbuckled sti€ness is zero. Also, the transverse pro®les seem nearly straight and parallel, but each incorporates the ``Poisson'' negative curvature (anti-clastic) necessary for zero transverse stress. [Deviation from the straight line joining the edges is 1(n/8a 2)w, i.e. 0.0004w approximately.] Turning now to the a = 1 case, Fig. 7, the ®rst thing to notice is that the initially ¯at plate buckles at P* 1 1.05. In other words, the Euler buckling stress

gives an underestimate of the true buckling stress. This is because Euler buckling incorporates a ``Poisson'' transverse negative curvature as mentioned in the previous paragraph, but in the current problem, the two loaded edges are held straight. For long plates these end e€ects are negligible, but for shorter plates they become a signi®cant sti€ening e€ect. Transverse curvature is very visible in the pro®les for the shorter plateÐFig. 7(b). The displacement functions used do not automatically provide symmetry about Y = 0.5; the fact that the calculated de¯ected shapes do have this symmetry, gives further support to the satisfactory operation of the method. Also, calculated values of y are extremely small, as should be the case. Imperfection±sensitivity is high, and of a similar magnitude for both aspect ratios. The a = 1 case,with D0 = 0.1t, was also analysed using a commercial ®nite element package, LUSAS [17] (the package also used for the analysis of outstands in Ref. [15]). The semi-loof shell element was used (QSL8), in an elastic non-linear analysis (total Langrangian). Because this case is symmetrical about both axes, only one quadrant of the plate was modelled, with symmetry conditions speci®ed along the two plate axes, and the actual boundary conditions along the two edges. A 2  2 element mesh over the quadrant gave very similar results to a 4  4 mesh, and both sets of ®nite element results were indistinguishable from the results produced by the proposed method and given in Fig. 7(a) and (b). LUSAS was

G.H. Little / Computers and Structures 71 (1999) 333±352

347

Fig. 6. Wide columns, a = 10: (a) load-shortening curves; (b) transverse pro®les at mid-length, D0 = 0.1t. 20 pro®les are shown, drawn from e* = 0.1, 0.2, 0.3, . . .,2.0 (rightmost curve).

mounted on a DEC 3000-M300X hardware platform, and so the current program was also run on it for this case. For the same number of load stages, and for the 2  2 mesh, LUSAS took 30 s of CPU time, whereas the proposed method required only 0.3 s, a factor of 100 times less. The savings in storage requirements were of

the same order. Bearing in mind that this particular case is one requiring, by its nature, only a few elements, the comparison would be even more striking for cases with a = 10 and/or without symmetry about the longitudinal plate axis, where many more elements would be required. (For long plates, however, a ®nite

348

G.H. Little / Computers and Structures 71 (1999) 333±352

Fig. 7. Wide columns, a = 1: (a) load-shortening curves; (b) transverse pro®les at mid-length, D0 = 0.1t. 20 pro®les are shown, drawn for e* = 0.1, 0.2, 0.3, . . . 2.0, (rightmost curve).

strip model would be signi®cantly more ecient than a ®nite element one.) 4.5. Deep beams (SSFF with pure bending in-plane) For the deep beam problem, both edges Y = 0 and Y = 1 are free as for the column problem, and the basic series for w in Eq. (9) is used, with the same set

of (m, n) terms as in Eq. (55). For these results ky is taken as follows: 2

3 p2 6 2…1 ÿ n† ‡ 7 a 6 12a2 7 ky ˆ 6 7  1=2 2 5 …1 ÿ n † 4 1 8a2 ‡ 2 …1 ÿ n† 3 p

…57†

G.H. Little / Computers and Structures 71 (1999) 333±352

349

Fig. 8. Deep beams, a = 10: (a) moment±rotation curves; (b) transverse pro®les at mid-length, D0 = 0.1t. 20 pro®les are shown, drawn for y* = 0.1, 0.2, 0.3, . . . 2.0 (rightmost curve).

which relates to a buckling mode corresponding to the ®rst two terms in the displacement function setÐ Eq. (55)Ðtogether with pure in-plane bending, i.e. a longitudinal stress distribution at buckling varying linearly across the width from +6 Md/td 2 at Y = 0 to ÿ6 Md/td 2 at Y = 1. Thus, Md represents the in-plane bending moment at buckling based on an assumed mode having a (general) straight line in the y-direction, and will be an approximation to the actual buckling moment. The rotation y is prescribed; the load P* is

prescribed to be zero. Results have again been calculated for the two aspect ratios, a = 1 and a = 10; the moment±rotation curves are given in Fig. 8(a) and 9(a), and transverse pro®les at mid-length in Fig. 8(b) and 9(b). For both aspect ratios 10 and 1, it can be seen that the approximate buckling moment is an overestimateÐby about 5% and 8%, respectively. Clearly, the buckling mode involves transverse curvature in both cases, but more so for the shorter plate, as

350

G.H. Little / Computers and Structures 71 (1999) 333±352

Fig. 9. Deep beams, a = 1: (a) moment±rotation curves; (b) transverse pro®les at mid-length, D0 = 0.1t. 20 pro®les are shown, drawn for y* = 0.1, 0.2, 0.3, . . . 2.0.

would be expected. The transverse curvature present in the longer plate is not visible in the pro®les in Fig. 8(b) (nor was it visible in the column pro®les for a = 10). This is simply because it is related to the longitudinal curvature, which is 100 times (a 2) smaller than for the shorter plate at the same displacement. The transverse pro®les [Fig. 8(b)] show the expected displacement and rotationÐmore displacement at the

compressed extreme ®bre. For the a = 1 case, the transverse pro®les become signi®cantly curved. Also, the displacement at the tension (Y = 0) extreme ®bre is small. For the longer plate, the post-buckled sti€ness is zero, but the imperfection±sensitivity is small. The shorter plate, however, shows a positive post-buckled sti€ness of about 14% of the full elastic value, but

G.H. Little / Computers and Structures 71 (1999) 333±352

imperfection±sensitivity is greater than for the longer plate.

5. Conclusions An ecient method has been proposed for analysing the large de¯ection behaviour of thin, elastic plates, simply supported along two opposite edges (the transverse edges), and with various types of support, including the free edge condition, along the other two edges. The method is based on the representation of lateral displacement by a series, satisfaction of the von Karman/Marguerre compatibility equation, and direct minimization of the total potential energy. In the present paper, only in-plane loading has been considered, but it is a straightforward matter to include lateral loading. The in-plane loading is applied along the two transverse edges, which remain straight but can rotate in-plane. In order to demonstrate the versatility of the method, results have been presented for plates supported along all four edges, and for outstands, wide columns, and deep beams. The results agree well with those of Yamaki for the fully supported plates, and are consistent with previous work on outstands and columns. For outstands, columns and beams, results for plates of both low and high aspect ratio, show clearly the transverse form of displacement present in each case. The degree of imperfection±sensitivity is seen to be greatest in columns, less in outstands, and less still in beams. For the latter two types, imperfection±sensitivity is greater in the shorter plates.

351

the strains mE at mid-thickness are given by the strain± displacement equations: 9 E ˆ m E ‡ zk k …a† > > > > > > where 8 > > 9 > > 1 2 2 > u;x ‡ 2 …w ;x ÿ w 0;x †; > > > < = > > > > 1 2 2 …b† mE ˆ > v;y ‡ 2 …w ;y ÿ w 0;y †; > > > = : ; u;y ‡ v;x ‡ w;x w;y ÿ w0;x w0;y ; …A:3† > > > and > > > > > > k T ˆ ÿ‰w ;xx w ;yy 2w ;xy Š …c† > > > > > where > > > > ; w ˆ w ÿ w0 …d† Equation (A.3b) above includes geometrical non-linearity. Eqs. (A.3a) and (c) correspond to classical plate theory (Love±Kirchho€ hypothesis), and include the usual assumption that the slopes w,x and w,y are small compared to unity. The stresses are elastic, isotropic stress: 8 9 2 < Ex = 1 1 Ey ˆ 4 ÿn :g ; E 0 xy

related to the strains by the usual stress±strain equations for plane ÿn 1 0

8 9 < sx = 0 sy ; 0 : ; txy 2…1 ‡ n†

…A:4†

where E is Young's modulus, and n is Poisson's ratio. The mid-thickness strains are similarly related to the stress resultants: 8 9 8 9 2 < Nx = 1 ÿn 0 < m Ex = 1 4 Ny : ÿn 1 0 ˆ …A:5† m Ey : g ; Et : ; Nxy 0 0 2…1 ‡ n† m xy

Appendix A

The stress resultants N are calculated from the force function F as follows, and this ensures that in-plane equilibrium is automatically satis®ed:

A.1. Appendix. Basic equations

NT ˆ ‰F;yy F;xx ÿ F;xy Š:

In the following, the comma subscript denotes partial di€erentiation once by each of the subscripts following the comma. Thus, for example:

Eqs. (A.3b), (A.5) and (A.6) are now combined to relate the displacements directly to the force function:   Et m Ex ˆ Et u;x ‡ 12 …w 2;x ÿ w 20;x † ˆ F;yy ÿ nF;xx ;   Et m Ey ˆ Et v;y ‡ 12 …w 2;y ÿ w 20;y † ˆ F;xx ÿ nF;yy ;   Et m gxy ˆ Et u;y ‡ v;x ‡ …w;x w;y ÿ w0;x w0;y †

v;x ˆ

@v @ 2 Mx and Mx;xx ˆ : @x @x 2

…A:1†

At any position (x, y) the following equilibrium equations must be satis®ed: 9 Rx ˆ 0 where Rx ˆ Nx;x ‡ Nxy;y ; …a† > > > > Ry ˆ 0 where Ry ˆ Ny;y ‡ Nxy;x ; …b† = Rz ˆ 0 where Rz ˆ Mx;xx ‡ 2Mxy;xy ‡ My;yy > > ‡ …Nx w;x ‡ Nxy w;y †;x > > ; ‡ …Ny w;y ‡ Nxy w;x †;y : …c† …A:2† At any position (x, y), the strains E at any level z, and

ˆ ÿ 2…1 ‡ n†F;xy :

…A:6†

…A:7†

By di€erentiation, eliminate u and v from Eq. (A.7) to obtain m Ex;yy

‡m Ey;xx ÿm gxy;xy ˆ N …w† ÿ N …w0 †

…A:8†

where the non-linear di€erential operator N … † is de®ned as follows: N … † ˆ … †2;xy ÿ … †;xx … †;yy :

…A:9†

352

G.H. Little / Computers and Structures 71 (1999) 333±352

From Eqs. (A.5), (A.6) and (A.8): 4

r F ˆ Et‰N …w† ÿ N …w0 †Š:

…A:10†

Eq. (A.10) is the von Karman [2]/Marguerre [3] compatibility equation for large de¯ections of isotropic plates.

References [1] Little GH. Ecient large de¯ection analysis of rectangular orthotropic plates by direct energy minimisation. Comput Struct 1987;26:871±84. [2] von KaÂrmaÂn T. Festigkeitsprobleme im Maschinenbau. EncyklopaÈdie der Mathematischen Wissenschaften 1910;4:349. [3] Marguerre K. The apparent width of the plate in compression. NACA TM No. 833 (1937). [4] Levy S. Bending of rectangular plates with large de¯ections. NACA report no. 737, 1942. [5] Little GH. An ecient computer program for the large de¯ection analysis of rectangular orthotropic plates. Comput Struct 1987;27:467±82. [6] Little GH. Large de¯ections of rectangular plates with transverse edges remaining straight. Comput. Struct. 1999;71(3):353±7.

[7] Yamaki N. Postbuckling behaviour of rectangular plates with small initial curvature loaded in edge compression. J Appl Mech 1959;26:407±14. [8] Ronalds BF. Local buckling strength of plate outstands. Int J Mech Sci 1990;32:725±44. [9] Ayrton WE, Perry J. On struts. Engr 1886;62:464, 465, 513, 514. [10] Timoshenko SP., Gere JM. Theory of elastic stability, 2nd ed. New York: McGraw-Hill, 1961. p. 362. [11] Rhodes J, Harvey JM. Plates in uniaxial compression with various support conditions at the unloaded boundaries. Int J Mech Sci 1971;13:787±802. [12] Rhodes J, Harvey JM. E€ects of eccentricity of load or compression on the buckling and post-buckling behaviour of ¯at plates. Int J Mech Sci 1971;13:867±79. [13] Rhodes J, Harvey JM, Fok WC. The load-carrying capacity of initially imperfect eccentrically loaded plates. Int J Mech Sci 1975;17:161±75. [14] Ronalds BF, Chapman JC. Elastic and ®rst yield analysis of hinged outstands under uniform compression. J Construct Steel Res 1991;20:191±220. [15] Louca LA, Harding JE. Torsional buckling of outstands in longitudinally sti€ened panels. Thin-Walled Struct 1996;24:211±29. [16] Dawe DJ, Wang S. Spline ®nite strip postbuckling analysis of thin rectangular laminated plates. Am Soc Mech Engineers 1996;73 (1):43±53 Petrol Div. [17] LUSAS 1998. A general purpose ®nite element system (version 12) User manual. FEA Ltd, Kingston, Surrey, UK.