Local and Overall Buckling of Thin-Walled Beams and Columns Using Finite Elements by Janine Diana Toneff B.Eng., McMaster University, 1982 A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies Department of Civil Engineering We accept this thesis as conforming to the required standard The University of British Columbia April 1986 © Janine D. Toneff 1986 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. I t i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of C I V I L , ^/ve?(N£a^N6> . The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date Apr! | ^ H8Q> ii Abstract Steel members with open thin-walled cross sections are used extensively in civil engineering structures. In addition to overall instability, (Euler buckling, torsional buckling, lateral-torsional buckling, etc.), the thin plates making up the cross section may themselves be susceptible to local plate buckling. The possible interaction of these two modes of buckling and its effect on member stability and strength is therefore of interest in the analysis of such members. The purpose of this thesis is to develop a tool for this type of analysis by adapting a finite element for a thin-walled beam-column of arbitrary cross section to include 'local' degrees of freedom to allow for cross section distortion. The formulation will involve geometric nonlinearities due to large displacements and rotations, but material behaviour will be limited to the linear elastic case. iii Contents Abstract ii Figures v Tables vi Acknowledgement vii 1 Background 1 1.1 Introduction 1 1.2 Literature Review 4 1.3 Theoretical Formulation 5 2 Theory 9 2.1 Incremental Procedure 9 2.2 Stress and Strain Tensors 11 2.3 Incremental Virtual Work Equation 12 2.4 Virtual Work for Thin-Walled Cross Sections . ' . . . . 14 2.5 General Method of Forumulating Stiffness Matrix . . . . 16 3 Formulation of Stiffness M a t r i x 10 3.1 Coordinate Systems 19 3.2 Incremental Virtual Work Equation 20 3.3 Generalized Coordinates 20 3.4 Stress Assumptions 24 3.5 Strain-Displacement Assumptions 25 3.6 Cross Section Displacement Assumptions 26 3.7 Stress Resultants 29 3.8 Incremental Forces 30 3.9 Shape Functions 31 3.10 Substitutions . 34 3.10.1 Linear Beam Term . 3 5 3.10.2 Linear Plate Term 36 3.10.3 Geometric Beam Term 37 3.10.4 Geometric Plate Term 38 3.10.5 Geometric Beam-Plate (coupled) Term 39 3.11 Equilibrium Equations 40 3.12 Element Matrices . 49 3.13 Transformation 57 4 Implementation 58 4.1 NISA's Thin-Walled Beam Element 58 4.2 Incorporating Superbeam 59 4.3 Eigenvalue Analysis 60 Contents iv 4.4 NISA 61 5 Results and Conclusions 66 5.1 Discussion of Results 66 5.1.1 Theoretical Calculations 66 5.1.2 Convergence 68 5.1.3 Angle Section 69 5.1.4 Wide Flange Section 71 5.1.5 Channel Section 75 5.1.6 Hat Section 76 5.2 Conclusions . 77 5.3 Further Study 78 References 80 Appendices 82 Al Symbols and Notation 82 A2 Component Matrices 84 V Figures Fig. 2.1 Rigid Body Movement 10 Fig. 3.1 Local Beam Coordinate System 19 Fig. 3.2 Local Plate Coordinate System 20 Fig. 3.3 Beam Degrees of Freedom 21 Fig. 3.4 Examples of Plate Nodes 23 Fig. 3.5 Plate Degrees of Freedom 23 Fig. 3.6 Rigid Body Beam Displacements 27 Fig. 3.7 Cancellation of Shear Terms 35 Fig. 3.8 Examples of Plate Angles, V»; 39 Fig. 3.9 Schematic Stiffness Matrix 50 Fig. 4.1 Buckling Solution 63 Fig. 4.2 Element Stiffness Routines 64 Fig. 5.1 Convergence of Plate Buckling Load . 68 Fig. 5.2 Plate Buckling Modes 69 Fig. 5.3 Angle Column Buckling 70 Fig. 5.4 Wide Flange Column Buckling (i) 71 Fig. 5.5 Wide Flange Column Buckling (ii) 72 Fig. 5.6 Wide Flange Lateral Buckling 73 Fig. 5.7 Wide Flange Local Buckling 74 Fig. 5.8 Channel Column Buckling 75 vi Tables Table 5.1 Hat Section Buckling 76 vii Acknowledgement The financial support of the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. Thanks are also due to Dr. Osterrieder for suggesting the topic and his initial guidance, Dr. Nathan for reading through the first draft of the theoretical formulation, and making suggestions, and Dr. Stiemer for additional financial support. CHAPTER 1 Background 1.1 I N T R O D U C T I O N For efficient design of structural members it is often desired to minimize the cross-sectional area of the member, while maximizing the radius of gyration (for columns) or moment of inertia (for beams). The intention is to economize on material and member weight while maintaining the desired strength. This reasoning has resulted in an increase in the use of members of thin-walled cross section in steel construction. Some examples of structural members in this category are: wide flange steel columns and beams, built up sections such as plate girders and box girders, hollow structural steel sections, angle sections (hot rolled) tubing and hat sections used in open web steel joists, channel sections for wall studs (cold rolled), etc. The behaviour of thin-walled structural members is not entirely the same as that of solid members. In pure flexure the same assumptions hold for both thin-walled and solid sections, i.e. plane sections remain plane and shear deformations are neglected. In torsion however, the case of thin-walled members becomes more complicated than that of solid sections. The torsional resistance of a thin-walled cross section is due to two different types of response; St. Venant torsion and warping torsion. In the case of solid sections warping deformations due to torsion are negligible and the torsional resistance is due only to the St. Venant term. The amount of torque taken by each of these two responses for 1 2 a thin-walled cross section depends on the shape of the section, the amount of restraint against warping and the type of torque applied. St. Venant torsion is due to "direct stress' or shear flow through a section. For a closed cross section this shear flow produces a moment about the shear centre and gives relatively high torsional resistance. For an open cross section the shear flow, which must still follow a closed path around the section, actually goes in both directions along each plate element. The resulting moment resistance is equivalent to a plate bending moment and much smaller than for a closed section. Warping torsion is a result of the stresses developed by the resistance of the cross section to warping. Both normal and shear stresses are developed along the warped section. Shear stresses due to warping contribute directly to a torsional moment about the shear centre, giving the warping torsional resistance. Normal stresses also have a component contributing to the moment about the shear centre since they are acting on the warped cross section. This second order effect contributing to the warping torsion is called the Wagner effect. Warping deformation of closed cross sections is usually negligible and these warping effects are ignored, but for open cross sections they can constitute the major part of the torsional resistance. Thin-walled beam-column members must often be analysed with respect to stability as well as strength. If a member is slender, overall buckling may be a concern. Overall member buckling may occur in different modes, depending on the cross section and type of loading, eg. simple Euler buckling, torsional or combined tosional-flexural modes under axial loads, lateral-torsional buckling due to bending, etc. The cross section of a thin-walled member may be said to be made up of a number of thin plates connected along their edges. These plate elements making up the cross section are susceptible to local plate buckling due to in-plane loading. This local buckling may be analysed independently of the overall member buckling. There may, however be some interaction between the local and overall member stability and strength and it would be useful to see whether including this interaction gives better results than separate analyses of the two buckling modes. 3 Much research has been done on the overall buckling of beams, including beams with thin-walled cross sections for which the cross section is assumed to remain rigid. There has also been extensive research done on plate buckling which can be applied to the local buckling of thin-walled cross sections. However, the interaction of the two buckling modes for thin-walled beam-columns has not been studied so extensively. The work that has been done in this area, deals with specific situations and cross sections only. Some generalized analysis of the interaction of local and overall buckling for an arbitrary cross section is still needed. A thin-walled beam element which has already been developed [1] will be the basis of the analysis, and local degrees of freedom will be added to allow local plate distortions. There will be coupling between the overall beam degrees of freedom and the local plate degrees of freedom which will allow observation of any interaction between the two buckling modes. The objective will be to discover under which conditions, if any, the interaction of these two buckling modes reduces the initial buckling load and/or ultimate strength and should be taken into account when analysing the stability of thin-walled members. The stability problem of thin-walled structural elements becomes quite com-plicated. Buckling may take place in the elastic or inelastic range of the material. There may be some yielding of the cross section before buckling occurs, especially with shorter member lengths, which reduces the member stability and changes the properties of the cross section. With steel sections there are usually residual stresses present in the cross section (due to welding, differential cooling in hot rolled members, cold bending in cold formed members, etc.). These residual stresses will have an effect on how and where plas-ticity develops through the cross section. Any strain hardening also affects the section response in the inelastic range. Imperfections, such as initial crookedness or dents, can reduce the ultimate strength of a short beam below its inelastic buckling strength. Initially some assumptions will be made to simplify the analysis with respect to these factors. 4 1.2 LITERATURE REVIEW The following is a brief overview of some of the work done on the topic of coupled local and overall buckling of thin-walled beam-columns to look at some differ-ent methods that have been used to approach the problem, and to provide a means of comparison and evaluation of the results of the present analysis. Only papers dealing specifically with the interaction of local and overall buckling of open thin-walled cross sections are considered here. More work has been done dealing with the two modes separately, (especially overall buckling of thin-walled cross sections), which is relevant to the general formulation presented here. References [18] to [21] give som general information on the analysis of thin-walled cross sectional members including buckling. The papers reviewed all deal with elastic deformations only. They ignore residual stresses and initial imperfections, unless otherwise stated. References [2] to [5] deal with local and overall buckling applied to I sections only. The work of Rajasekaran [2] was looked at in detail and in many aspects provided a basis for the formulation of the element used in the present work (i.e. the virtual work formulation using plate and beam degrees of freedom). References [3] and [4] both use the finite strip method for analysis. In reference [3] Hancock deals only with bending about the major axis and does not include warping. In reference [4] Hancock et al. include warping and some axial load - moment interaction. Akay et al., in reference [5], use the finite element method with plate elements for the web and one-dimensional elements for the flanges. They include analysis of a plate girder and a gabled frame. References [6] to [10] deal with channel cross sections only. Reference [6] analyses stiffened and unstiffened channels using the finite element method. The effective width concept is utilized, resulting in nonlinearities so that an iterative solution is used in the post local buckling range. Reference [7] deals with U-shaped beams used for composite box girders and gives results of tests performed on welded steel U-beams. The finite strip 5 method is used for analysis due to bending loads only. The effects of residual stresses and nonlinear material behaviour are included in this analysis. Tso and Ghobarah, in reference [8] model the channel cross section as a system of plates with springs at the plate intersections and use the matrix transfer method for their analysis. Torsional buckling is not considered. Reference [9] claims an exact analytical matrix solution which gives the interaction between local and overall modes, including warping and torsional effects. Stresses for this analysis must be constant in the longitudinal direction. Reis and Branco in reference [10] use a Galerkin solution for overall buckling and a Rayleigh-Ritz approach to the local buckling. Loading is combined bending and torsion. Reference [11] uses the finite strip method to determine the post locally buckled equilibrium path and performs tests and analysis on box columns. 1.3 T H E O R E T I C A L F O R M U L A T I O N Objectives There are two methods of approaching a stability problem. One is the eigen-value analysis, by which the critical buckling load is obtained in essentially one step. The other is to follow the load-deflection path as the load is increased to obtain the critical load. This second method is more laborious but it allows the complete load-deflection history to be studied and also allows observation of any post-buckling strength. Since post-buckling strength can be significant for thin-walled members and information about the total load-deflection path is useful, the load-deflection solution is intended rather than simply an eigenvalue analysis. The main objective of the theoretical formulation is to derive the stiffness matrix for a thin-walled beam element with local degrees of freedom. The stiffness matrix will be derived using both beam and plate degrees of freedom and will have both linear and geometric stiffness terms. Since a thin-walled open cross section is being considered, torsion and warping effects must be included. The formulation is kept as general as possible 6 so that any extensions to the analysis (eg. including inelastic behaviour) may be made easily in the future. General Assumptions and Formulation 1. Assume homogenous, isotropic material properties throughout the element and ignore any variations in properties due to residual stresses. 2. For the time being, all deformations are assumed to take place in the elastic range. When this is the case, the deformations due to local buckling and those due to overall buckling may be directly superimposed. 3. Most structural analyses are based on the assumption that deflections and rotations are small. When displacements are small, equilibrium calculations may be performed using the original geometry of the structure, i.e. displacements and rotations are small enough that the system of forces remains essentially the same in the deformed configuration as in the original. If, in addition, strains are assumed small, only the linear terms of the strain displacement relation are retained for a sufficiently accurate description of the behaviour. These two assumptions greatly simplify structural analysis. In buckling analysis, however, displacements can be large (even in the elastic range for slender members) and the effect of these displacements can no longer be ignored. Equilibium equations must be written with respect to the deformed configuration. This is then called a geometrically nonlinear problem (as opposed to a materially nonlinear problem resulting from a nonlinear stress-strain relationship). In the present analysis the strains are assumed small, but second order terms in the strain-displacement relationship will be retained to allow for large rotations. 4. The virtual work formulation will be used to derive the stiffness matrices for the finite element rather than minimization of potential energy. The potential energy formulation is valid only for the elastic range of deformations, since it depends on conser-vation of energy in the system. In anticipation of the possibility of extending the present formulation into the inelastic range, the more general virtual work formulation will be 7 used. The more common formulation of virtual displacements rather than virtual forces will be adopted. 5. Nonlinear analyses, in general, can be approached using incremental or it-erative solution procedures, or by combining the two methods. An incremental approach is employed by dividing the load-deflection path into a series of steps which are assumed small enough that the equilibrium equations may be linearized within each step. By following the load-deflection path in this manner, the configuration may be observed at intermediate points before failure. However, because of the linearization, an 'unbalanced load' will arise at each step. This error is removed at each step by iterating until the unbalanced load is negligible. With a purely incremental approach this error would accumulate, and although convergence is assured, the accuracy of the solution is not known. With a purely iterative approach convergence is not guaranteed, and information is only obtained for the total loading condition. In this study a combined incremental and iterative approach will be used. 6. An incremental virtual work equation will be derived. The approach taken to derive this equation depends partially on the reference system used to describe the deforming body. In the Lagrange description the deformed body is described as a function of a known previous configuration. In the finite element application, the Total Lagrangian formulation uses the original undeformed configuration as reference. The Up-dated Lagrangian formulation, refers the deformed body to the last known equilibrium configuration. In the Eulerian description the deformed body is described as a function of the current deformed configuration. With the Total Lagrangian formulation the elastic element stiffness matrix, must be augmented by two additional matrices; the geometric 'stress' stiffness matrix, due to the second order strain displacement terms, and the initial displacement matrix, due to large displacements. With the Updated Lagrangian formulation however, if the load increments are small enough, the initial displacement matrix can be neglected. The 8 updated Lagrangian approach will be adopted here since it avoids the calculation of an additional matrix at each step. 9 CHAPTER 2 Theory 2.1 I N C R E M E N T A L P R O C E D U R E Tensor Notation Tensor notation is commonly used to represent the components of stress and strain at a point in a body. The nine independent stress and strain components on any infinitesimal element are represented in tensor notation (in cartesian coordinates x, y, z) as follows: - strain tenson stress tenson Movement of a B o d y in Space Referring to figure 2.1, - Position 0 is the undeformed, original configuration. - Position 1 is the deformed configuration; the configuration after a number of displacement increments. This configuration is known from previous increments and is the reference po-sition for the Updated Lagrangian formulation. - Position 2 is the incremented configuration; an unknown configuration, a small increment away from configuration 1. 10 x 2 x 3 F I G U R E 2.1 Rigid Body Movement. Coordinate Systems Cartesian coordinates are shown as mz t- in general, for the coordinate in the » direction. The left superscript indicates the position of the reference coordinate system. This coordinate system can be set up in position 0, 1 or 2. Positions 0 and 1 are known, and coordinate systems set up in these positions are known as Lagrangian coordinates. Position 2 is unknown and a coordinate system set up in this position is known as an Eulerian coordinate system. Notat ion Tensors, coordinates and displacements will use the following notation: JJV.y - stress tensor *^c,-y - strain tensor ™Xi - coordinate in i direction JJ*u,- - displacement in the i direction The superscript is the current position (incremented position). This is usually the position at which quantities such as stress and strain are sought. The subscript is the reference position, ie. the position at which coordinate system is defined. 11 2.2 S T R E S S A N D S T R A I N T E N S O R S Green (Green-Lagrange) Strain Tensor: « 6 i 3 G - »'i + • i " " + n tt*»« tt*»y * - tensor defined w.r.t. the original (unstrained) configuration - symmetric - strain tensor in Lagrangian coordinates Almansi (Enler) Strain Tensor: - tensor defined w.r.t. in the deformed (strained) configuration - symmetric - strain tensor in Eulerian coordinates Cauchy (or Enler) Stress Tensor: {JJr,y -refers to in current configuration -gives 'true' current stress -associated with the Almansi strain tensor -symmetric Lagrange (1st KirchhofF-Piola) Stress Tensor: {|*Tt-y -refers to original undeformed configuration -unsymmetric Kirchhoff (2nd Kirchhofif-Piola) Stress Tensor: £V,y -refers to original undeformed configuration -associated with the Green strain tensor -symmetric -obtained by transformation from the Cauchy stress tensor using inverse deformation gra-o - n dients to make it symmetric; %<T}1 = ^ JJ,z/M- Qmxj)k ™rki 12 For the present analysis the best combination of coordinate system and ten-sors must be chosen to use with the Updated Lagrangian formulation of the incremental procedure. According to Fung [12], if the original state of a body is a "natural state, such as the homogeneous, stress-free state of an elastic body, the kinematic equations are simple ... in the Lagrangian description... whereas... for the Eulerian description ... the kinematic equations are complicated". Since we are only dealing with the static case for now, the Green strain tensor is more appropriate than the Almansi. Since the stresses are related to strains in the formulation, it is desirable to use a stress tensor which corresponds with the chosen strain tensor, in this case the Green strain tensor. Both Lagrange and Kirchhoff stress tensors refer to the original configuration, as does the Green strain tensor. However, the Kirchhoff tensor is symmetric and therefore more appropriate. Thus the Green strain tensor and Kirchhoff stress tensor will be used in the Updated Lagrangian formulation to derive the required stiffness matrices from the incremental virtual work equation. 2.3 INCREMENTAL VIRTUAL WORK EQUATION The principal of virtual work using virtual displacements, can be expressed as: 6Wi = 6We where SWi = J aSedV V SWe = J PSudS (2.2) S Here body forces are being ignored and P represents the surface tractions. Stresses will be represented by the second Kirchoff-Piola stress tensor, cr,y, and strains will be represented by Green's strain tensor, e,y = jfu^ y+uy,,- +uk„,uk,j). 13 For the incremental equation the total stresses and displacements (in position 2) must be expressed in terms of the quantities at the last known configuration (position 1) plus the incremental quantities. 2 i . . ff.y = <Tij + (7,-y where u,- and <r,y are the incremental quantities. U p d a t e d Lagrangian Approach lffij = + *»/ K = \ui + Ui The left subscripts indicate that the reference system for the updated formulation is defined in the last known configuration. This reference system moves with the deforming body so that the displacement w.r.t. this coordinate system is actually just the incremental displacement (ie. Ju,- = 0). There may, however be a stress at position 1 resulting from proceeding load increments, so that is not necessarily 0. Let the left subscripts be understood from now on in the Updated Lagrangian formulation and omit them from the notation so that we have 2u$- = ^ and V,y = ^.y + <7,y and the internal virtual work becomes: 62Wi = j2<Tij82eiidV v = J(V.y + <7,y) • 1/2 5(ti,-,y + U y „ - + U * , , - uktj ) dV V = 1/2 J l<Tij * ( t i , . , y + t t y „ - + « * „ • U j t . y ) dV + 1/2 J ^(5ti,-,y +Sujfi ) dV V V + l/2J&ij6{ukiiuk,i)dV v If the last term is neglected, the equation is linearized'within the increment. (This means that we will now be using a tangent stiffness for the increment and it is necessary to iterate 14 to obtain equilibrium in the unknown configuration). Ignoring body forces, the external virtual work can be expressed as 1/2 J 1 a s ; « ( « l „ u i , y ) ^ + l/2 J ^jSiti^+Uj^dV M j PSudS (2.3) V V s 2.4 V I R T U A L W O R K F O R T H I N - W A L L E D C R O S S S E C T I O N S To allow coupling between local plate buckling and overall beam buckling modes, dis-placements and stresses from both of these sources must be included in the virtual work equation. Since for the elastic case the principal of superposition is valid, the incremental displacements at any point on the cross section can be expressed as the superposition of the overall beam column and the local plate buckling displacement increments. The stress increments are also a superposition of those due to beam behaviour and those due to plate behaviour. S If position 1 is in equilibrium, then v s and the linearized incremental virtual work equation becomes Using Green's strain tensor and expressing the displacement and stress in-crements as sums of beam and plate displacements and stresses as follows: u t = uf + uf the virtual work equation becomes: V V (2.4) 5 15 With some manipulation of the dummy indices the first integral becomes /(*?•+*£•) •*(*.? +«.-.; ) ^ V and then expanding gives: J «,,? «,,f 5ti,,f + 0 $ . *u,,f) dV If plate displacement and stress increments are assumed to be due to plate flexure only, and the beam displacement and stress increments are assumed to be constant over the thickness of individual plates, the coupling terms in this integral disappear i.e. the integration, over the plate thickness, of the product of plate stress and beam strain, or of beam stress and plate strain, will be zero. Therefore the first integral is only: v Expanding the second integral of equation 2.4 and again using some manipulation of the dummy indices for the cross term gives: v The final form of the virtual work equation will have five terms contributing to the internal virtual work. These five terms will result in five terms in the element stiffness matrix as follows: 6{LB) = / of-Suifj dV gives linear beam stiffness V S(LP) = f crfjb'uif- dV gives linear plate stiffness v S(GB) = j / 1o-ij6(u/.,f uk*j ) dV gives geometric beam stiffness V 6(GP) = § / l^ijS{uk,f ukf) dV gives geometric plate stiffness V 6(GBP) = / l<T i j S ( u k , f ukl?)dV gives geometric coupled beam-plate stiffness V so that the virtual work equation can be written as 6{LB) + 6{LP) + 6(GB) + 6{GP) + 6{GBP) = 6We , (2.5) 16 2.5 GENERAL METHOD OF FORMULATING STIFFNESS MATRIX The following steps outline the basic formulation to be used to obtain the stiffness matrices for a thin-walled beam element with geometric nonlinearities. (1) Incremental Virtual Work Equation SWi = 8We See updated Lagrangian formulation and equation 2.5 for terms included in this specific case. (2) Generalized Coordinates/Degrees of Freedom Choose the degrees of freedom appropriate to describe the overall beam de-formations and the local plate deformations. Let the vector of these generalized coordinates be s = Bp + Bp (including both beam and plate d.o.f.) (3) Beam and Plate Displacement Assumptions The displacements of any point on the beam cross section must be expressed in terms of the generalized coordinates. Let I = J s, where § is the displacement of any point on the section and J is a matrix defining the geometric relation between the displacement of this point and the displacements and rotations given by the generalized coordinates. (4) Finite Element Displacement Assumptions/ Shape Functions The shape functions will give the displacement of any point in the element in terms of the nodal d.o.f. Let u = A a, where u is a vector of displacements at any point in the element and A is a matrix of shape functions describing the deformations (5) Strain Displacement Relation/ Kinematic Equations Strains are given as derivatives of displacements, e = Ln where c is the vector of strains at any point in the element and L is a matrix of derivatives; 1st derivatives are used to give axial strains, and 2nd derivatives are used for bending strains (curvatures). 17 (6) Material Law/ Constitutive Relation Using a linear elastic material law, the expressions for incremental and total stresses respectively are: a = Ee, and a = Ee. (7) Stress Resultants/ Generalized Forces Stress resultants may be expressed by the general form F = JA Bo dA, where a is a general stress and B is a geometric factor allowing bending and torsional moments to be included as generalized forces. Incremental stresses have their corresponding incre-mental forces. (8) Combining Steps (1) to (7) (a) Linear Terms Let the linear internal virtual work be 6WiL and the linear virtual strain terms be Sei. MiL = J 6eLTB&dV v -JSeLT(J B&dA)dx = ^6eLTtdx using the following substitutions 6e = LSn = LAJcSs P = C s where C is the appropriate stiffness for axial, bending or torsional d.o.f. EA, EI, or JG respectively, yields 6WiL = J^a{AJ)TLTCadx if P is a vector of incremental surface tractions, the external virtual work, 6We is given by <5Bt P and the virtual work equation becomes: 6aT{J{A3)TtrCdx)a = 6aT~P 18 Cancelling the incremental virtual displacements, gives the form K ^ s = P , so that we get the linear element stiffness matrix, KL = J{AJ)TLTCdx (2.6) (b) Geometric Terms Let the geometric internal virtual work be 6WiG and the geometric virtual strain terms be given by 8kG. 6WiG = J 6eGT<rdV v = 16€'GT(J B<rdA)dx = J^eGT¥dx F is the vector of (known) total forces. The geometric stiffness comes from the second order terms in the strain tensor and these can be represented as follows: to = (£«)2 = u r L r L u = s r ( A J ) r L r L ( A J ) s with this substitution and the same expression for the external virtual work as before, the virtual work equation for the geometric terms becomes: 6aT{ j f ( A J ) T L r F L ( A J ) dx) S = SBTP Again cancelling virtual displacements gives the form K ^ s = P and the geometric element stiffness, K g , is given by KG = ^ ( A J ) r L T F L ( A J ) dx 19 CHAPTER 3 Formulation of Stiffness Mat r ix This chapter will go through the steps given in section 2.5 in more detail to formulate the stiffness matrix for a thin-walled beam element with local distortion capability. 3.1 COORDINATE SYSTEMS L o c a l B e a m Coordinate System y,v i z, u z,w \ [ I FIGURE 3.1 Local Beam Coordinate System. • I = element length - x, y, z are cartesian coordinates (right handed system) - u, w, it; are displacements in the directions of z, y, z respectively - dimensionless longitudinal coordinate = A = x/l - beam coordinate system has origin at the centroid of beam node i; node j is 20 Local Plate Coordinate System a,r 7,< 1 / 1 F I G U R E 3.2 LocaJ Plate Coordinate System. - b = plate width - I — plate length = element length - a,/3,7 are cartesian coordinates (right handed system) - r,s,t are displacements in the directions of a,/?, 7, respectively - a is in the same direction as the longitudinal beam coordinate, x - £ = (5/b is the dimensionless coordinate in the transverse plate direction (/? direction) 3.2 INCREMENTAL VIRTUAL WORK EQUATION From section 2.4 the incremental virtual work equation is of the form: 6{LB) + 6{LP) + 6{GB) + 6(GP) + S[GBP) = SWe where the left side of the equation represents the internal virtual work from the linear and geometric, beam, plate and coupled beam-plate terms. 3.3 GENERALIZED COORDINATES B e a m Degrees of Freedom Thin-walled beam finite elements have already been developed by a number of researchers ([2], [6], etc.) . They almost all include the following seven generalized coordinates which 21 will describe translation in the 3 coordinate directions, bending about the 2 transverse beam axes, twisting about the longitudinal beam axis and warping of the cross section. Axial translation is referred to the centroid but all other cooordinates are referred to the shear centre. The 7 degrees of freedom for each beam node shown in figure 3.3 are: 1. u; axial translation of centroid 2. w; translation of shear centre in the y direction 3. t/; slope in the y direction (rotation about z) 4. to; translation of shear centre in the z direction 5. u/| slope in the z direction (rotation about y) 6. 9] rotation about the longitudinal beam axis through shear centre 7. 0'; rate of change of 0 w.r.t. the longitudinal axis. FIGURE 3.3 Beam Degrees of Freedom. For small incremental displacements v and w, v' = || and w' = — ^ ~ so that u/ is in the positive y direction. Plate Degrees of Freedom Not very much work has been done regarding the local deformations of thin-walled beams and the choice of nodal variables to describe these deformations is not so clear. 22 It is generally observed that when local buckling occurs in a short section of a thin-walled beam, the corners remain straight, ie the corners between plates making up the cross section do not buckle themselves. The corners may rotate without a change of the the angle between the adjoining plates [11]. These observations apply to most practical cases, except where the cross section is made up of a large number of plates and/or the angles between the plates are large so that the section approaches a shell type of structure. In these cases there may be significant changes in the corner angles when the section buckles. As a result of these observations, corner rotations are chosen as generalized coordinates to describe local buckling deformations. The corner angle is assumed to remain rigid during the deformation. We can define 'edge plates' as plates which are connected to the cross section along only one edge and have one free edge (eg. simple stiffeners, flanges of I sections). It will be assumed that these edge plates rotate rigidly about their connecting edges, ie. the edge plates remain straight in the transverse direction at any cross section). This has been observed in practice, and therefore rotations do not need to be specified at free edges. 'Internal' plates are then defined as plates with 2 longitudinal edges sup-ported by other plates. These plates will have corner rotations specified at both of their edges to determine the deformation across their width. The nodes defining the cross section will be called 'plate' nodes (as opposed to 'beam' nodes which define the 2 ends of the beam elements). 'Internal' plate nodes can now be defined as plate nodes connecting more than one plate, and 'edge' nodes as plate nodes touching only one plate. The number of plate nodes depends on the cross section being analysed and the detail required (it may be possible to obtain a more accurate description of deforma-tions if more nodes are added on the straight parts of the plates). The corner rotations will describe the plate deformations in the transverse direction at any cross section. However, the deformations are localized and therefore the rotations used to describe these deforma-tions willvary along the length of the element. To describe this variation along the length 23 (a) 1,2 internal nodes 3-6 external nodes a internal plate b-e edge plates (b) 1-6 internal nodes 7,8 external nodes a-e internal plates f,g edge plates F I G U R E 3.4 Examples of Plate Nodes. of the element, the derivative of the comer rotations w.r.t the longitudinal coordinate will be included as additional plate degrees of freedom. Therefore at each internal plate node there will be 2 plate d.o.f., 6 and 6\ where 6 is the rigid rotation of the node and 6' is its derivative w.r.t the longitudinal coordinate. In the program, internal nodes must be numbered first, then edge nodes. ft K 4>i 6i! 4i <& F I G U R E 3.6 Plate Degrees of Freedom. In figure 3.5 then, there are 8 local generalized coordinates for (internal) plate A and 4 generalized coordinates for (edge) plate B. 24 Element Generalised Coordinates The total number of generalized coordinates per element is 14 + 4(n) if there are 14 beam d.o.f. and n internal plate nodes in the cross section. The element d.o.f. can be set up in vector form as follows: v = w = e = u = teu/)0. {fa 4>i 4>j <t>))m Here m = 1, n and primes imply derivatives w.r.t. the longitudinal coordinate. 3.4 STRESS ASSUMPTIONS Beam Stresses From basic theory of slender beams the assumptions are: 1. The cross section remains rigid (see note in plate stress assumptions). Plane sections do not necessarily remain plane since warping of the section is being included. 2. Fibres remain normal to the neutral axis after deformation; shear deforma-tions are ignored. 3. Poisson's effect is ignored. With these assumptions the only stress acting for pure beam bending is axx. However, in order to satisfy equilibrium, there must be shear stresses associated with flexure; rxy from bending about the z axis and TXX from bending about the y axis. All normal stresses acting perpendicular to the beam axis are ignored, as are the shear stresses in the y-z plane. Therefore the applicable stress tensor for general beam bending, using the local beam coordinates of section 3.1 is: 25 Plate Stresses From bending theory for thin plates, we have the assumptions: 1. Initially straight lines normal to the middle surface remain straight and nor-mal to the middle surface during deformation. 2. Stresses normal to the surface of the plate are ignored. With these assumptions and using the local plate coordinate system of section 3.1 the applicable plate stress tensor is: V o o o 3.5 STRAIN-DISPLACEMENT ASSUMPTIONS B e a m Strains The strains required in the virtual work equation will be only those corre-sponding to the non-zero beam stresses given in the previous section, that is exx, i x y, and Ixz' For slender beams it can be assumed that axial strains remain small, even for large out-of-plane displacements. It can also be assumed that variations of axial dis-placement across the section are small. Therefore all products of derivatives of axial displacements, u, may be neglected. Using Green's strain tensor, we are left with the following strain-displacement expressions: exx «tt„+-(t/lx2 + u),x2) Ixy = lyx « «!» +v>x +*>i* » i f + f « Ixz = Ixx « «iz V,z +*>>x W,z Plate Strains Only the strains saat fap and e^, corresponding to the non-zero plate stresses of the last section need to be considered. Strains in the middle surface are assumed 26 small compared to bending strains. This means that, although we are taking geometric nonlinearities into account in the beam deformations, out-of-plane plate displacements are still assumed to be small. Therefore the products of in-plane displacement gradients may be neglected compared to those out of the plane. Using Green's strain tensor, the relevent plate strain-displacement expressions are: eaa w r i o +2 3.6 CROSS SECTION DISPLACEMENT ASSUMPTIONS It is necessary to express the displacement of any point on the thin-walled cross section in terms of the generalized coordinates. In doing this, some simplifying assumptions are made about the deformations. For the updated Lagrangian formulation, it is assumed that incremental rotations are small enough that expressions involving this quantity may • • • • be linearized, i.e. the incremental rotation, 9, is small enough that sin0 « 9 and cos0 « 1. The beam generalized coordinates describe the rigid body motion and warp-ing of the cross section. The movement of any point on the cross section can be described in terms of these coordinates. In addition, the plate generalized coordinates describe flex-ural deformation of the plates making up the cross section. The movement of any point through the thickness of the plate may be described as a function of these plate generalized coordinates. Note that there is no translation of internal plate nodes with respect to each other, and in this sense the section still behaves as 'rigid' body. 27 F I G U R E 3.6 Rigid Body Beam Displacements. B e a m Displacements For a general cross section the following quantities are defined: - O is the origin of local beam coordinate system which is also the centroid of the section - S is the shear centre or centre of rotation - vs and u>5 are rigid body translations of the shear centre - $s is the rigid body rotation about S (positive counter-clockwise) - y and z are local coordinates of any point P on the cross section; these coordinates remain constant during rigid body beam deformations. - ey and ez are eccentricities of the shear centre from the centroid; these also remain con-stant during beam deformations. The following expressions derived from the geometry of rigid translation and rotation of the section, can be used to describe the movement of any point, P, in the y and z directions, due to beam bending: vp = vs - (z - et) sin 9S - {y - e,)(l - cos 6S) wp = ws - (y - ep) sin 9s-{z- ex){\ - cos 9S) 28 Linearizing w.r.t. 8$ yields the commonly used expressions: VP « «S - (* - cx)8s (3-2) u/p « u;5 - (y - e„)05 (3.3) The translation of P in the longitudinal direction arises from 4 different sources: 1. rigid body translation in the axial direction 2. bending about the y axis 3. bending about the z axis 4. warping The contribution from each of these sources is as follows: 1. axial translation of the centroid = u0 2. For plane sections the longitudinal displacement at P, due to bending about z, is proportional to the distance of P from the neutral axis in the y direction. For a small incremental rotation, |j, the displacement of P in the axial direction is — y vtx 3. Similarly, for bending about the y axis, a small rotation ^ produces an axial displacement —z w,x 4. The axial displacement due to warping of a thin-walled cross open section is given by 0s'un, where 6S' is the twist and w„ is the normalized sectorial coordinate (see appendix 1 for definition). The total displacement in the x direction from these sources is then, ;• •• ^••^•^->.t '^.^*^+.w,iV (3-4) Equations 3.2 to 3.4 give the displacement of any point P on the cross section in terms of the 7 generalized beam coordinates. Plate Displacements 29 Local displacements of points on the plates making up the cross section must be expressed in terms of the plate generalized coordinates, ttm. When suitable shape func-tions are chosen, these #m will define the out-of-plane displacement, t, which will in turn define the in-plane plate displacements due to bending (since in-plane plate displacements due to pure flexure are functions of the out-of-plane displacement derivatives). The dis-placements in the three plate coordinate directions are: r=-itia (3.5) a = ~lt,p (3.6) t = t The derivatives of equations 3.2 to 3.6 will be used in the virtual work expressions for the strain terms. 3.7 STRESS RESULTANTS From basic strength of materials 'stress resultants' or 'generalized forces' can be defined by integrating the various stresses over the cross sectional area. Beam Stress Resultants For beam stress resultants the integration is performed over the entire cross-sectional area. P = J axx dA axial force A Vy = J rxy dA shear force in y direction A Vt = / Txz dA shear force in z direction A. My = /oxxzdA bending moment about y A Mx — J<rxxydA bending moment about z A Tw = f(Txg(y — ey) — Txy{z — et)) dA torsional moment due to warping A Mw = JA oxxu dA bimoment KP = IA a*x((y ~ e*)2 + (2 ~ O 2 ) d A Wagner coefficient 30 Notes: 1. the torsional moment, Ms has 3 components; Mx = Tw + Tsv + T9. - Tw is the moment due to warping shear stresses (= —EIU6'S') - Tsv 13 the St. Venant moment due to non-uniform shear stresses on the cross section (= JG9'S) - is the moment due to the second order Wagner effect. (Normal stresses due to warping have a component in the y-z plane if the cross section is warped) (= Kp9s) 2. The St. Venant moment for an open thin-walled section is due to the twisting moments along the plate edges, resulting from the variation in shear across the thickness of the plate. If beam stresses are assumed to be constant through the thickness of the plates making up the cross section, there can be no St. Venant torsional moment. The St. Venant torsional moment is therefore introduced artificially to the virtual work equation later in the formulation. 3. Using incremental rather than total stresses gives the corresponding incremental stress resultants. Plate Stress Resultants For plate stress resultants the integration is performed over the plate thick-ness and the resultants are expressed per unit width of the plate. A/2 Ma = J aaa • 7 <f7 bending moment about -A/2 A/2 = / • 7 dr\ bending moment about a -A/2 A/2 Map — / rap • 7 ^7 twisting moment -A/2 where h is the plate thickness and the incremental stresses will again give corresponding incremental stress resultants. 3.8 INCREMENTAL FORCES The geometric stiffness terms were derived using total stresses and therefore total forces, 31 which are known from the previous increments. The linear stiffness terms contain in-cremental forces which are unknown. Constitutive relations are used to express these incremental forces in terms of geometric section properties, Young's modulus and the in-cremental strains. From elementary beam and plate theory the incremental forces may be ex-pressed as follows: P = EAu'Q M„ = -EIyws" Mz = -EIzvs" ' II Mu = —EIw9g • i Tsv — J G^s Ma = -D(t,aa+vttf)$) Mfi = -D(i4)aa+t,^) Ma? = -D(l - u)t1aP 3.9 SHAPE FUNCTIONS In section 3.3 the generalized coordinates were chosen for both the beam and the plate deformations. 'Shape functions' in terms of these generalized coordinates, along with boundary conditions, will define uniquely the deformed shape of the element. The choice of shape functions is the next step in the formulation. Polynomials are one common choice of shape function for use in finite element analysis, since they can be made to closely fit various curves, they are easily differentiable, and convenient in programming. Polynomials will be used in the present formulation to describe both beam and plate deformations. B e a m Displacement Shape Functions A common displacement shape function used for beam bending elements are the (cubic) Hermitian polynomials, which describe the shape of the beam in pure bending 32 in terms of the end displacements and rotations. Since (transverse) displacement and slope have already been chosen as nodal variables at the two beam nodes, we have the 4 constants necessary to determine uniquely the cubic polynomial and these Hermitian polynomials can be used as shape functions. They may also be used to describe the variation of rotation about the longitudinal axis, 0$, since we also defined 9g and its derivative at each node as beam degrees of freedom. node and therefore only a linear shape function can be defined uniquely. Since axial displacements are small, a linear polynomial should be sufficient to describe the variation in axial displacement along the element. The Hermitian shape function can be represented in vector form, in terms of the dimensionless longitudinal coordinate, A, as follows: The generalized coordinates specified in section 3.3, along with these shape functions, sat-at least, the finite element solution will converge to the exact theoretical solution as the element mesh is refined. Plate Displacement Shape Functions At any beam node we have the plate node rotations and the rate of change of these rotations in the axial direction as the plate generalized coordinates. In the longi-tudinal direction the variation of plate nodal rotations may be described using Hermitian polynomials in the same way in which the variation of 0§ was represented. In the axial direction we have specified only the displacement at each element and the linear shape function can be written as: isfy the convergence requirements for a finite element. That is, for the beam deformations 33 In the transverse direction we have two nodal variables to define the shape of an internal plate, and one for an external plate. It is convenient to use polynomials as shape functions to describe plate de-formations, to be consistent with the beam description. Since there is no translation of internal plate nodes relative to each other, and a cubic type of curve is appropriate to describe the buckled shape of an internal plate, Hennitian polynomials corresponding to rotations will be used to describe the out-of-plane plate deformation in the transverse di-rection. Since edge plates rotate rigidly, their displacement can be described in terms of a linear shape function in the transverse direction. The plate shape functions for internal plates, in terms of the dimensionless transverse plate coordinate, £, in vector form, are: and for external plates: {/!> = {« Therefore <f>m, the rotation of node m at any point along the element length is given as a cubic function of A and then the transverse plate displacement, t, is given as a cubic or linear function of this <f>m, for internal or external plates, respectively. The plate deflections can then be described in total as follows: if #m and * n are the vectors of plate generalized coordinates for two adjacent plate nodes, m and n, then the rotation of node m at any point along the element length will be given by <f>m = (7*3){*m}; similarity for node n, <f>n = (n3){<&„}. Then the actual out-of-plane plate displacements for an external plate at node m will be: «c = ( / i ) ^ = (/i){n 3){» m} = <^i){»-> and for an internal plate between nodes m and n: 34 The displacements at any point on the element can then be expressed in terms of the shape functions and generalized coordinates as follows: beam displacements: «o = K H * o > «S = M R } W S = K K w J h = (n3){e8} out-of-plane plate displacements: * e = ( * l > { * m } 3.10 SUBSTITUTIONS Now that the background information and assumptions have been given, the stiffness matrix for a thin-walled beam element, including local deformations can be derived using the steps of the general formulation of section 2.5. Substitutions will be made in each term of the virtual work equation (equation 2.5) individually. The general procedure for each of these terms will be as follows. Given the generalized coordinates of section 3.3, substitute first the results of sections 3.4 and 3.5, then the beam and plate displacements from section 3.6. Next insert the stress resultants of section 3.7, and for the linear terms, substitute again for the incremental forces given in section 3.8. Finally insert the shape functions for all of the terms. Note: All displacement quantities are understood to be incremental from now on, even if the dots are omitted. All forces and stresses are assumed to be total unless dots are shown. 35 3.10.1 Linear Beam Term 6(LB) = J afyuij dV v V = J ~ ™s" - Ws" + u*"s) dV+J rxyS{-vs' + u,1y8's + vs' - {z - et)8's) dV V V + J Tx,6{-WS' + u„ 8' + ws> + (y -t9)8s) dV v = f {P6uQ' - My6ws" - Mz6vs" + MJ0"s - TJ8'S) dx JQ The shear stress terms cancel, effectively eliminating the warping torsional moment when the following observations are made: z j y FIGURE 3.7 Cancelation of Shear Terms. dutn dun ds . . - -rama — -(y - ey) dz da dz 36 Also introduce a term JQ Tsv69ls d x *° mclude the St Venant torsional mo-ment. 6{LB) = ^ {EAUQ'SUQ' + EIvws"6wsn + EIzvs"6vs" + EIJS66S + JG9'S69'S) dx ^{5n 0 ) {n 1 ' } (n 1 ' ) {n 0 } + - / { « 5 w 8 ) { n 3 " } ( n 3 " ) { w 8 } + ^(5v 8 ){n 3 "}(n 3 "){v 8 } + ^(5e8){n3"}(n3"){eg} + ™(sea){n3'}M{&a}y\ 3.10.2 Linear Plate T e r m 6(LP) = J a£6«,,f dV v = J ea6r,a+01)88,0+Taf}6(r,p+8,a)dV v V = - J Ma6t)aa +Mp6t,pii +2Maf)6ttafi dA A = D J t,aa 6t,aa + v t m Stiaa +vt,aa 6t,fip +t,fifi Btm +2(1 - u)tmp 6t,afi dA A This intergral must be performed for each plate making up the cross section, so that this term actually becomes a sum over the N plates in the cross section. The integration over the plate area can be replaced by a double intergral over the plate width and length. N 1 b 6(LP) = Y,D J f (t,aa 6t,aa +t,fifi6t,fifi + v{t,aa 6t,fifi +tififi 6t,aa) , = 1 o o + 2{l-v)t,af)6t,ali)dl3da It has been shown that there will be different shape functions applied to internal and external plates. For brevity of notation, use the general shape function 'F' to apply to all plate deformations for the time being, until the actual stiffness matrices are calculated, at which time the internal and external plate shape functions will have to be dealt with 37 separately. Then the linear plate term can be expressed as: + K i -lt, 3.10.3 Geometric B e a m T e r m 6{GB) = llllvloij6{xiJuk*)dV = 1/2 (Txx6{v,l +w,2x ) d V + Jv rxyS{v,x v,y +w,x w,y) dV + / rxz6{v,xv,z+wixw,z)dV Jv = 1/2 J axx6[(vs' - { z - e2)B's? + {ws' + (y-ey)6°s)2} dV + J rxyS[(ws' + (y - ey)9's)9s\ dV + J Txz6{(Vs>-(z-ez)9's)(-9s)) dV = 1/2 J^(Txx6\vs'2 + ws'2 - 2zvs'B's + 2ezvs'9's + 2yws'9's - 2eyws'9's + (z- ez)2e's2 + (y- ey)26's2} dV + J Txy6[ws'9s + (y - ey)9's9s] dV + J rxzS[-vs'9s + (z- ez)9's9s] dV = 1/2 / P6[vs'2 + tos'2 + 2ezvs'9's - 2eyws°6's\ dx Jo + / \-My6(vs%) + Mz6{ws%) + l/2KpS(9's2)] dx Jo + f Vy6(ws'9s)-VzS{vs'9s)dx Jo + fv\rxy{y-ey) + rxz{z - ez)]6(9s9°s) dV The last term can be neglected since it is of higher order. Expanding the second order virtual terms and then substituting shape functions yields: 6{GB) = j f 1 (y{*v8){n3'}(n3'){v8} + y(Sw8){n3'}(n3'){w8} + ^{5e 8 ){n 3 ' }(n 3 ' ){v 8 } + £^(5v8){n3'}(n3'){e8}' 38 - ^(5G s ) {n 3 ' } (n 3 ' ) {w 8 } - (^5w8){n3'}(n3'){es}) dA - [ ( ^ { 6 * s ) M M M + ^ < K > W}(n3'){e8}) dA + fQ (^(*e8){»3'}(n3'){w8} + (^5ws){n3'}(n3'){e5}) dA + J* (vff(5ea){n3}(n3'){w8} + Vy(Sw8){n3'}{n3){Ba} ~ V,(5e8){n3}(n3'){v8} - V,{5v8){n3'}(n3){08}) dA 3.10.4 Geometric Plate T e r m 6(GP) = 1/2 j lOij6{uk{ u*,f) dV = 1/2 / <Txx6{tJ)dV Jv Again the total stiffness must have a contribution from all of the plates so that 8{GP) is represented as a sum over the N plates in the cross section. The volume integral can be represented as an area integral times the plate thickness (assuming constant plate -thickness). The beam stress is the total stress, due to axial force bending moments and warping, up to the last known equlibrium configuration. The geometric plate term then becomes 39 3.10.5 Geometric Beam-Plate (coupled) T e r m 6{GBP) = J V^tt*/ uk,p) dV v = / * « * ( « * „ % , / ) dV v The term tik,xBuk,ap represents the coupling of beam and plate displacements. The re-peated k subscript implies that the two displacements must be in the same direction. Neglecting again the product of axial displacements and retaining only out-of-plane plate displacement derivatives leaves u k , x B u k , a p = v,x ty,xp+u>B tz,xp, where ty and tg are the components of the out-of-plane plate displacements in the local beam y and z directions. For a plate in the cross section with arbitrary orientation, the components of its normal displacement, X, can be expressed as ty — tsin^ and tz = —tcos^, where rjt is the angle from the positive beam y axis to the plate at the lower plate node number (internal nodes numbered first), measured positive counter-clockwise. plate i> a 270 -t 0 b 180 0 t c 0 0 -t d 180 0 t e 0 0 -t plate V» a 60 0.866t -0.5t b 0 0 -t c 300 -0.866t -0.5t d 180 0 t e 0 0 -t F I G U R E 3.8 Examples of Plate Angles, t/>(. Summing over the plates in the cross section and expressing the volume as the plate area times the thickness, gives 1/ . A >B ty,p + w B tz,p)d/3dx 40 (v>S +{y~ cos fcl dfi dx + {5v8){n3'}(F'){*i} ~ (z-et)[6*t{F'}M{ea} - {z-e,){«0B){iis,>(^{»i>] + cos^[(5*i){n(n3'){w8} + (5w8){n3'}{F'){*i> + (y - e,)(0»l){l*}{n,'){O>} + (y - e„){*e.){n3'}(^ ){*i}] d^A 3.11 EQUILIBRIUM EQUATIONS Gathering terms according to their virtual displacements yields five equilibrium equations. The right side of these equations is the external virtual work, ie. the virtual displacements multiplied by their corresponding real nodal forces. When the virtual displacements are effectively cancelled, we are left with equilibrium equations of the form K A = P, where K is the stiffness, A is the incremental displacement and P is the real nodal (generalized) force. The five equilibrium equations before cancelling virtual displacements are shown below. The (L' subscripts indicate terms contributing to the linear stiffness matrix. All other terms contribute to the geometric stiffness. (1) (2) 41 + (6w3)(V;Mi/lViMi/l) (3) + [ P 6 Z 7 ^ fQlMM dxj {v8} - p * Z ^ i j\n3'}M <*A] { W 8 } + f * M M dx] {e8> + f\n3}M <*A] {W8} - [v, j\nz)M dX [cosfc(y - ef) jf 1 j\n3'}(F>) d£ dA]) {*;}} {v8} } = {6BB)(Mi Mt/l Mi Mi/l)J (4) < * » i > / f l jf + jf jf { ^ } < F , ^ « i A +fer/oi{^};F">rfedA]+fe/70i{F">^)rf^A + [ 2 ( i - v ) ^ r / o l { ^ } ( ^ > d ^ A ] ) x w -8in^(z-ez)f f {F'}{n3V^A{e8} + cos^ / / {F'}{n3V^A{ws} ./(WO «/0 •'0 + 42 + cos V,(y - e9) J" f \ F ' } M d(dX{&a} = MiJl M>px Mplf (5) For the geometric stiffness terms, the generalized forces must be expressed in terms of the nodal forces. From the displacement assumptions and basic strength of materials relations we can obtain the variation of these generalized forces along the length of the element in terms of the nodal values (ie. at x = 0 and at x = /). A linear shape function was assumed for u. Since P = E A ^ , P is constant along the length of the element. A cubic polynomial shape function was assumed for the transverse beam displacements, v and w, and the rotation, 9. Since the bending moments and bimoment are functions of the second derivatives of these quantities, (Afff = —EIyw", Mz = —EIzv"f and Afw = —2?/u0"), they can be taken to vary linearly between the two beam ends. Using the linear shape function in terms of the longitudinal coordinate, A, the expressions for the bending moment and bimoment along the element are as follows: Ms = (1 - A)M„(0) + A Ay/) Mz = (l-A)M,(0) + AAf,(/) A/u> = (l-A)A/w(0) + AMu;(0 The Wagner coefficient, Kp, also appears in the geometric beam term and must be ex-pressed in terms of nodal values. The normal stress on the beam section due to warping is given by a™ = Eun9" and the Wagner coefficient is defined in Appendix 1. Therefore, if 9 is a cubic function, Kp will be linear and its variation can be expressed as: Kp = (l-\)Kp(0) + \Kp(l) The shear forces required for rotational equilibrium can be expressed in terms of the moments at the nodes. v _ Mz(l) - Mz(0) 43 M„(Q - M,(0) / The five equilibrium equations are revised by cancelling the virtual displacements and substituting the above expressions for Myi Mx, Afw, Kp, Vy, and Vt ( let My(0) = M ,^ My(l) = My, etc.) Also now separate the internal and edge plates with their dif-ferent shape functions F3 and 2^ . The summation is divided into two parts; one for the external plates (represented as J^ ) and one for the internal plates (represented as £ ) . ep ip *j- ['{«/}(»/) dX{n0} = (P° P»)a (1) Pe -M° r1 M* - M* fl + z l y I MMdX{e^-^-j^JQ x{n3'}(n3')dx{es} Ml-M?. fl j—2. jf {n3'}(n3)dA{ej v + , -((Ml - M°y)j + {Ml - M?)j- + (JU* - jr) ^ sin A £ f\ A { « , ' } { 0 ^ ^{SJ + ((A*J - AfJ) j - + {Ml - M?)± + (Mi - M3)£j ^ sin * J J A{n3'}(F3') ^ d\{*$ = (V; A^/lViMl/lf (2) 2?/ r"1 P f 1 y y o K } ( « . " ) 4 A { w , } + ^ y o { n ^ ^ v A i w j _Pey-M° J * M M d X { Q a } + M^-Ml j \ M M d X { B n } + ^^JQlM(n3)dx{ea} 44 = {V?M9/lV>Mtv/l)T (3) ^ f\nz"}(n3") dX{Ba} + ^ f\n3')M dX{08} + S JQlMM dX{ea} + jf1 A{n3'}<n3'} <*A{08} Pe -M° rl M1 - r1 + ^ - r ^ Jf {n3'}<»3') <fA{v8} - ^ L _ J ^ Jf A{n3'}(n3') dA{v8> - JT (n3'}(n3') <*A{w8} + ^ L_^L jT A{n3'}(n3') dA{w8> + ^^/ o 1 {n 3 }(n 3 VA{w 8 } + ((AiJ - AiJ)^ + (JilJ - Ai?)^ + (Mi - Af^)£) [- sin - e.) jf' jf' A{n3'}(F/) dfilA{» i} + cos UV ~ e,) £ ^ A W X O <*A{*i}] } + ? ^ { g * f + ¥ + ¥ ) [ - * « - - o j [ , j [ W ' . - e « ( « s ) / 1 j [ 1 { B » ' > W ) ^ ^ { » i } ] + cos 45 + cos * ( y - ey) jf * jf' A{n3'}(F3') ^ dX{9{^ } = ( A ^ Ai2// Ml Ml/if (4) + 5 ll jf{*'«}{Fi>,) didX+k I! f ) + fjn* / ' / V i W > ^ « t t { * . ) - " i n f c(*-0 /' f" {F^M d^d\{0a} J o Jo Jo Jo + cw*j[ j[ W}(«3V^A{w8} + cos^(y-« y)^ jf W}{n3V^A{G8}] + ((AfJ - AfJ)-^ + (A^-Af?)-^ + (Afj - Af*)-^ [jf* jf' A{F1'}(F1') ^ ^ { ^ J + fin ft Z"1 / " ' A W X ^ V ^ A i v J - s i n ^ - c J C C A{Fi'}(n3') d£ dA{08} Jo Jo Jo Jo + a»tof f AWXnaV^AiwJ + cosV.Cy-^) / / A{F1,}(n3')d£dA{08} Jo Jo Jo Jo + ! > { £ j f W W ) dtdX + il j f / V s , « }<F3>«> dtd\ ip * +k I! f!{F»« }{F3")d*dx+k I! I! w><** > w J o J o J o J o -rCosV.^ jf W}{n3')d^A{w8> + cos0,(y-e i f)^ jf {F3'}{n3') ^ dA{Os}] + ((AfJ - + (AfJ-AfJ)jU (AfJ - Atf )£) jf 1 A{F3'}<F3') d*<*A{*J 46 + » i n * f l flX{F3'}MdidX{Y9}-8m^{z-eg) [* /'X{F3'}{nz')dfdX{Sa} Jo Jo Jo Jo + o o • ^ Y V , M ^ ' } { ^ ^ , ) ^ r f A { w . } + c o • ^ ( I f - « r ) / , / X{F3'}(nz')dfdX{&s} Jo Jo Jo Jo = (M° x M%/1 M*px A ^ / j f (5) Integration of the products of the shape functions gives the component stiff-ness matrices which will be fit into the total element stiffness matrix, according to their generalized coordinates. Care must be taken to adjust signs for the rotations and moments in the stiffness matrix since the nodal forces should be made to follow the usual finite ele-ment sign convention. The following sign changes will have to be made when forming the element stiffness matrix. w -u/(0) w'3 = -w'{l) Ai j = - M y ( 0 ) Mi = ~Mz(l) The five equilibrium equations can be rewritten in terms of the component stiffness ma-trices shown below. The right sides of the equations are the same force vectors as in the previous equations. The numbered component stiffness matrices can be found in appendix 2. ^%i]x{iio> = P (1) IT [* 3L { V s } + f M { V- } +T 1 M {e«> - My[~[K2] + [Jf4] + [K21]]{BS} - M>[\KA) + [K21)]{GB} + £ *f sin * { j [K5] {9J -Mi-[\Kb\ - [K7\] {©i} + M>j [K7] {•,} + Mtj- [[K6] - [K7]] {•,} - M>± [K7] {9,} } ip * v + M / ^ [KS] {©J + M i ^ [[JT,] - [Kt]] {»!> - M / | - [jr8] } ^ [*L { w->+ 7 W <w'>" W {e-> - Mi [-[K2] + [jy + [*21]] {e8> - Af/[[Jf4] + [Kn\] {OJ + ET C O S * {IM { e » } " N " ™ ] + MjlL [K7]m + MIJ-[[JST,] - [K7]\ - M>£ [ j f T ] } +ET c o s * { i W - * j f N - M + Mi j- [K9] {9,} + Mij- [[jr.] - [jr,]] - Mij- [jr.] { * j } ^[[^]-TO]{0.} + ^ [^4]{0 8> + ^ [K2] {V 8} - Aij [-[JT2] + [JT4] + [tf21]r] {v8} - Mi [[KA] + [K2lf] {•.} " ^rN(w 8 } - M*[-[tf2] + [JSTJ + |/f2i]r]{w8} - Mi[[KA] + [*21]r]{w8} -Mij{- sin - O + cos *(y - e , ) ) [[JT,] - [ tf7]] {*i> + sin - + cos *(y - [*7] {*i> v + Mtj-(- sin * ( « - O + cos fc(y - c y ) ) [[Kb] - [K7]] {*J - Af/|-(- sin *(* - «,) + cos &(y - [K7] {*i} } 48 + E {I (~ 8 i n ~° + 0 0 8 ^ '(y " c»}) M - M ^ ( - sin fi{M - ez) + cos*(y - [[JT,] - [K8\] {*{} + M ' f (- sin *(s - e.) + cos 4(y - [jT8] {•,} + A i : ^ ( - s i n ^ - e J + cos^(y-cff))[[A-6]-[if8]]{*i} . - Af/|-(- sin - «,) + cos * ( y - e,)) [jf,] } = M , (4) E ^ W w + UKl0\{*i} + w\[Kn] + ™ r]<**>+ V K ] { * i } } + ( - s i n i>i(z - « , ) + c o s fc(y - « , ) ) [ t f 5 ] T { © s } ) - JM{i - [if 14]] { * i } + sin ^ [[Kbf - [K7f] {vj + c o s * [[K5]T - [K7]T] {w8} + ( - s i n * ( , ~ O + c o s *(y - ev)) [[KS]T - [K7f] {ej) + M'j- ([tf 1 4]{ftj 4- s i n * [ K 7 ] T M + c o s * [ j r 7 ] T { w B } + ( - s i n - eJ + c o s ^ (y - «„)) [/^ {©j) + Ml\j- ([[Ku] - [K1A]] { • , ) + s i n * [[Kb]T - [K7]T] {v8} + c o s * [ [ t f 5 ] r - [K7f] { w 8 } + ( - s i n *(s - O + c o s f c ( y - e„)) [[Ksf - [K7)T] {0.}) - M i f ([ifu] { * i > + dn * [JTT] *W + c o s [ ff7] r{w8} + ( - s i n rf>i(z - O + c o s f t ( y - ey)) [ /^{Oj ) J + E D { % [KA <*i> + [ * « ] <*i> + W [ l * n ] + l ^ i 7 ] r ] ( * i > + [*is] } + ( W ^ + ^ ^ M V . } + COSV,[if6]r{w8} 49 + (- sin 1>i{z - eg) + cos tf,(y - [jr6]r{e8}) - Mij-(j[K19]-[*„]]<#.} + sinfc[[*.]r - [if8]r]{v8} + cos^[[if6]r - [#8]r]{w8} + (-sin *(* - O + cos *(y - «„))[[/f6]r - M T ] {©.}) + JMjJj- ([jTao] {*i} + sin * [K&] T { V 8 } + cos * : [ j r , ] T { w i > + (- sin - cj + cos MV ~ [^{©j) + J4J- ([[*»! - [*2o]] + * i * [[if«]T - [K8f] {v8} + cos * [[JT,]r - [A,]r] {w8} + (- sin -et) + cos*(y - e,)) [[JT6]r - [*8]r]{©,}) - M/|- ([iSTao] {*i} + sin * [ * » ] T { T . ) + cos * [jT,] r{w8} + (- sin - e.) + cos V,(y - e,)) [K8]T{&,}) } = M p (5) 3.12 ELEMENT MATRICES Before the stiffness matrices in equations (1) to (5) of the last section are actually imple-mented in a finite element program there are a few adjustments to be made. Up to this point in the analysis the plate degrees of freedom have been given in the order of the beam nodes, i.e. for a cross section with two internal plate nodes, m and n, <f>T = ((4? 4? V 43')m{4i 4? 4> p\). In the final element stiff-ness matrix it is necessary to group all of the plate degrees of freedom from each beam node together, so that the plate degrees of freedom will be ordered as follows: 4T — {{4m 4>n 4n)* {4m $n)* )• This means that some of the rows and columns will have to be interchanged in the matrices involving internal plates when these matrices are actually implemented in the program. Another detail to be considered is the 'compatibility* of units in the general-ized force, displacement and stiffness terms. Some of the stiffness terms must be multiplied through by the element length so that the force units come out as required. 50 Since the element cross section is arbitrary, the number of plate nodes on the section and therefore the overall size of the element stiffness matrix is variable. The number of plates contributing to each node will also be variable, as indicated by the summations in the equilibrium equations (see figure 3.4 for examples). Figure 3.9 shows a schematic of how the different parts of the element ma-trices are put together to form the element linear and geometric stiffness matrices in local coordinates. The linear, geometric, beam, plate, and coupled parts come from the ap-propriate terms in the equilibrium equations. All of the terms involving plate degrees of freedom are built up from the summations over the internal and external plates as appro-priate. The actual matrices themselves are given in the next few pages. Symbols used to simplify certain expressions are given before each matrix as they apply. -i r LB LB beam LP LP plate LB LB beam LP LP plate LINEAR STIFFNESS GB GBP GB GBP beam GBP GP GBP GP plate GB GBP GB GBP beam GBP GP GBP GP plate GEOMETRIC STD7FNESS FIGURE 3.9 Schematic Stiffness Matrix. -12 + t o I H 2 1 t o 1 3-"12 + t o 1 S > -t o 1 t o s i 5 -• M 1 7 -12 1 t o 1 + t o mm & 7 t o 1 -12 + t o + t o t o 1 who + n +» S S 2h II 0) 52 Linear Plate Stiffness External Plate Terms: A = D B = D b^ P lb b is the width of the external plate attached to plate node m. This matrix will be calculated for each external plate connected to this node. LPE = f 1 4A + \B 2A + &B I \ »A+A* 4^4 + J B sym. il* + A* / Internal Plate Terms: 1401s 2520 4s 90016 Here 6 represents the width of the plate between nodes m and n. K 'K 'K < IK #i IK ft. r 16C + 3744£ BC + 528£ - 1 6 C + 1296£ 8 C - 3 1 2 £ - 1 2 C + 1872£ -6C + 264£ 12C + 648£ -6C - 1S6£ +288F + 144C7 + 144f + 12G -288F - 144C +24.F + 12G -12F - 36G -36F - 3G +72F + 36G -6F - 3C J*C + 96£ - 8 C + 312£ \C~ 72£ -6C + 264E - 4 C + 48£ 6C + 1S6£ - 2 C - 36£ +32/W6G - 2 4 f - 12G -8F - 4G -36F - 3G - 8 F - 4 G +6F + 3G +2F + G *i I6C + 3744£ - 8 C - S 2 8 £ 12<? + 648£ 6(7+ 156 E -12(7 + 1872£ 6C - 2G4£ +288F+ 144G -U4F- 12G +72F + 36G +6F + 3G - 7 2 f - 36G +36F + 3G wL J fC + 96£ - 6 C - 156£ - 2 C - 36£ 6C - 264£ - 4 C + 48£ +32F + 16G -6F - 3G +2F + G +36F + 3G -8F - 4G *i 16C + 3744£ +288F + 144G 8C + S28£ +144F + 12G - 1 6 C + 1296£ -288F - 144G 8 G - 3 1 2 £ +24F+ 12G ^ C + 96£ +32F + 16G - 8 C + 3 I2£ -24*" - 12G \C - 12E -8F - 4G symmetric 16G + 3744£ +288f + 144C -8C - S28B - 1 4 4 F - 12G i J fC + 8«£ +32F+16G j Geometric Beam Stiffness u' P lu/"' l/ 1*-' e> M> /»" \ 1 if 36P 36Pe, ZP 3Pe, +3A<; -36P - 3 6 P e , - 3 A / J + 33 A / / 3 P 3P«, -3M} vf 36P -36Pet +33MJ - 3MJ -3 J» -3Pt, +3MJ - 3 6 P 36Pe, -3JMJ + 33Atf -3P -3Ptt -3M! 8' 18JTJ +18*/ 3Pe, 3Pe, -3Af* - 6Mj 3K'f -36Pe, -33Mi + 3MJ 36Pe f - 3 3 M J + 3 M i - 1 8 K J - 1 8 A ; 3Pet -6M\-3M', 3P«, +6Mi + 3M', 3KJ Iv/ 4 P 4Pe, -3M) + M! 3 P - 3 P e , -3AfJ - 6 M / -P -P'. +Mt i/ 4 P We, + 3 M J - . M J - 3 P - 3 P « , +3JWJ + 6M,J -P -P'. -Mi «*' 3K', + K'f - 3 P e , - 3 M * 3Pe , -3K'f -Pe, -MI -P', +M', vf v> 36.P 36Pe, +3Mi-33M', -3P - 3 P e , +3A<jJ w> symmetric 36P - 3 6 P e , +3Mi - 33Af/ 3 P 3Pe f + 3 « / V 18.K* + 18K} - P . , +6JWI + 3 M / - 3 P e , -6 JM; - 3Jtf/ -3Ki Ivf' 4 P APt, -Mi+3M! It/' 4P 4Pc, Mi-3M', 18" Ki + 3Ki j Geometric Plate Stiffness External Plats Termi • _ As A r t . M'y A / . ' - T ( - ^ - - — ) P. = kbP A 1 A b and h refer to the width and thickness, respectively of the external plate connected to node m. This matrix will be calculated for each of the external plates connected to node 36^ + 18A/J + ISA// 3^ + 3A// -36^ - 18A/J - ISA// 1 \ 3^ + 3A/J 45 + 3Aii + A// -3^ - 3A// -J^-JM/ 36^ + 18A/J + ISA// -3$-3A/J sym. <^ + Afl' + 3A// j Internal Plate Terms: P hbP • _ U J g i A//y / 1 / , ~ / , J where the quantities A and ft apply to the plate between nodes m and n. «& t*'i ( - 4 8 f 45 -365 "35 365 - 3 5 +24A/"' + M A / ' +4A/' -24AT - 24A/' +4AT - J8AT - ISA/' -3A/J +18AT + ISA/' - 3 A T ¥ 5 "45 -15 "45 35 5 +4A/*' + }A/' - 4 A / ' - | A T - }A/' - 3 A/' -3AT - M> +3A/' +JA/- + JA/' «£ 485 "45 365 tP 3Z -365 3* +24AT + 24A/' -4A/" + 18A/"' +18A/' +3A/' - 1 8 A T - ISA/' +3A/*' i«P TI p J 3* -«5 +JAT + 4A/' -3A/« +{AT + {A/' +3AT - A T - 3A/' 485 +24AT + 24A/' 45 +4A/' -485 -24AT - 24A/' 45 symmetric *5 +4AT + jA/» - « 5 -4A/» - J 5 -JAT-JA/' «£ 485 +24AT + 24A/' -45 -4A/* k ¥ 5 +}AT + 4 A / ' J Geometric Beam-Plate Stiffness External Plate Terms m « - - i r + ~ r m't = — y-„ hb . , C„ = y CO8 0 SCe = y(-sin^(z-e,) + cos^(y-« v)) 6, A and ip are the width, thickness and angle, respectively, of the external plate connected to node m. This matrix will have to be calculated for each external plate connected to node m. 4L 'K f u*' Se(36% + 18mj + 18m'.) Se(3$ + 3m>) Sc(-36% - 18m* - 18m1.) 5e(3$+3mj) to'' c.n c.n c.n <?.(") e( SC.?) s e n s e n SC.(») to*' C e(3§ + 3m{) (7 , (45 + 3m', + mi) e{-3%-3m{) to"' s.n s.n s.n SCC(") s e n s e n u> Se(-36% - 18m* - 18m£) Se{-3$-3m>) St(36% + 18m; + 18m£) 5«(-3$-3m') w> c.n c.n c.n B> SCt{') s e n s e n s e n to"' C e ( 3 5 + 3m') C « ( - 3 5 - 3 m ' ) + m* + 3mJ) Iv" ' Se{') s.n s.n s.n 18" { s e n s e n s e n s e n J Internal Beam-Plate Terms: m _ - _ + _ '» *• . M>z M'y "7 7~ S = y >in^ C = y C O » * As SC = —(- sin « ( » - « , ) + cos «l(y - «,)) where A, 6 and ^ apply to the internal plate between nodes m and n. *L «?*. K «.'• »•' S(36j + 18m'' + 18m>) S[3$ + 3m>) S(-36j - 18m'' - 18mJ) S(35 + 3m') S(-3e5 + 18m1' - 18m') S[-3$-3m>) 5(36$ - 18m'' + 18m') S(-35+3m') C(') C O C O C O co C O C O co «•' *?(*) S C O S C O S C O S C O S C O S C O sen C(3j J- 3m») C(4$ + 3m'' + m') C(-3$ - 3m>) C(-3J[ -3m') C(-45+3m'"-m') C ( 3 j + 3m') C ( 5 - J m ' + Jm»-) is so S(') so so so so so so 1ST sen S C O S C O sco S C O S C O S C O S C O «' v> S(-36j - 18m'' - 18m') S ( -3$-3m') S(365 + 18m'' + 18m') S(-35 -3m ; ) S(365 - 18m'' + 18m') S(3j + 3m') S ( -36 j + 18m' - 18m') S(3j - 3m') v> C O co co co C(') C(") C(-) C O *»' sen S C O S C O S C O SC(-) sco SC(*) S C O C(3$ + 3m') C ( - 5 - J m ' - J m » C( -3$ -3m' ) c(45+ m '+ 3 m ' ) C ( - 3 j + 3m') C ( $ - J m ' + Jm') C(35-3m') C(-4$ + m' -3m' ) so so so so so so S(-) S(') is" I S C O S C O S C O sco S C O SC(') S C O S C O J 57 3.13 TRANSFORMATION The element stiffness matrix has now been derived in local coordinates. This stiffness must eventually be transformed to some global coordinate system so that all the element stiffnesses can be put into the global structure stiffness matrix. In anticipation of the fact that the program NISA is going to be used to implement this element, the transformation used in reference [10] is used here. A rather complicated transformation procedure is used, based on the angle fi, which is defined as the angle between the local z axis and the intersection of the plane of the cross section and the plane defined by the local z and global z axes. The same transformation is used for <f> and <p' degrees of freedom as is used for the 0 and 9' degrees of freedom. 58 CHAPTER 4 Implementation The next step in the process is to use this new beam element in a finite element program. It was decided to use the nonlinear finite element program NISA, devel-oped at the University of Stuttgart, and use the thin-walled beam element developed by P. Osterrieder, for this program, as a starting point for this new element. NISA itself will be described briefly in section 4.4. First there will be a short decription of the beam element currently available with NISA, and then a general description of the changes made to accomodate this new "superbeam" element with its plate degrees of freedom. 4.1 NISA'S THIN-WALLED B E A M ELEMENT Osterrieder developed his thin-walled beam element from a regular beam element by adding a 9' degree of freedom at each end of the element (the derivative w.r.t. the longitudinal axis of the rotation about this axis). This allowed the effects of warping displacement (which is proportional to 9) to be taken into account. The thin-walled beam element then has 14 degrees of freedom per element (9' as the 13th and 14th d.o.f.). Cubic shape functions were used for the transverse displacements and for axial rotation, and a linear shape function for the axial displacement. An incremental solution using the Updated Lagrangian approach is used for the large deflection/large rotation problem. Material nonlinearities are dealt with by using 59 a bilinear stress-strain relationship with kinematic hardening. Reference [1] (in German) gives more detail on the development of this element, including some numerical results. 4.2 INCORPORATING SUPERBEAM The element subroutines of the previously existing thin-walled beam element were modified and augmented to form the new superbeam element. The objective was to add this new element to the NISA element library without making changes to the main program or storage. If the number of plate degrees of freedom is the same for all elements the only new input data required is Poisson's ratio. The number of plate degrees of freedom is then calculated as the total number of d.o.f. minus 7 (beam d.o.f.) for all nodes. It may be desirable eventually to be able to vary the number of plate d.o.f. for different nodes. The input routine in NISA had to be altered slightly to allow for more than 7 d.o.f. per node and in the main routine the superbeam element was simply added to the list of available elements. Parts of the element routines dealing with material nonlinearity were re-moved for the sake of clarity. The existing eleme routines were then used to calculate the beam stiffness terms. New routines were added to change the order of the beam degrees of freedom (specifically to put the 6' d.o.f. with the rest of their respective nodal d.o.f.), to calculate the linear and geometric plate stiffness terms and coupled beam-plate terms, and then to place the stiffness terms in the correct positions in the element stiffness matrix according to their d.o.f. Some problems arose in dealing with the dynamic storage when implementing this new element. The easiest solution was to create another one-dimensional array for temporary storage of element matrices and vectors, which now have variable dimensions, depending on the cross section being analysed. Initially some limitations have been set on the capabilities of the element subroutines. As already mentioned, the element will deal with geometric nonlinearities 60 only, and not with material nonlinearity. Only one material definition is allowed. The number of beam degrees of freedom is fixed at seven. The plate degrees of freedom must be even and equals the total d.o.f. - 7. The maximum number of plate degrees of freedom per node is 12 (6 internal plate nodes) and the maximum number of plates making up the cross section is 10. These last two limitations can quite easily be extended if more complicated cross sections are being analysed. 4.3 EIGENVALUE ANALYSIS The development of the element stiffness matrices using an incremental procedure was done with the intention of an incremental solution to the nonlinear problem of thin-walled beam stability. The analysis will however first be limited to a linear eigenvalue analysis, which gives only the critical load for a 'perfect' structure. This type of analysis assumes that the initial linear stiffness matrix does not change. The total stiffness is obtained by adding to this a multiple of the 'initial' geometric stiffness. With an instability problem there are two infinitesimally close configurations which are both in equilibrium with the given loading at the point of buckling, ie. there is more than one equilibrium configuration at the critical load level. The problem can be set up as follows (see ref. [13]): [Kt + \Kg\{D} = [K, + \Kg\{D + 6D} = R (4.1) where K\ is the initial linear stiffness, Kg is the geometric stiffness, A is the load factor (and eigenvalue), D and 6D are the total and incremental displacements, respectively, and R is the load vector. Configurations D and D + 6D are both in equilibrium with R. From this equation the following must be true, [KI + \K,]{6D}=0 This equation is then solved for A. Kg requires values for internal forces to be calculated and therefore one incremental load step is needed obtain the forces before the eigenvalue 61 analysis is performed. The eigenvalue, A, indicates the intensity of a given load distribution, that will cause instability, so that the critical buckling load is given as Ait*. The buckled shape can be calculated by putting the solution for A back into equation 4.1 and solving for the displacements. These displacements give the mode shape, but the actual magnitudes of the displacements are meaningless since the system is now unstable. The aim at this point is to solve the linear eigenvalue problem using NISA and see whether the coupled analysis will give buckling loads significantly different from those obtained when analysing local and overall buckling separately. 4.4 NISA NISA (Nonlinear Incremental Structural Analysis) is a finite element program designed for analysis of large nonlinear structural problems. The program has several elements available in its library, including a truss element, beam elements, and plate/shell elements. The program is modular in its construction, with a relatively short main routine. It uses an efficient scheme of dynamic storage, overwriting storage space where possible. The block storage system facilitates out-of-core solution of a system of equations, which is more efficient for large problems. Banded matrices are stored in a one-dimensional array using the skyline method, avoiding unnecessary operations on zeros. The equation solver is based on Gaussian elimination. The coefficient matrix is triangularized and then the solution vector obtained by back-substitution. The program has a restart capability so that an incremental analysis can be stopped part way through, the intermediate data stored, and then restarted again from that point in the analysis. The eigenvalue solver is also designed for efficiency in dealing with large problems. The eigenvalue solution process will be described in some detail in the next section. Total Lagrangian and Updated Lagrangian formulations are available for some elements, however the beam element uses only the Updated Lagrangian (U.L.) method. 62 For incremental analyses different iteration procedures are available, such as Newton-Raphson, modified Newton-Raphson, and variations of displacement controlled (constant arc length) methods. Dynamic analysis is partially implemented in NISA, but is not yet opera-tional for the beam element. Figure 4.1 gives and outline of the buckling solution in NISA. Flags set in the program or input by the user determine the flow as shown. Input Phase : read and store input data (boundary conditions loads, buckling solution parameters, etc.) Solution Phase: set storage parameters for solution INITD: initialize displacements set solution Bags STATS: static analysis LOADEF: effective loads ASSEM: assemble stiffness matrices needed for solution NISAEL: calls element subroutines linear and geometric stiffness COLSOL: solves system of equations for displacements print out displacements STRESS: calls element subroutines to calculate and print element forces ASSEM NISAEL: geometric stiffness from element routines NIBUCL/BUCSS: buckling analysis set storage for buckling solution SSPACE: subspace iteration to find smallest eigenvalues WRTBUC: print out eigenvalues and buckling modes FIGURE 4.1 Buckling Solution There are two other areas that merit some explanation; one is the element subroutine itself and the other is the subspace iteration method used by NISA to solve the eigenvalue problem. The set of element subroutines is relatively straightforward, and is shown in figure 4.2. 63 Input Phase: read element input data, initialize displacements ,'mternal forces to 0 Solution Phase: set solution storage parameters for Superbeam element TEAM: element stiffness UPDATE: update element orientation angle, nodal coordinates and transformation matrix for U.L. formulation FORCE DISLOC: element local displacements from global displacements GEOPLST: geometric plate stiffness terms LINPLST: linear plate stiffness terms beam stiffness terms, linear and geometric solve for element forces P = KA UNBALFOR: unbalanced forces (for incremental solution) dear element stiffness matrices LINBMST: linear beam stiffness terms SWITCHBM: adjust order of beam d.o.f. PLACBM: place beam terms in element stiffness matrix LINPLST SETUP: set up arrays needed for plate calculations loop over internal plates PLAPROPL: plate properties LINPLINT: linear plate stiffness PLACINTL: place plate terms in element stiffness matrix loop over external plates PLAPROPL: plate properties LINPLEXT: linear plate stiffness PLACEXTL: place plate terms in element stiffness matrix GEOBMST: geometric beam stiffness SWITCHBM: adjust order of d.o.f. PLACBM: place beam terms in element stiffness matrix GEOPLST SETUP loop over internal plates PLAPROPG: required plate properties GEOPLINT: geometric plate terms GEOBPINT: geometric coupled terms PLACINTG: place geometric plate K in geometric element K loop over external plates PLAPROPG: plate properties GEOPLEXT: geometric plate terms GEOBPEXT: geometric coupled terms PLACEXTG: place in element K add geometric stiffness to linear stiffness to get element stiffness matrix in local coordinates TRANSF: transform to global coordinates ADDBAN: add element (global) stiffness to structure stiffness FIGURE 4.2 Element Stiffness Routines 64 note: It is also possible to calculate only the linear or only geometric stiffness when re-quired. The stress calculation in the element routine involves simply calling the routines UPDATE and FORCE to calculate element forces. The subspace iteration procedure is quite complicated and some background knowledge is required. It is required to solve the equation K,x = -XKgx (4.2) for the smallest eigenvalue, A, and the displacement vector, x. The subspace iteration method is just one of many algorithms available for solving the standard eigenvalue prob-lem. It was developed by Bathe (see ref. [14]) and aimed at the efficient solution of aout-of-core systems and systems with large bandwidths". It is also efficient for the simul-taneous calculation of more than one eigenvalue. The subspace iteration method is based on the inverse iteration method of solving an eigenvalue problem. This method solves the problem iteratively by first assuming a A and z on the right side of equation 4.2, solving for z on the left side, then using this new z in the right side again, and iterating in this way. The computer implementation of this method however, becomes a little more involved. The program chooses q starting iteration vectors, where q is greater than the number of eigenvalues to be calculated, and uses "simultaneous inverse iteration on these q vectors and a Ritz analysis to extract the best eigenvalue and eigenvector approximations from these q iteration vectors... the iteration is equivalent to iterating with a ^ -dimensional subspace" (a subspace being defined by q linearly independent (basis) vectors). The starting iteration vectors are chosen to obtain convergence as fast as possible. After the inverse iteration is performed, the Ray lei gh minimum principle is invoked to find the best eigenvalue and eigenvector lying in the subspace spanned by the Ritz basis vectors. A Gaussian elimination type of solution is used to solve for the eigenvector, z, at each iteration. Convergence is reached v/hen the 65 change in the eigenvalues in successive iterations is less than a specified tolerance. After convergence of these 'best* eigenvalues, the eigenvalue is checked using the Sturm sequence property. 66 CHAPTER 5 Results and Conclusions 5.1 DISCUSSION OF RESULTS After checking the element stiffness matrices, some simple analyses were performed for further checking of the new element. Examples with plate degrees of freedom fixed were checked against results from the previously existing thin-walled beam element in [1]. Adding plate d.o.f. should always decrease the critical load since it is making the element more flexible. When the beam d.o.f. are fixed, plate buckling should be observed independently. Symmetry should be shown in mode shapes and element forces where ex-pected. Convergence to the theoretical solution should be observed as the finite element mesh is refined. After these checks were made, some examples of different cross sections were analysed for the interaction of local and overally buckling. 5.1.1 Theoretical Calculations This section will review some of the theoretical calculations used for com-parison with the finite element results. From equilibrium considerations of a beam column member in the y, z and 6 directions, three differential equations are obtained which can be used to find the buckling loads. If the shear centre and centroid coincide, these three equations are uncoupled and each gives a critical load for buckling in the appropriate mode, ie. Euler buckling about the 2, and y axes and torsional buckling. The expressions for these three uncoupled critical 67 loads are: P x2EIz 1*2 P0 = (JG + **EIU. A I2 }IV + JZ If the shear centre and the centroid are not at the same point, the three equations are coupled and a resulting coupled mode may give a lower critical buckling load than any of the uncoupled modes. centre from the centroid, along the major axis, so that major axis buckling is coupled with torsional buckling. When solving the resulting quadratic equation, this coupled strong axis-torsional mode gives buckling loads slightly lower than pure torsional buckling, for short members (torsional buckling is the governing uncoupled mode for this case). As the member length increases, the governing buckling mode changes from this coupled torsional mode to weak axis Euler buckling. For a simply supported beam with equal end moments, the critical moment to cause lateral buckling is given by (ref. [18]): where / is the moment of inertia about the weak axis. From the fourth order differential equation for plate buckling under constant in-plane loading (in one direction), an expression for the critical plate buckling stress is obtained. D is the plate bending stiffness, 6 is the plate width, and is a plate buckling factor which depends on the plate boundary conditions, dimensions and (for most cases) the buckled shape. Although k actually varies slightly, it can be taken to be constant at a minimum For angle and channel cross sections there is an eccentricity of the shear When a beam is bent about its strong axis it may fail by lateral buckling. 68 FIGURE 5.1 CONVERGENCE OF PLATE BUCKLING LOAD 1.7-1 NUMBER OF ELEMENTS value for plates with length-to-width ratios greater than one. Assuming that the plates are simply supported along their connecting longitudinal edges, k can be taken as 0.425 for external plates (one longitudinal edge simply supported and the other free) and k=4 for internal plates (both longitudinal edges simply supported) (see refs. [15], [16], etc.). For both of these cases the ends of the plates are assumed to be simply supported. 5.1.2 Convergence Convergence for plate buckling of this element is shown in figure 5.1. An angle section is used since it provides a simple case of two external plates connected along their common edge. The figure shows monotonic convergence, from above, to the theoretical solution as the number of elements is increased. It can be seen that convergence to the theoretical solution is quite rapid and the finite element solution is very good even with only two elements. 69 For cross sections with internal plates the plate buckling critical loads fluctu-ate slightly for long member lengths and more iterations are required to obtain convergence of the eigenvalues. For a plate simply supported along one edge and free along the other, the buckled shape is one half sine wave, but for a plate simply supported on both edges, the buckled shape is a number of waves, depending on the plate dimesions. This shape cannot be accurately represented with just a few elements. Convergence of the coupled beam-plate solution was also difficult for long members for this reason. Convergence of the finite element solution as the mesh is refined was not always obtained for long members with internal plates. 5.1.3 A n g l e Section Figure 5.2 shows the first three mode shapes for plate buckling of an angle section. Boundary conditions were set to force symmetry about the centreline. This plot 70 FIGURE 5.3 ANGLE COLUMN BUCKLING O 4 ( /) ( / l I— to < <J Tors ion aJ-f lexur&l H3-Q—'^# © — A a Legend THEORETICAL -PLATE THEORETICAL —TORSIONAL COUPLED THEORETICAL -EULER O F.E. -LOCAL • F.E. -OVERALL A F.E. -COUPLED "A A - A 10 20 30 L/B 40 50 60 represents the actual variation of the plate node rotation, <f>, or (for the case of external plates), the perpendicular displacement of the free edge of the plates. The plate buckling loads are very close to those calculated theoretically, indicating that the element gives a good representation of external plate buckling. A cruciform cross section was also tested and gave correspondingly good results (equivalent to two angle sections). Figure 5.3 shows the different types of buckling for the angle cross section. Since the shear centre and centroid do not coincide there are two distinct theoretical overall buckling modes; a coupled torsional-flexural mode for short columns and pure Euler for long columns. The finite element solution also shows this change in governing buckling mode. The finite element solution for the torsional coupled mode is significantly lower than the theoretical solution, however, and the reason for this is unknown. On the other hand, 71 FIGURE 5.4 WIDE FLANGE COLUMN BUCKLING 50-i 0 5 tO 15 20 25 30 35 40 45 L/Bw the finite element torsional buckling stresses are very close to the plate buckling stresses, which would be expected, since for this cross section, torsional buckling about the shear centre is equivalent to buckling of the two plates about their common edge. Because of this fact, when the coupled local-overall buckling mode is analysed, the 9 and 4> degrees of freedom are redundant in that they both represent rotations of the cross section about the shear centre, 9 with stiffness due to St Venant and warping torsion, and <f> with stiffness due to plate bending. Because of this overlap the buckling stresses for the coupled analysis are lower than they should be. 5.1.4 W i d e Flange Section For a doubly symmetric wide flange cross section, torsional buckling is not a governing buckling mode. Figure 5.4 shows a plot of critical stress vs. dimensionless length (member length/ web depth: l/Ba) for an I-section with axial loading. Theoretical 72 FIGURE 5.5 WIDE FLANGE COLUMN BUCKLING 60 50 1 ft. _ l 30 Legend D THEORETICAL -PLATE (flange) - 1 THEORETICAL -EULER 1 O F.E. -PLATE 1 • F.E. -EULER 1 A F.E. -COUPLED 1 t , r r 3 1 to — 9 1 X 1 - L*-Lo e — l 4 11 0 ' © \ \ - I I r ! 4J 1 $ 0 5 10 15 20 25 30 35 40 45 L/Bw buckling stresses are shown for column buckling and flange plate buckling, along with the finite element results of column, plate and coupled analyses. It can be seen that the finite element results follow the theoretical results exactly for the overall buckling case. For plate buckling, although a little on the high side, the computer model gives the same general behaviour as the theoretical curve. The theoretical plate buckling stresses were calculated for the flange being simply supported along one side. In an I-section the flanges are actually restrained somewhat by the web, and cannot rotate freely. Therefore higher plate buckling stresses, as shown by the computed results, would be expected. The computed results are lower than the theoretical values calculated using the plate buckling coefficient for a fixed-free plate, (k = 1.277, ref. [15]) however, indicating that the flange/web junction is modelled somewhere between a fixed support and a simple support for the flange plate. The coupled finite element analysis follows the plate buckling curve for short 73 FIGURE 5.6 WIDE FLANGE LATERAL BUCKLING CONSTANT MOMENT 450-E E 400-350-300-Legend THEORETICAL -PLATE (top flange) THEORETICAL -LATERAL O F.E. -PLATE • F.E. -LATERAL A F.E. -COUPLED 3 O 250-< O 200-150-100-50-—r-15 20 25 L/Bw 30 columns and drops slightly as it makes the transition to Euler buckling for longer column lengths. This indicates that the interaction between local and overall buckling is negligible, except for the range of column lengths where the buckling mode is in transition between the local and overall modes. Figure 5.5 shows the same type of plot except now for a section with smaller flanges. The theoretical plate buckling stresses for the narrower flanges are higher, but the finite element results still correspond fairly well, being slightly high because of the restraint of the web. The Euler load is lower because of the reduction in weak axis moment of inertia, but the interaction pattern in general is the same as in the previous figure. For the longer columns the plate buckling mode shapes show more than one wave in the longitudinal direction. Theoretically the edge plates (flanges) should show only one half wave in their lowest buckled mode shape for any column length. In this case 74 FIGURE 5.7 WIDE FLANGE LOCAL BUCKLING 25000-o a. 5 in 20000-15000-cc t— t/1 _ J < u t-5 IOOOO-o 5000-Legend o f.e. l o c a l theoretical - f lange theoretical - web flange 0.4 0.5 0.6 Bf/Bw 0.7 0.8 0.9 1.1 however, there is interference of the web buckling modes (in which the number of half sine waves in the buckled shape depends on the plate length) which forces the flanges into a similar pattern. The same holds for other cross sections in which the internal plates will affect the buckling of the external plates. The results of the analysis of an I-section with equal end moments is shown in figure 5.6. The cross section is the same as in figure 5.5. The same type of results are plotted as previously, and it can be seen that there is still no significant interaction. Ref. [4] claims that web buckling is more significant when bending is involved than for pure axial loading. This would suggest perhaps, that there should be more interaction of buckling modes for the case of applied moment. A plot of critical buckling stress vs. flange width/web depth ratio is plotted in figure 5.7. The web depth actually remains constant as the flange width is varied, for 75 FIGURE 5.8 CHANNEL COLUMN BUCKLING L/Bw a column length of 400 mm. The figure shows that the model represents flange buckling rather than web buckling, for the local buckling of and I-section. For sections with very-small flanges the model shows very high stresses required to buckle the flanges, rather than the lower stresses which would be required to buckle the web. For these sections, the buckling stress should actually drop down, as the local buckling mode changes to represent web buckling. The results imply that web buckling is not well represented by this element. In ref. [2], Rajasekaran makes a similar observation and suggests that using a cubic function for web deformation (and therefore only a linear variation of curvature across the web), will not represent the fundamental web buckling shape. 5.1.6 C h a n n e l Section Figure 5.8 shows yet another plot of critical stress vs. dimensionless length, this time for a simple channel cross section with axial loading. Here coupled torsional-76 Section Thickness E ayidd G'uh (exper.) crCT (P.E.) (in) (ksi) (ksi) (ksi) (ksi) 0.1046 17110 62.6 62.6 170.0 0.1345 19090 71.7 75.4 300.8 0.1046 23040 65.1 65.7 143.8 n 0.1345 19270 69.3 71.1 199.7 T A B L E 5.1 Hat Section Buckling strong axis buckling is the governing overall buckling mode for short columns. Again there is the problem that the finite element analysis underestimates the critical load when torsional buckling is involved. In the coupled analysis there is a significant reduction in critical stress, due to interaction of local and overall modes in the region where coupled torsional-flexural buckling is the governing overall buckling mode. Reference [5] quotes Bijlaard and Fischer when considering the possibility of interaction between local and overall buckling; "The interaction effect is negligible for box sections... the same conclusion applies to common sizes of H and channel sections, but not to sections for which torsional instability is an important factor, such as T and angle shapes". The channel section used in figure 5.8 has relatively wide flanges so that torsional instability is more of a factor than for narrower channel shapes. This may account for the fact that there is more significant interaction. Reference [6] shows some results for channel sections of different dimensions with both stiffened and unstiffened flanges. Their coupled analysis give results between the local and overall results, not below the local buckling loads as obtained in this study. 5.1.0 H a t Section Several hat sections were analysed for plate buckling only, using the super-beam element. The numerical results were compared to experimental results obtained from 77 a testing series given in reference [17], and are given in table 5.1. It can be seen that the finite element results do not correspond, in magnitude, to the experimental results, but the trend is the same, ie. higher buckling stress for the thicker sections and for sections with stiffened flanges. The main reason for the large discrepancy is that the hat sections tested in the lab all buckled after the yield stress was reached. This reinforces the statement that elastic buckling is often not sufficient to model real situations, in which some plasticity, especially for local buckling, is usually involved. The sections tested had very rounded corners, rather than well defined angles, and the dimensions of the sections were quite different from those specified by the manufacturer. Moduli of elasticity measured in the lab and used for the analysis were also different from those specified (partially due to cold forming). Some of the plate angles (^ ,-, see fig. 3.8) in the sections tested were relatively large (eg. 125° for the lips in the stiffened sections). Because of this, and the fact that they were quite rounded, cross sectional distortion, in the form of changes in these angles, was observed in some cases. The finite element model of plate nodes rotating rigidly is not a good representation for the large, poorly-defined corner angles of these sections. 5.2 CONCLUSIONS More analyses should be performed before judging the value of this element as a model for the interaction between local and overall buckling for thin-walled cross sections. From the results shown in section 5.1 however, some initial observations and tentative conclusions may be made. Some of these are reinforced by the work of other researchers. In general only a few elements are needed to give accurate results. This was also observed by Rajasekaran [2] and others. With long members however, there are problems in representing the local buckling shapes with only a few elements. Since only an eigenvalue analysis was done, any post-local buckling strength was not accounted for. Initial buckling, whether local, overall, or combined, meant collapse. 78 Actual interaction of the local and overall modes would be better observed by doing a load-deflection analysis and observing the deflections throughout the loading process. From the eigenvalue analysis it appears that the interaction of local and overall buckling is small for cross sections which do involve torsional buckling in their governing overall buckling modes. For these sections there is a small reduction in critical load due to interaction, for members with lengths in the range in which the governing buckling mode is in transition between plate buckling and overall buckling. Reference [2] also observed that local and overall buckling are "relatively detatched phenomena* for wide flange sections in the elastic range. For sections involving torsional buckling there is more significant interaction of these buckling modes, as observed in ref. [15]. The beam element does not model torsional buckling to the same degree of accuracy as it does other types of buckling. There seems to have been some problems in this regard, with the element as it existed previously in NISA. The buckling of external plates is modelled very well. When flanges or stiffen-ers of sections with well defined corners buckle, this new element will give a good prediction of the buckling behaviour. The buckling of webs and internal plates is not represented quite as accurately, perhaps due to the assumed shape function for the internal plates. This new element is easily incorporated in any analysis that used the pre-viously existing thin-walled beam element, since it followed the same procedure in its formulation. It would be relatively simple to add the effects of local buckling to any analysis where these beam elements are used. 5.3 F U R T H E R S T U D Y There are many aspects of this study which could be expanded upon to obtain more practical results. Firstly, it would be useful to use the incremental load-deflection approach to analyse the problem and observe the development of any interaction of buckling modes 79 for these simple members. This would them take into account post-local buckling strength and also the effects of initial imperfections and curvatures. An important extension to the analysis would be to include inelastic material behaviour, since buckling most often does not take place until some plasticity has developed in the section. This would allow more practical situations to be analysed and comparison with experimental results would be more meaningful, in the ranges where interaction is important. Residual stresses should also be taken into account, since they are significant in the behaviour of thin-walled sections. More results from other studies, both experimental and numerical, should be found to compare directly to analyses done in this study. This new element could be incorporated in the analysis of frames, trusses, open web steel joists, etc. and the results evaluated. Joining elements with different cross sections, different numbers of plate degrees of freedom, at different angles, etc. should be investigated. Also warping compatibility and the effects of warping restraint should be looked at in more detail. Although it appears that there is much more work to be done, the element as it is developed thus far gives reasonable results, is easy to use and is a good beginning for the analysis of the interaction of local and overall buckling of thin-walled cross sections. 80 References [l] Osterrieder, P., Traglastberechnung von Raeumlichen Stabtragwerken bei Grossen Verformungen mit Finiten Elementen, PhD Thesis (in German), Institut fuer Baustatik der Universitaet Stuttgart, July 1983. [2] Rajasekaran, S., Murray, D.W., Coupled Local Buckling in Wide Flange Beam Columns, A.S.C.E. Journal of the Structural Division, vol. 99 ST6, June 1973., pp. 1003-1023 [3] Hancock, G.J., Local, Distortional and Lateral Buckling of I Beams, A.S.C.E. Journal of the Structural Division, vol. 104, No. STll, Nov. 1978., pp. 1787-1798 [ 4] Hancock, G.J., Bradford, M.A., Trahair, N.S., Web Distortion and Flexu-ral Torsional Buckling, A.S.C.E. Journal of the Structural Division, vol. 106, No. ST7, July 1980., pp. 1557-1571 [5] Akay, H.U., Johnson, CP., Will, K.M., Lateral and Local Buckling of Beams and Frames, A.S.C.E. Journal of the Structural Division, vol. 103, No. ST9, Sept. 1977., pp. 1821-1832 [6] Wang, S.T., Pao, H.Y., Torsional-Flexural Buckling of Locally Buckled Columns, Computers and Structures, vol. 11, 1980., pp. 127-136 [7] Fukumoto, Y., Kubo, M., Buckling in Steel U-Shaped Beams, A.S.C.E. Jour-nal of the Structural Division, vol. 108, No. ST5, May 1982., pp. 1174-1190 (8] Ghobarah, A.A., Tso, W.K., Overall and Local Buckling of Channel Columns, A.S.C.E. Journal of the Engineering Mechanics Division, vol. 95, No. EM2, April 1969., pp. 447-462 (9] Wittrick, W.H., Williams, F.W., Initial Buckling of Channels in Compres-sion, A.S.C.E. Journal of the Engineering Mechanics Divison, vol. 97, No. EM3, June 1971., pp. 711-726 [10] Reis, A.J., Branco, F.A., Lateral and Local Stability of Thin-Walled Sections under Eccentric Loading, Thin- Walled Structures, Conference Proceedings, Glasgow, April 1979., pp. 705-717 [11] Graves-Smith, T.R., Sidharan, S., Elastic Collapse of Thin-Walled Columns, Thin-Walled Structures, Conference Proceedings, Glasgow, April 1979., pp. 718-730 References 81 [12] Fung, Y.C., Foundations of Solid Mechanics, Prentice Hall, 1965. [13] Cook, R., Concepts and Applications of Finite Element Analysis, John Wiley and Sons, 2nd edition, 1981. [14] Bathe, K-J, Finite Element Procedures in Engineering Analysis, Prentice-Hall Inc., 1982. [15] Johnston, B. (ed.), Guide to Stability Design Criteria for Metal Structures John Wiley and Sons, 3rd edition, 1976 [16] Timoshenko, S.P., Gere, J.M., Theory of Elastic Stability, 2nd edition, McGraw-Hill Kogakusha Ltd. 1961. [17] Toneff, J.D., Local Buckling in Cold Formed Members, C.E. 511 course project, U.B.C. 1983. [18] Narayanan, R. (ed.) Beams and Beam-Columns - Stability and Strength Applied Science Publishers, 1983. [19] Chen, W.F., Astuta, T., Theory of Beam Columns, vol. I: In Plane Be-haviour and Design, vol. II: Strength and Stability, McGraw Hill, 1976. [20] Bleich, F., Buckling Strength of Metal Structures, McGraw Hill, 1952. [21] Murray, N.W., Introduction to the Theory of Thin-Walled Structures, Clarendon Press, Oxford, 1984. 82 APPENDIX A l Symbols and Notation Symbols The following definitions apply to symbols used in the text: E Young's modulus ^ = 2(1+1/) shear modulus v Poisson's ratio A total beam cross-sectional area Iy moment of inertia about y axis Iz moment of inertia about z axis h plate thickness D = plate bending stiffness ui = — f rds sectorial coordinate ua — j fAudA — u normalized sectorial coordinate r perpendicular distance from shear centre to tangent on section 3 tangential coordinate along the cross section 7W = J^u^hds warping constant J = \/ZJ$hzds St. Venant torsional constant for an open cross section <r" normal stress due to warping a distance from shear centre to a point on the cross section Kp = JA aPo™ dA Wagner coefficient My moment about y axis Ms moment about z axis Vy shear in y direction Vg shear in 2 direction Mx torsional moment about x axis Afw bimoment S shear centre vo,vSiwS>9s rigid body beam degrees of freedom <f> plate node rotation r, a, t plate displacements in plate coordinate system Notation Some of the general notation used in this thesis is as follows: - { ) indiates a row vector - { } indicates a column vector -» as a right subscript represents a plate number in the cross section - J2 indicates a sum over all the internal plates in the cross section - J2 indicates a sum over all the external plates in the cross section ep (see section 3.3 for definition of internal and external plates) - i,j as right superscripts indicate beam node i or j 84 APPENDIX A2 Component Matrices / o 1 W } ( n 1 V A = ( _ 1 1 - 1 ) J\n3'}(nz')d\ = ± 36 3 -36 3 \ 4 -3 -1 30 | 36 -3 4 J -15 •15 3 -3 -1/2 3 0 [\n3"}(n3")d\ = Jo (12 6 -12 6 \ 4 - 6 2 12 -6 4 / / o 1 A { n 3 ' } { n 3 V A = ^ /18 3 -18 0 1 -3 -1/2 18 0 3 [ f*M(Fx')did\=l-[K2] j0lf0lM{F3')d^x = ^o f 36 3 -36 3 -36 -3 36 -3 3 4 -3 -1 -3 -4 3 1 -36 -3 36 -3 36 3 -36 3 V 3 -1 -3 4 -3 1 3 -4. f! I!xM{Fi')d(dx==\[K'} / o 7 . ' A { n 3 ' > ( W A = 3 o O / 18 3 -18 ^ 0 3 1 -3 -1/2 -18 -3 18 0 0 •1/2 0 3 -18 -3 -3 -1 18 3 0 1/2 18 0 3 1/2 -18 0 0 -3 f1 /1{F1,«}( JP1 > a)rferfA = (0) Jo fl /Vi>«}(^ i'V£<*A = [0] Jo jfJj[Vi*>to*}^«[*»] f jQ1X{Fl'}{Fl')d^dX = l[KA] Jo Jo /16 8 -16 8 -12 -6 12 -6 \ 16/3 -8 8/3 -6 -4 6 -2 16 -8 12 6 -12 6 1 16/3 -6 -2 6 -4 140 16 8 16/3 -16 -8 16 8 8/3 -8 16/3/ fl fl{F^}(F3^)d(dX Jo Jo /3744 528 1296 96 312 3744 2520 -312 1872 264 648 -72 264 48 156 -528 648 156 1872 96 -156 -36 -264 3744 528 1296 96 312 3744 -156\ -36 -264 48 -312 -72 -528 96 ) Jo Jo ( 144 132 -144 12 -36 -33 36 -3 ^ 12 16 -12 -4 -3 -4 3 1 -144 -12 144 -132 36 3 -36 33 12 -4 -12 16 -3 1 3 -4 -36 -33 36 -3 144 132 -144 12 -3 -4 3 1 12 16 -12 -4 36 3 -36 33 -144 -12 144 -132 I -3 1 3 -4 12 -4 -12 16 ) Jo Jo /144 12 16 -144 -12 144 12 -36 -3 36 -3 \ -4 -3 -4 3 1 -12 36 3 -36 3 16 -3 1 3 -4 144 12 -144 12 16 -12 144 -4 -12 16 J f1 ['wm^dtdx Jo Jo /48 4 -48 4 -36 -3 36 -3 \ 16/3 -4 -4/3 -3 -4 3 1 48 -4 36 3 -36 3 16/3 -3 1 3 -4 48 4 16/3 -48 -4 48 4 -4/3 -4 16/3/ f1 [lx{F3'm')dtdx Jo Jo /24 1 4200 4 -24 0 -18 -3 18 0 \ 4/3 -4 -2/3 -3 -1 3 1/2 24 0 18 3 -18 0 4 0 1/2 0 -3 24 4 4/3 -24 -4 24 0 -2/3 0 4 /
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Local and overall buckling of thin-walled beams and...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Local and overall buckling of thin-walled beams and columns using finite elements Toneff, Janine Diana 1986
pdf
Page Metadata
Item Metadata
Title | Local and overall buckling of thin-walled beams and columns using finite elements |
Creator |
Toneff, Janine Diana |
Publisher | University of British Columbia |
Date Issued | 1986 |
Description | Steel members with open thin-walled cross sections are used extensively in civil engineering structures. In addition to overall instability, (Euler buckling, torsional buckling, lateral-torsional buckling, etc.), the thin plates making up the cross section may themselves be susceptible to local plate buckling. The possible interaction of these two modes of buckling and its effect on member stability and strength is therefore of interest in the analysis of such members. The purpose of this thesis is to develop a tool for this type of analysis by adapting a finite element for a thin-walled beam-column of arbitrary cross section to include 'local' degrees of freedom to allow for cross section distortion. The formulation will involve geometric nonlinearities due to large displacements and rotations, but material behaviour will be limited to the linear elastic case. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-07-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0062925 |
URI | http://hdl.handle.net/2429/26334 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1986_A7 T66.pdf [ 4.71MB ]
- Metadata
- JSON: 831-1.0062925.json
- JSON-LD: 831-1.0062925-ld.json
- RDF/XML (Pretty): 831-1.0062925-rdf.xml
- RDF/JSON: 831-1.0062925-rdf.json
- Turtle: 831-1.0062925-turtle.txt
- N-Triples: 831-1.0062925-rdf-ntriples.txt
- Original Record: 831-1.0062925-source.json
- Full Text
- 831-1.0062925-fulltext.txt
- Citation
- 831-1.0062925.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062925/manifest