C f NUMERICAL SIMULATION OF THE NON-LINEAR TRANSIENT RESPONSE OF SLENDER BEAMS By BRYAN RUSSELL FOLZ B.Sc, Simon Fraser University, 1982 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF MASTER OF APPLIED SCIENCE in T H E 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 © Bryan Russell Folz, 1986 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y o f B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be g r a n t e d by the head o f my department or by h i s o r her r e p r e s e n t a t i v e s . I t i s understood t h a t c o p y i n g or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be allowed without my w r i t t e n p e r m i s s i o n . Department of /?/ y / £ £/-</£~/^/<£'dr^//^^ The U n i v e r s i t y of B r i t i s h Columbia 2075 Wesbrook P l a c e Vancouver, Canada V6T 1W5 Date A/K/t- /33j& 1Q\ ii Abstract A simple numerical solution strategy for predicting the transient response of slender ductile beams, exhibiting geometric and/or material non-linear behaviour, is presented in this study. In the theoretical development of the problem the governing non-linear equa-tion of motion, in variational form, for the bending and stretching of a Bernoulli-Euler beam is established. The numerical solution procedure is then initiated by employing the assumed displacement version of the Finite Element Method with 1-dimensional 6 -DOF beam elements. Elastic-plastic strain-hardening of the beam material is conveniently ac-counted for by means of the "mechanical sublayer model". Visco-plastic material behaviour is included in the analysis through a simple strain-rate dependent constitutive relation-ship. The equations of motion for the spatially discretized beam are integrated time-wise by means of the central difference method. A variety of examples are then solved and the results compared with solutions from other sources. In general, the numerical solution strategy yields an efficient and accurate modelling of the non-linear transient response of slender ductile beams. iii Contents Abstract ii Figures iv Tables v Acknowledgements vi 1 Introduction 1 1.1 Background 1 1.2 Purpose and Scope 4 2 Formulation and Analysis of the Problem 6 2.1 Introduction 6 2.2 Kinematics of a Bernoulli-Euler Beam . . . . . . . . . . 6 2.3 Variational Equation of Motion 10 2.4 Finite Element Formulation of the Equations of Motion . . . . 13 2.5 Modelling the Material Behaviour of the Beam 20 2.6 Solution Procedure Over the Time Domain 29 2.7 Summary 35 3 Numerical Examples 37 3.1 Introduction 37 3.2 Example Exhibiting Material Non-linear Behaviour 37 3.3 Example Exhibiting Geometric Non-linear Behaviour . . . . 42 3.4 Examples Exhibiting Geometric and Material Non-linearities . . . 44 4 Conclusions 56 4.1 Summary 56 4.2 Future Research Areas 58 References 59 Appendix A : Stability Criterion for the Central Difference Method 63 iv Figures 2.1 Bernoulli-Euler beam kinematics 8 2.2 Typical beam finite element 14 2.3 Shape functions 16 2.4 Uniaxial stress-strain behaviour of a ductile material under quasi-static cyclic loading conditions . 22 2.5 The mechanical sublayer model 23 2.6 Effect of strain-rate on the yield stress 25 2.7 Idealization of strain-rate behaviour for an elastic-plastic strain-hardening material . 2 5 3.1 Simply supported beam subjected to a step load 39 3.2 Dynamic small displacement elastic and elastic perfectly-plastic responses of a simply supported beam 39 3.3 Strain history of a simply supported beam . . . . . . . 41 3.4 Comparison of numerically predicted response histories of an elastic-plastic simply supported beam 41 3.5 Clamped beam subjected to a point step load at midspan . . . . 4 3 3.6 Large deflection elastic dynamic response of a clamped beam . . . 43 3.7 Impulsively loaded clamped beam 46 3.8 Finite element modelling of the impulsive load 46 3.9 Comparison of predicted and experimental responses of an explosively loaded clamped beam 48 3.10 Strain history of an explosively loaded clamped beam 48 3.11 Deflection profiles of an explosively loaded clamped beam . . . . 49 3.12 Axially constrained pin-ended I-beam subjected to a uniformily distributed step load of finite duration ' 51 3.13 Dynamic large displacement elastic-plastic response of an axially constrained pin-ended I-beam 51 3.14 Deflection profiles of an axially constrained pin-ended I-beam . . . 53 3.15 Strain history of an axially constrained pin-ended I-beam . . . . 53 3.16 Interaction between axial force and moment over time for an axially constrained pin-ended I-beam 54 V Tables 2.1 One dimensional Gauss-quadrature 27 2.2 Flowchart for the solution algorithm 36 vi Acknowledgements The author wishes to thank Dr. M.D. Olson and Dr. D.L. Anderson for their advice and guidance throughout the execution of this work. As well, the financial support of the Canadian Department of National Defence, through a contract from the Defence Research Establishment Suffield, is gratefully acknowledged. Finally, the author is deeply appreciative of the support and encouragement offered by his wife, Ruby, during the preparation of this thesis. CHAPTER 1 INTRODUCTION 1.1 BACKGROUND The response of structures exhibiting large deformations and inelastic ma-terial behaviour under transient dynamic loading conditions is a subject which warrants considerable attention. The susceptibilty of structures to exceptionally high transient pres-sure pulses which are capable of causing failure is an ever present concern to the structural analyst. Potentially catastrophic loading conditions can be caused by explosions, impact, or earthquakes, thus affecting many areas of structural design. For example, the possib-lity of explosions occuring in storage containers and piping networks in the chemical and nuclear industries presents a very real danger. In such a situation the facility in question must be "fail safe", so as to ensure human safety and protect the environment. In the mil-itary arena, much of the military hardware must be able to adaquately absorb the energy released by the detonation of high explosives. In the transportation industry, vehicular crashworthiness during impact is another important area of concern. Also, offshore struc-tures and ships must be capable of sustaining the impact imparted by large wave forces. Lastly, the potentially destructive transient forces caused by the violent ground motion of J INTRODUCTION / 1.1 2 earthquakes must be acknowledged in the design of various structural systems throughout the world. In each of the examples cited above, the particular high intensity transient load may cause significant plastic deformation of the structure. Nonetheless, it is imperative that the damaged structure retain to a required degree a level of serviceability, so as to con-tinue to perform its intended function. The ability to accurately predict the deformation profile and the strain history throughout the structure is necessary in order to determine the ensuing damage to the various structural elements and the eventual failure of the en-tire structure. It thus behooves the structural engineering community to investigate both the geometric and material non-linear response of structures to transient dynamic loading conditions. To this end, much attention has been directed towards developing methods of analysis that adequately model the non-linear dynamic response of various structures. In an effort to reduce the complexity involved in modelling a complete structure many investigations have been limited to considering the response of individual structural com-ponents (i.e. cables, beams, plates, shells, etc.), which comprise the fundamental building blocks for entire structures. Realizing the immensity of the subject at hand this study is restricted to the non-linear dynamic analysis of thin ductile beams. The difficulties that arise from seeking a closed-form analytical solution to this non-linear beam problem, which involves time dependent finite deformations and inelastic material behaviour, has motivated many investigators to apply various assumptions and approximations to simplify their beam model. A considerable amount of fruitful research has been given over to evaluating the response of idealized rigid-plastic beams to impulsive loading . A classical example of this type of analysis is contained in the work of Symonds and Mentel [l] .f However, application of this rigid-plastic theory to such problems is only valid when the initial kinetic energy imparted to the structure is much greater than the elastic strain energy available within the structure. Furthermore, analytical solution procedures based on this f Numbers in brackets refer to references given at the end of this study. INTRODUCTION / 1.1 3 theory become unwieldy when the problem is extended to take account of various types of beam cross-section, strain-hardening, strain-rate dependence, and forcing functions of arbitrary distribution and time history. An alternate approach has been to accept the complexity of this non-linear beam problem, in a more general form, and then utilize a numerical prediction technique in order to obtain a solution. One such solution procedure is based on numerically evaluating the governing differential equations of motion, over the spatial domain of the problem, by application of the Finite Difference Method (FDM) and then employing a direct integration time operator over the temporal domain. This method has proven to be both accurate and efficient when applied to a particular problem. A noteworthy application of the FDM, to the problem at hand, is attributable to the work of Witmer et ai. [2]. In recent years, a more popular numerical solution technique to problems of structural dynamics has been to utilize spatially the Finite Element Method (FEM) in conjunction with a time-wise integration scheme. There exists a number of notable works that have applied this type of solution procedure in order to numerically simulate the non-linear transient response of beams [3-7]. As it pertains directly to the research undertaken in this study the solution technique developed by Wu and Witmer [4], merits special acknowledgement. The popularity of the FEM over the equally computationaly efficient FDM is attributable to the versitilty with which the F E M can be applied to structures with complicated geometric shapes, boundary conditions, and material properties. It is these features of the F E M which have made possible the development of general purpose finite element analysis computer programs. A number of published results for the non-linear beam problem have been made by using such general purpose programs as ADINA [8] and NONSAP [9]. An important feature in each of the aforementioned numerical procedures is the proper choice of a time integration scheme. There exists a number of well known time integration operators, that can be used in conjunction with a spatial solution procedure; namely, the Houbolt, Newmark, Wilson, and central-difference operators. A considerable INTRODUCTION / 1.2 4 amount of research has been given over to evaluating the efficiency and accuracy of these integration operators, when employed in non-linear problems [10,11]. An extensive literature review pertaining to the dynamic plastic response of structures has been compiled by Jones [12-15]. The interested reader will find in these references a comprehensive overview of much of the research that has been undertaken in this important area of structural engineering. 1.2 P U R P O S E A N D S C O P E In general terms, this study is concerned with the numerical simulation of Bernoulli-Euler beams undergoing large deformations and elasto-viscoplastic material be-haviour in response to transient dynamic loading conditions. The motivation behind this study has been to develop an efficient and accurate solution algorithm which can be em-ployed in various parametric studies relating to the non-linear response of slender ductile beams to air-blast pressure pulses. This will provide the analyst/designer with a numeri-cal tool that will aid in establishing design criterion for beams susceptible to such loading conditions. The solution procedure presented herein employes the Principle of Virtual Work in concert with the assumed displacement version of the Finite Element Method to formulate, in a variationally consistent manner, the differential equations of motion for the spatially discretized beam. These equations are then integrated time-wise by means of the central-difference method. The work undertaken in this study is restricted to considering only ini-tially straight beams, which possess at least one axis of cross-sectional symmetry and undergo planar (two dimensional) deformations. The beam kinematics are assumed to be of Bernoulli-Euler type, and are limited to producing small strains and moderately large displacements. The finite element modelling of this displacement field is adequately ac-complished by utilizing one-dimensional two node six degree-of-freedom beam elements. A INTRODUCTION / 1.2 5 versitile material modelling procedure is adopted in this study, which is capable of repre-senting constitutive behaviour ranging from linear elastic to elastic-plastic with kinematic strain-hardening. As well, visco-plastic behaviour is accounted for by means of a simple strain-rate dependent constitutive relationship. The effect of velocity dependent damping forces has not been incorporated into this analysis. Only the natural damping that arises from the plastification of the beam cross-section is taken into account. Finally, the flexibil-ity of this solution procedure is enhanced by being able to accommodate forcing function of arbitrary distribution and time history. Based on the solution algorithm developed in this study, a computer program has been developed to provide numerical predictions. In order to evaluate the accuracy and versatility of this solution scheme, numerical results have been generated for several examples manifesting geometric and/or material non-linear behaviour. All the examples selected are those for which either alternate analytical, experimental, or numerical solutions exist, so that comparisions can be made and conclusions drawn. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.2 6 CHAPTER 2 FORMULATION AND ANALYSIS OF THE PROBLEM 2.1 INTRODUCTION In this chapter a solution algorithm is presented which is capable of nu-merically simulating the transient dynamic response of Bernoulli-Euler beams undergoing geometric and material non-linear behaviour. Every effort has been made herein, to present a formulation which is mathematically consistent and faithful to the underlying principles of continuum mechanics, as found in the references of Fung [16] and Washizu [17]. It is hoped that this treatise elucidates the full scope of the required analysis, without being pedantic. 2.2 KINEMATICS OF A BERNOULLI-EULER BEAM A first step towards advancing a solution procedure is to establish the type and extent of beam deformations to be accounted for in this study. To this end, the strain-displacement relationship for a Bernoulli-Euler beam, undergoing large deformations, is derived. The terminology "large deformation" as used in this study, is employed with the proviso that even though the maximum lateral deflection experienced by the beam can FORMULATION AND ANALYSIS OF THE PROBLEM / 2.2 7 be large in comparision with its height it must still be small relative to the longitudinal dimension of the beam. The following assumptions are made regarding the kinematics of this non-linear beam formulation: 1. In the undeformed configuration the centroidal axis of the beam is straight and coincident with the z-axis. Furthermore, ( i — z) is a plane of symmetry for the beam cross-section and the applied loads. 2. The Bernoulli-Euler Hypothesis is valid. That is, normal sections remain plane, undistorted, and normal to the axis of the beam in its deformed configuration. 3. The Green strain measure is everywhere small, i.e. |e,- -| -C l.Also, the mid-surface strain and the square of the slope of the beam axis are small, i.e. \du/dx\ < 1 and {dw/dx)2 « 1. 4. The beam is slender so that terms of order h? can be neglected in the strain-displacement relationship, [h is a parameter representative of the cross-sectional height of the beam). As a direct consequence of Assumption 1 the deformation experienced by the beam is torsion-free and confined to the (x — z) plane [18]. Consider a segment of the beam, as illustrated in Figure 2.1, before and after deformation occurs. Each material point comprising the beam is specified by being A A A referenced to a fixed Cartesian coordinate system, for which i j , i 2 , and i 3 form the triad of unit base vectors. The position vector of a generic material point, °PQ, located on the undeformed beam axis is given by \ °T0 = (2.2.1) Another material point, °P , lying in the same cross-sectional plane as °P0, but not coinci-dent with the centroidal axis is located by the position vector °r = °r0 + y i 2 + z\z. (2.2.2) \ The superscript "(•) and subscript (-) 0 are used in this formulation to denote, respectively, a quantity in the undeformed configuration and that a quantity is associated with the middle surface of the beam. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.2 8 F I G U R E 2.1 Bernoulli-Euler beam kinematics. Upon deformation of the beam, the position vectors locating these material points, now denoted as P0 and P, are specified by the expressions: T O = °T0 + U O (2.2.3) r = °r + U, (2.2.4) where the u 0 and u are the resulting displacement vectors. Expressed in terms of the Cartesian base vectors these displacement vectors have the following form: u 0 = u(x)ii + w(x%, (2.2.5) u = u(x,y,z)i1 + v(x,y,z)i2 + w{x, y, z ) i 3 . (2.2.6) In equation (2.2.5) u(x) and w(x) represent, respectively, the transverse and longitudinal mid-surface, (z = 0), displacements of the beam. By utilizing the Bernoulli-Euler Hypothesis of Assumption (2), the position vector, r , can also be expressed as r = r 0 + zn, (2.2.7) FORMULATION AND ANALYSIS OF THE PROBLEM / 2.2 9 where n is the unit normal vector associated with point Pc on the deformed centroidal axis of the beam, and is determined from the relation, A = {*'o x h)/\*o x k\- (2.2.8) In the above equation and throughout the present chapter the prime denotes differentiation with respect to the variable x, namely, (•)' = d()/dx. With the introduction of Assumption (3) into the analysis Equation (2.2.8) has the following simple form: n = -u;'i1+i3. (2.2.9) Now the appropriate combination of Equations (2.2.2)-(2.2.4), and (2.2.7) yields the displacement field under which a Bernoulli-Euler beam deforms when the motion is confined to the (x — z) plane: u = uQ + z(n - i3). (2.2.10) Alternatively, this displacement field can be written in component form by combining Equations (2.2.5), (2.2.6), (2.2.9), and (2.2.10) to yield: u(x, y,z) — u(x) — zw'(x) v{x,y,z) = 0 (2.2.11) w(x,y,z) = w(x). With the displacement field now established attention turns to the resulting non-linear strain-displacement relationship. For a 3-dimensional continuum, referenced to a Cartestian coordinate system, the Green strain tensor is given by e -• = -t J 2 3ut- duj duk duk {i,j = 1,2,3). (2.2.12) dxj di{ dxi dxj The strain-displacement relationship given above can be applied to the Bernoulli-Euler beam kinematics considered in this study by identifying ux = u, = v, and u3 = w, as given in Equation (2.2.11), and acknowledging the restrictions imposed by Assumptions FORMULATION AND ANALYSIS OF THE PROBLEM / 2.3 10 (3) and (4). This results in the existence of only one non-vanishing Green strain tensor component, namely, €U = e + ZK, (2.2.13) where, du 1 (dw\2 J d2w dx2' In Equation (2.2.13) e represents the mid-surface strain experienced by the beam and K represents the curvature of the beam upon deformation. Hereafter, is denoted simply as e. The appearance of only one non-zero strain component is a direct consequence of Assumptions (l)-(4). It is to be noted, that the restricted class of beam kinematics adopted in this section will greatly simplify the analysis which is to follow. 2.3 VAR IAT IONAL E Q U A T I O N OF M O T I O N In this section the non-linear variational equation of motion for a Bernoulli-Euler beam is established. A Total Lagrangian Approach (TLA) is adopted, with the undeformed beam configuration referenced to a fixed Cartesian coordinate sytem, as illus-trated in Figure 2.1 . In order to obtain the equations of motion the Principle of Virtual Work is invoked, which for a 3-dimensional continuum, states that S W i n l - 6 W e x t = 0, (2.3.1) where <5W,n( and 6 W e x t are to be identified, respectively, with kinematically admissible variations in the internal virtual work and external virtual work of the continuum. For a deformable body undergoing finite deformations these quantities take the following form: ext = j dV, ay = J X-6adV + J T-5a, ^ (2.3.2) FORMULATION AND ANALYSIS OF THE PROBLEM / 2.3 In these equations °p, <V, and °S represent, respectively, the mass density, the volume, and the surface area of the continuum in its undeformed configuration, 5,y is the 2nd Piola-Kirchhoff stress tensor, e,y is the Green strain tensor, X is the body force vector, T is the vector of surface tractions, and u is the displacement vector. The symbol 6() is used to denote a generalized displacement variation in a given quantity consistent with the kinematical constraints imposed on the structure. A proof that the 2nd Piola-Kirchhoff stress tensor, and the Green strain tensor, e,y, form a suitable pair of energy conjugate stress and strain measures for use in the Total Lagrangian Approach of the Principle of Virtual Work is given in reference [19]. equilibrium deformation state to exist for a structure under a prescribed set of body forces, surface tractions, and boundary conditions. In order to extend this formulation to the dynamic realm d'Alembert's Principle is invoked. By so doing, the external virtual work expression becomes The contribution to the external virtual work arising from the self weight of the continuum has been ignored in the above equation. Specialization of Equation (2.3.3) to the beam theory of this study results in the following representation for the external virtual work: In the above equations °A and °L represent, respectively, the cross-sectional area and the length of the beam. The applied load vector, with the units of force/unit length, is denoted by f = L/to/wJ- Consistent with the T L A each of these quantities is given with reference to the undeformed configuration. Here, and throughout the remainder of this chapter, a superimposed dot denotes differentiation with respect to time, namely, (") = d{ )/dt. Note, The Principle of Virtual Work provides a necessary condition for a static (2.3.3) (2.3.4) FORMULATION AND ANALYSIS OF THE PROBLEM / 2.3 12 that the influence of rotatary inertia has been neglected in external virtual work expression of Equation (2.3.4). The justification for this is based on the fact that in general rotatary inertia constributes a smaller amount to the overall dynamic response of beams than shear deformation, which is being neglected in this study. Under the kinematical assumptions made in Section 2.2 it was determined that the only non-vanishing component of the Green strain tensor was en. Thus, the internal virtual work expression of Equation (2.3.2), when applied to a Bernoulli-Euler beam, reduces to SWint = J Sn5endV. (2.3.5) °v In the solution procedure which is to be developed it is advantageous to have the internal virtual work stated in terms of the true (Eulerian) stresses, <r,y. This is accomplished by application of the transformation equation relating the 2nd Piola-Kirchhoff stress tensor to the Eulerian stress tensor; In the above equation p is the mass density in the current configuration and 5,-y is the Kronnecker-delta.f For this equation to apply to the Bernoulli-Euler beam theory of this study it is necessary to make the identification u x = u — zw'. Then by consideration of the kinematical assumptions of Section 2.2, and assuming °p « p, which is valid in the presence of small strains, the above equation of transformation reduces to Su « cru, hereafter simply denoted as a. Thus, the internal virtual work expression for the beam becomes 6 W i n l = J a 5 e d V = j <r(& + Z5K) dV. (2.3.7) ay ay \ The Kronnecker-delta, <5,-y, is a second order tensor with the following functional form; FORMULATION AND ANALYSIS OF THE PROBLEM / 2.4 13 By introducing stress-resultants over the beam cross-section Equations (2.3.7) can be rewritten as 6Winl = j (JV& + M6K) dx, (2.3.8) with the resultant axial force, M, and resultant moment, M, given by M = J crdydz M = J azdydz. (2.3.9) °A °A A procedure to determine M and M is deferred to Section 2.5, where a non-linear con-stitutive relationship of the form a = /(e, e) will be considered. Finally, the variational equation of motion is obtained by equating the r.h.s. of Equations (2.3.4) and (2.3.8); j [M6e + M6K + (ftt - °p°Au) 6u + (fw- °p°Aw) 6w] dx = 0. (2.3.10) °L From this variational formulation of the problem it is possible to derive the complete set of governing differential equations of motion, and the natural boundary conditions, for a Bernoulli-Euler beam. However, the presense of geometric and material non-linearities in the problem makes the analytical solution of these partial differential equations, in general, impossible. Alternatively, a computationally efficient numerical solution strategy can be initiated by applying the assumed displacement version of the Finite Element Method to the variational formulation of the problem just presented. It is this approach that will be developed further in the next section. 2.4 F IN ITE E L E M E N T F O R M U L A T I O N OF T H E EQUAT IONS OF M O T I O N By means of the Finite Element Method the spatial domain of the beam will now be discretized into a number of subdomains, more commonly termed finite elements. Let NE represent the number of finite elements comprising the finite element mesh. In each finite element an approximate displacement field will be assumed. This displacement FORMULATION AND ANALYSIS OF THE PROBLEM / 2.4 14 F I G U R E 2.2 Typical beam finite element. field will be appropriately chosen so as to satisfy the conditions of completeness and com-patibility imposed by the variational formulation of the problem [20]. The deformation state of each beam element will then be characterized in terms of a finite set of generalized nodal coordinates. To this end, linear and cubic interpolating polynomials are chosen to approx-imate, respectively, the axial and transverse mid-surface displacements within each beam element. With reference to a local coordinate system as shown in Figure 2.2, this assumed displacement field has the following form: u{tj) = ao + diT} w(r)) = 60 + 6 ^ + b2t]2 + 63T? , where a0, alt 6 0,... ,6 3 represent as yet unspecified parameters. By choosing a one-dimensional two node finite element, as illustrated in Fig-ure 2.2, an appropriate set of generalized nodal coordinates, corresponding to the assumed FORMULATION AND ANALYSIS OF THE PROBLEM / 2.4 15 displacement field, is given by d,= dwi dw2 T (2.4.2) ft) * ft)'1 dV ' " " dt] J The subscript () e, as introduced in the above equation, associates a given quantity with a typical element in the finite element mesh. By expressing the unspecified parameters of Equation (2.4.1) in terms of these generalized nodal coordinates the assumed displacement field can be written as f u | = K 0 0 N2 0 O l where the shape functions Ni=i^^, are given by the following algebraic equations [21]: - - - f t ) MS NA = I? | \ - 2 * •= • [ - f t ) * ft)']' ' with le representing the undeformed length of a typical beam element. From the perspective of classical interpolation theory Ni &N2 are Lagrange polynomials and iV3, 7Y 4, J V 5 , & N6 are first order Hermitian polynomials for two node line elements. Their functional form is illustrated in Figure 2.3 . These shape functions constitute a set of linearly independent base functions by which the deformation profile of the beam element is modelled. With the assumed displacement field now established the resulting approxi-mate Green strain measure within each element can be given in terms of its nodal gener-alized coordinates. Substituting Equation (2.4.3) into Equation (2.2.13) yields the desired result ee = ee + CKE (2.4.5) FORMULATION AND ANALYSIS OF THE PROBLEM / 2.4 16 F I G U R E 2.3 Shape functions. where, with, B 3 = KE = B 3 d e dN, n n dN2 n n o, drj ' dr) ' ' cfy ' --32AT3 -d2N< -d2Nb -d2Ne For the spatially discretized beam the Principle of Virtual Work takes on the following form; (2.4.6) FORMULATION AND ANALYSIS OF THE PROBLEM / 2.4 17 where, U 6W?nt = J {^SH+MSK) dr} o Weext = I [(/. - YAu) 6u + (fw- °p°Aw) 5w) dr, o + / t u ^ i + fun^i + A*i#i + / « 2 ^ 2 + fwi^ + M2692. In the above external virtual work expression /„,-, & Mi (i = 1,2) are the inter-element nodal forces that arise due to the discretization of the beam. Note, that by ensuring that the assumed displacement field that is introduced into this variational formulation is at least CQ continuous in the axial displacement variable u and Ci continuous in the transverse displacement variable w then the above mentioned inter-element nodal forces will vanish in the summmation process of Equation (2.4.6).f The assumed displacement field of Equation (2.4.1) was chosen so as to provide this degree of inter-element compatibility. Since these quantities in the end provide no contribution to the virtual work expression they will hereafter be ignored in this formulation. Upon substituting the approximate displacement field, for a typical element in the finite element mesh, into the above external virtual work expression the following equation results; 6Weext = -5dTm ed e + 5 d J f e ) { 2 A 7 ) where, U m e = JyANTNdrj o u fe = JNTfdr,, o with N representing the matrix of shape functions given in Equation (2.4.3). The quantity m e is termed the consistent mass matrix for the beam element. Based on the assumed f The symbolism C m is used to signify that the derivatives through order 171 of a given displacement field are continuous. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.4 18 displacement field chosen for the beam element the consistent mass matrix has the following explicit form; °P°Ale 420 140 0 0 156 22/e All 70 0 0 140 0 0 156 -22/, All 0 0 54 -13/ e 13/e -3/2 (2.4.8) The quantity fe in Equation (2.4.7) is termed the external force vector for the beam element. The explicit form of this vector is dependent on the specified loading function f(*,0 = L/.(*.0./«(*.0J-Attention now turns to representing the internal virtual work for a typical element in the finite element mesh. To achieve this end, the approximate displacement field of Equation (2.4.3) is substituted into the internal virtual work expression to yield; where, = 6 d T p , ( d ) = 6 d T (pj + p$"d e ) , pf = y (B f t f+ BjAt) dt) (2.4.9) PeL = j BTB2M dr,. (2.4.10) (2.4.11) The quantity pe(d) denotes the internal force vector for a typical element. This vector takes account of the internal resistance developed within the beam element in response to the imposed dynamic loading. Even if the response causes geometric and/or material non-linearities to develop within the beam element this internal force vector continues to register the internal resistance that results to oppose the external loading conditions. The notation, pe(d), is adopted here to emphasize the functional dependence that this vector has on the element generalized coordinates, as a result of taking into account geometric non-linearities in the analysis. Evaluation of pe(d) is accomplished by resorting to numerical integration techniques, as will be explained later. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.4 19 At this point in the formulation of the problem it is advantageous to trans-form the element generalized coordinates, as given for a typical element by Equation (2.4.2), to a global set of generalized coordinates for the discretized beam as a whole. This is accomplished by means of a connectivety transformation of the form d e = C e D , (2.4.12) where D is the global vector of generalized coordinates and C e is the connectivity matrix associated with the finite element in question. Since the local coordinate system for each element is aligned with the global coordinate system for the beam the C e = 1 JVE are simply Boolean matrices. The Principle of Virtual Work, as expressed by Equation (2.4.6), has the following form when the global generalized displacements constitute the primary unknowns; NE 6BT (cTemeCeD + C j p f + C f p f L C e D - Cf f e ) = 0. (2.4.13) e=l Kinematically admissible variations in the components of the generalized displacement vec-tor are arbitrary and independent so that 5D can be eliminated from the above equation. Thus, the governing differential equations of motion, based on the spatial discretization of the Bernoulli-Euler beam by the Finite Element Method, are expressible in final form as f M ' D - f ' P ( D ) ='F, (2.4.14) where NE M = C T m e C e = Global Consistent Mass Matrix. e=l NE 'P(D) =^2 (CePe + C f P f L C e D J = Global Internal Force Vector. c=l NE 'F = ^ cefe = Global External Force Vector. e=l f The notation '(•) is introduced here, and nsed throughout the remainder of this chapter, to acknowledge the time dependent nature of various quanthes. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 20 In this formulation just presented, the assembly of the global equations of motion from the individual element matrices was accomplished by means of the connectiv-ity matrices Ce. This method of assembly is mathematically sound but computationally burdensome. In pratical application of the Finite Element Method efficiency is gained by use of the Direct Stiffness Method [22] for the assembly process. 2.5 MODELLING THE MATERIAL BEHAVIOUR OF THE B E A M A first step in solving the governing differential equations of motion presented in the previous section is to develop a procedure to determine the state of stress throughout the beam over the time domain. This will make possible the numerical evaluation of the element internal force vectors, tpe(d). For ductile materials subjected to high intensity transient loading condi-tions it is desirable to consider a constitutive law that incorporates elastic-plastic strain-hardening, and strain-rate sensitive material behaviour. It will be assumed in this study that the inelastic phenomenon of plastic yielding can be fully accounted for by considering only the uniaxial material behaviour of the beam, thus eliminating the need to incorporate a 2-dimensional yield criterion and associated viso-plastic flow rule into the analysis.!/ This assumption will greatly reduce the computational effort required to determine the state of stress, at any given time, within the beam. For the material modelling undertaken in this study it will be assumed that the constitutive law is expressible as a = f{e)g{e)t (2.5.1) where /(e) and g(e) represent, respectively, the strain-hardening and strain-rate behaviour of the material under consideration. Experimental evidence has shown that the uncoupling of strain-hardening and strain-rate effects, as expressed by Equation (2.5.1), is not realized f In accordance with Bernoulli-Euler beam theory the shear strain, "7, is everywhere zero, but from equilibrium considerations there may exist a shear stress distribution, T , within the beam. As a result of the assumption made above any interaction between (T and T in determining plastic yielding of the material has been neglected. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 21 by all materials. However, in the interest of mathematical simplicity the above functional form of the constitutive law is employed. Under quasi-static cyclic loading conditions a kinematic strain-hardening rule, which takes into account the Bauschinger effect, has been found to model the actual behaviour of ductile materials very well. Figure 2.4 illustrates the characteristics of this constitutive behaviour. This type of material strain-hardening can be efficiently accounted for numerically by utilizing the "mechanical sublayer model" [23,24].f With this material model the beam is assumed, for computational convenience only, to be composed of several material sublayers each of which experiences the same strain. By assigning different con-stitutive properties to each sublayer a composite behaviour can be obtained which exhibits all the essential characteristics of a ductile material undergoing cyclic loading conditions, as illustrated in Figure 2.4 . To apply this method, it is first necessary to approximate the experimental uniaxial tension stress-strain curve by (m+1) piecewise linear segments with defining coordinates ((xk,ek), (k = l , . . . , m ) , as illustrated in Figure 2.5a . At every point in the beam the material is assumed to consist of m equally strained "sublayers" of elastic, perfectly-plastic material with a common elastic modulus , i£ , and appropriately different quasi-static yield stresses, <jok = Eek (k = 1,..., m). The stress-strain curves for each of the sublayers is shown in Figure 2.56 . The state of stress, (rk, associated with the A;-th sublayer, at any time in the analysis is uniquely determined by consideration of the stress history and the current value of strain occuring at the point in question. By applying apppropriately chosen material weighting factors, Ak, to each sublayer the actual stress at any material point, corresponding to a strain e, is given by m * = ]£4b**(€) (2-5.2) k=l where, Ak = ^ , f The mechanical sublayer model is synonymous with the "overlay model", as it is sometimes identified as in the literature. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 22 F I G U R E 2.4 Uniaxial stress-strain behaviour of a ductile material under quasi-static cyclic loading conditions. with, El = E, Ek = °k~ ak~x (* = 2,...,m), E m + 1 = 0. ek - €k-i As given above, E^ represents the slope of the fc-th segment of the piecewise linear stress-strain curve. With a suitable number of elastic perfectly-plastic sublayers it is possible to model the actual uniaxial material behaviour of the beam to a high degree of accuracy. The procedure to obtain the state of stress at a given point in the beam at a given time in the analysis, employing the mechanical sublayer method, is as follows. For a given point in the continuum the k-th. sublayer stress ' - A <o"jfc a* * i m e * — At, and the strain increment A'e = 'e — l~Ate at time t are taken as known from the previous analysis. A trial stress value, denoted as TRIAL to the stress ta/s, at time t is computed by assuming that the material point in question undergoes a completely elastic deformation between time FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 23 F I G U R E 2.5 (a). Piecevtise linear approximation to the uniaxial tension stress-strain curve of Figure 2.4 . (b). The mechanical sublayer model is illustrated by assuming that the material is comprised of 3 "sublayers" of elastic perfectly-plastic material with a common elastic modulus and appropriately different yield stresses. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 24 t - At and t, so that TRIAL = t _ A t a f c + EAl€. (2.5.3) The correct value of lck is then obtained by consideration of the following conditional statements: If - aok < TRIAL < aok then t a k = TRIAL. (2.5.4) If TRIAL < -aok then lak = -crok. (2.5.5) If TRIAL > aok then l<rk = croJt. (2.5.6) This analysis is carried out in turn for each of the sublayers comprising the material model. The actual stress at the point in the beam under consideration at time t is then obtained by use of Equation (2.5.2). As demonstrated above, the mechanical sublayer method is computationally efficient and ideally suited to model material behaviour in the presence of cyclic loading conditions. Furthermore, it is very adaptable to the strain-rate theory given below. It has been observed experimentally that the yield stress of many materials increases with an increase in the applied strain-rate. Under high intensity dynamic loading conditions, which causes a continuum to experience high strain-rates, the instantaneous yield stress of the material, <ry, can be significantly greater than the quasi-static value, o~0. For strain-rate sensitive materials an approximate accounting of this phenomenon is possible by employing the Cowper-Symonds relationship [25]; = 1 + <7„ e D i (2.5.7) where D and p are experimentally determined material parameters. The functional form of this relationship is shown in Figure 2.6 for two common strain-rate sensitive structural materials: namely, 6061 —T6 aluminum alloy, and mild steel. An example of the idealized uniaxial stress-strain behaviour for such a strain-rate dependent material is illustrated in Figure 2.7 . FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 25 2.5-b i b 2 < OC 1.5 C/l (/) UJ a: i — (/> q i 0.5 y _ = 1 + D / / / / / Legend D = 6500/sec. p = 4 (mild steel) D = 40/soc. p=5 (6061-T6 aluminum) 0 1 i i i i i i i 11 II i i i i i 11II 0.001 0.01 0.1 1 10 STRAIN RATE e F I G U R E 2 . 6 Effect of strain-rate on the yield stress. I I 11 I— l-TTTTTT too 1000 F I G U R E 2 . 7 Idealization of strain-rate behaviour for an elastic-plastic strain hard-ening material. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 26 With a method to predict the state of stress point-wise within the beam now established, attention turns to presenting a procedure to evaluate the resultant axial force, and moment, *M, along the beam axis for each element of the finite element mesh. Evaluation of these integral expressions, as found in Equation (2.3.9), is best accomplished by resorting to a numerical integration technique. The optimal choice of an accurate and efficient numerical integration procedure is dependent upon the cross-sectional shape of the beam under consideration. In general, if the domain of integration permits Gauss-quadrature to be employed, it is a very good method to utilize. An indication of the efficiency and accuracy associated with this integration scheme is offered by the fact that it is capable of integrating exactly a polynomial of degree (2n—1) with only n Gauss evaluation points. Table 2.1 offers a brief summary of some of the details of Gauss-quadrature. In order to illustrate the procedure it will be assumed that the beam cross-section is rectangular; of width b and height h. The integrals of Equation (2.3.9), have the following form when rewritten with respect to the natural coordinate £ = 2^//i; l 2 1 ^ = b ± j t c d i tft-^L J tai d^ (2.5.8) - l - l Utilizing Gauss-quadrature and the material modeling procedure presented in this sec-tion it is possible to numerically evaluate the above integrals, resulting in the following expressions: I L 9 / m y=i \k=i j=i \k=i where Wj is the weight factor associated with the j'-th depthwise Gaussian evaluation point. With a procedure now established to determine the resultant axial force/^/, and moment, lM, pointwise along the beam element axis it is possible to evaluate the FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 27 T A B L E 2.1 One dimensional Gauss-quadrature. I = +1 n -1 »=i weight factor where W,- = & = coordinate of the tth integration point n = total number of integration points En = error = 0{d2n f / d?n) n i 1 (linear) 1 0 2 2 (cubic) 1 2 -+ 1 + 1 3 (quintic) 1 2 3 0 + v/l5/5 ->/» 5/5 8/9 5/9 5/9 4 (septimal) 1 2 3 4 + 0.861 136 31 -0.861 136 31 + 0.339 981 04 -0.339 981 04 0.347 854 85 0.347 854 85 0.652 145 15 0.652 145 15 internal force vector tp e(d) for each element. Again, employing Gauss-quadrature this vector, as given by Equations (2.4.9)-(2.4.11), can be evaluated approximately as lVe(d) = \JZWi (tf(r)i)%(rii) +B^r,i)B2(vi)tde%(Vi) + B ^ r ? , ) , (2.5.10) «=i where W,- is the weight factor associated with the i-th spanwise Gaussian integration point, located at tj,- = le(£j + l)/2 along the beam element axis. From the above discussion regarding the numerical integration of the internal force vector it is seen that an integration grid of size p x q is required to be specified, where p and q represent, respectively, the number of spanwise and depthwise integration points within the element. A determination of the number of depthwise integration points, to FORMULATION AND ANALYSIS OF THE PROBLEM / 2.5 28 accurately evaluate the resultant axial force and moment, will first be considered. For a beam of rectangular cross-section undergoing purely elastic deforma-tions only two depthwise Gaussian evaluation points are required to exactly numerically integrate these quantities. However, as plastification of the material extends through the beam cross-section a higher order of Gauss integration is required to capture the non-linear distribution of stress. For thin beams, of rectangular cross-section, Wu and Witmer [26] found that 4 depthwise Gaussian points were sufficient to give an accurate representation of this non-linear stress distribution. This same order of Gauss integration will be em-ployed in this study for beams of rectangular cross-section. It is also desirable to consider beams which do not have rectangular cross-sect ions but nonetheless are commonly used in structural applications. Of particular interest is the non-linear transient response of I-beams and T-beams. For beams possessing such cross-sectional shapes, a 4 point Gauss rule is used for numerical integration through the web and a simple mid-point integration rule in the flange(s). For purely elastic deformations this integration scheme will again provide the exact value for the axial force and moment. However, as plastification of the cross-section occurs this integration scheme is not completely capable of reproducing the stress distribution through the thickness of the beam. In particular, it is to be noted that yielding of a flange will only take place when the stress level at the mid-point of the flange reaches |cr0|, upon which the entire flange is assumed to have yielded. Though this assumption introduces some error into the integration scheme it is believed to be minimal since most commerically fabricated I-beams and T-beams have relatively thin flanges. Next, a determination of the appropriate number of spanwise Gaussian eval-uation points is undertaken. The quantity which imposes the strictest condition on the number of spanwise integration points is clearly seen to be l p ^ L . With increasing finite element mesh refinements it is expected that the resultant axial force within any element, under static loading conditions, will approach a constant value. In order to ensure that the numerical integration scheme can exactly represent l p N L when the resultant axial force FORMULATION AND ANALYSIS OF THE PROBLEM / 2.6 29 is constant a minimum of 3 spanwise Gaussian evaluation points are required within the element.f For thin beams, Wu and Witmer [27] found that 3 spanwise Gaussian evalu-ation points, in the numerical integration scheme of the internal force vector, produced sufficiently accurate results. This same order of Gauss integration will be employed in this study for the spanwise evaluation of the element internal force vectors. Note, that by application of the Finite Element Method and the material modelling procedure just presented, the numerical simulation of the dynamic response of a Bernoulli-Euler beam is reduced to solving a system of 2nd order differential equations, which are functions of time, as expressed by Equation (2.4.14). 2.6 SOLUTION PROCEDURE OVER THE TIME DOMAIN The next step in the solution procedure is to solve the differential equations of motion over the time domain. Solution techniques to problems of structural dynamics are basically of two types; modal superposition and direct integration. Modal superposition, however, is most often restricted to use in linear problems. Direct integration, on the other hand, has been effectively used in both linear and non-linear problems alike. Of the many different direct integration methods developed the central difference method is employed in this study. Many investigators have concluded that the central difference method is the most efficient time integration operator, when applied to problems that involve short durations of loading, which result in the structure being excited throughout its frequency spectrum [28,29]. The positive aspects and the limitations associated with this particular numerical method of time integration will be acknowledged as the underlying theory of this method is presented and applied to the problem under consideration. By application of the Finite Element Method, in Section 2.4, the spatial domain of the beam was discretized into NE finite elements and the following 2nd order differential equations of motion resulted; M * D + ' P ( D ) = J F . (2.6.1) f In point of fact, 3 spanwise Gaussian evaluation points yields an exact representation of the internal force vector when there exists a linear variation in the resultant axial force along the beam element. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.6 30 In using the central difference method to solve the above equations the temporal domain of the problem is discretized into equal time intervals of duration At and an approximate solution is sought by employing this integration scheme at each of the time intervals, (0, At, 2 At, ...,t, t+At,T). The central difference method affords a way to express the acceleration vector of the generalized coordinates in terms of the generalized displacement vector, as shown by the following equation; 'D = ^ ± 1± + 0 ( A t ) 2 , (2.6.2) where O(At)2 indicates that the truncation error in this approximation is of the order {At)2. It is advantageous at this point to relate displacements and displacement increments in the generalized coordinates by the following expressions: 'AD = 'D - ' ~ A ' D (2.6.3) 'D = ° D + A t A D + . . . + ' " A ' A D + 'AD. (2.6.4) In terms of displacement increments, Equation (2.6.2) can be approximated as „ t + A t A D — 'AD 'D = — 2 . (2.6.5) Based on the central difference quotient of Equation (2.6.5) the dynamic equations of motion can be written as t + A ' A D = 'AD + A t 2 M _ 1 ('F - 'P(D)) . (2.6.6) Equation (2.6.6) establishes a reccurence relationship whereby the unknown displacement increments at time t + At are given in terms of known quantites, by considering dynamic equilibrium of the beam at time t. This type of solution method is appropriately called an explicit integration scheme, and clearly is computationally simple. To start this algorithm it is necessary to have an expression that relates the displacement vector, evaluated at time At, to known quantites given at time t = 0. This is accomplished by considering a Taylor series expansion of the displacement vector about time t = 0, A ' A D = At°D + }-(At)2oD -f O(At) 3. (2.6.7) It FORMULATION AND ANALYSIS OF THE PROBLEM / 2.6 31 Substituting Equation (2.6.1), evaluated at time t = 0, into the above equation results in A ' A D = At°D + ^ ( A * ) 2 M _ 1 ° F , (2.6.8) it under the assumptions that ° D = 0 and the initial stress distribution in the beam is zero. The latter assumption implies that ° P ( D ) = 0. With the integration scheme now complete, the numerical simulation of the dynamic response of the beam can be carried out over the entire time domain by employing Equations (2.6.4), (2.6.6), and (2.6.8) in conjunction with the material modelling procedure and the finite element formulation developed in the previous sections of this chapter. The computational simplicity associated with use of the central difference method is clearly evident from consideration of the reccurence relationship of Equation (2.6.6). By application of this equation the deformation increments in the vector of gener-alized coordinates, occuring at any time step, are directly obtainable by performing simple arithmetic operations on matrices of known value. Thus, to advance the solution one time step relatively little computational effort is required when compared with that needed by an implicit method of integration which requires the solution of a non-linear system of equations at each time step.f The computational efficiency of this time integration algorithm is further enhanced by utilizing a diagonal mass matrix in place of the consistent mass matrix. With a diagonal mass matrix the inversion process required in Equations (2.6.6) and (2.6.8) becomes trivial. The following diagonal mass matrix has been proposed in order to model the translational inertia of the beam element more simply; °p°Ale m « = ~420" 210 210 all 210 210 all (2.6.9) where a = 17.5 is a gradient scaling parameter [31]. Another important consequence of coupling a diagonal mass matrix with the central difference method is the improved f A complete solution algorithm for an implicit integration scheme can be found in reference [30]. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.6 32 accuracy of the numerical integration over the time domain, as compared with a central difference formulation that utilizes the consistent mass matrix [32]. The improved accuracy associated with this coupling manifests itself because whereas the diagonal mass matrix under estimates the frequency response of the structure the central difference integration scheme distorts the frequencies upward. These compensating errors assist in ensuring that an accurate representation of the dynamic response of the structure is achieved. The main shortcoming associated with use of the central difference method is that the algorithm is conditionally stable, even when the response of the structure is restricted to being linear. This implies the existence of a critical or limiting time step size, denoted as At c r , which when exceeded causes the displacement vector to grow without bound as the analysis proceeds. It is well known that the stability criterion for the central difference method is directly related to the magnitude of the largest natural frequency, (jJmax, occuring in the discretized system, such that,t Atcr = 2/umax. (2.6.10) The form of this stability criterion holds true whether the structural response is linear or non-linear. Of course, when the response is non-linear wmax represents the maximum natural frequency of the non-linear system. However, the presence of geometric and/or material non-linearities in the response generally does not introduce any severe distabilising mechanisms into the integration scheme. So that from a practical perspective, an estimate to A£ c r is obtained by reducing the result determined by using Equation (2.6.10), based on a linear eigenvalue analysis, by 10—20%. At this point, one may think it necessary to perform a complete eigenvalue analysis on the discretized beam in order to obtain the stability limit. This is not the case. An upper bound on the highest natural frequency of the structure, which results in a conservative estimate of the critical time step, can be simply obtained by determining the maximum unconstrained natural frequency of the smallest element f A proof establishing the stability criterion for the central difference method, when used as the integration scheme for a linear structural system, is given in Appendix A. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.6 33 in the finite element mesh.f This useful result is established in a theorem contained in reference [33]. The procedure is made even more attractive in that the eigenvalue analysis for the Bernoulli-Euler beam element can be carried out analytically. The linear eigenvalue problem, which is of the following form; (ke - u 2 m e ) d = 0, (2.6.11) is applied to the smallest element in the finite element mesh. In the above equation k e represents the linear stiffness matrix for the element, and is obtained by evaluating the following integral expression; k e = j (Bl + cB3)T E ( B j + <rB3) dV. (2.6.12) Upon integration, the element stiffness matrix has the explicit form; k e n All 0 0 12/ 6//, Alt -Al\ 0 0 Al\ 0 12 -6IL 0 -12/ 6//e 2IH 0 0 12/ -6//e All2 (2.6.13) where A and / represent, respectively, the cross-sectional area and the moment of inertia of the beam element. In turn, the linear eigenvalue problem is solved by employing first the consistent mass matrix of Equation (2.4.8) and then the diagonal mass matrix of Equation (2.6.9). Combining these results with the stability criterion of Equation (2.6.10) yields conservative bounds on the critical time step size when either a consistent mass matrix or a diagonal mass matrix is utilized in the integration scheme: Consistent Mass Matrix -A * , Diagonal Mass Matrix A t = f ' e / C v ^ c J , if/ e>10v /7r G; \ l2/ (lOv/21 cjr rG), otherwise. = fh/cL, \ile>As/3rG; \ l 2 l (4\/3cjr r G ) , otherwise. (2.6.14) (2.6.15) \ The maximum unconstrained natural frequency refers to the numerically largest frequency obtained, based on an eigenvalue analysis of the beam element with all of its degrees of freedom unconstrained. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.6 34 In the above equations, ci = \JE/p is the speed at which longitudinal waves prop-agate along the beam axis, and rG — \JljA is the radius of gyration of the beam cross-section. Typically, for most structural metals (i.e. steel, alumimum, titanium), cL « 200,000in/sec. From the stabilty analysis performed above a number of important results are immediately apparent. First, it is to be noted that using a diagonal mass matrix in conjunction with the central difference method optimizes the critical time step size, thus improving the efficiency of the integration scheme. This benefit, coupled with the aforementioned improved accuracy associated with this particular combination validates employing the diagonal mass matrix in concert with the central difference method in the solution algorithm. Next, consideration of Equation (2.6.15) reveals that the limiting time step size is a function of the spatial discretization of the beam model. Thus, refining the finite element mesh requires a more restrictive time step to be employed. In particular for h ^ Ay/Srg any mesh refinement results in a time step reduction proportional to the square of the element length. For this reason extreme mesh refinements in any one portion of the spatial domain is undesirable as the critical time step size for the entire problem will be governed by the size of these elements. A further step to ensure stability of the integration scheme over the temporal domain of a problem is achieved by performing an energy balance calculation at each time step of the analysis [34]. Any significant loss in energy balance ( 5% or more ) indicates an instability. The critical time step size given by Equation (2.6.15) is based on an assumed constant beam stiffness. This is no longer the case when geometric and/or material non-linearites are present in the structural response. In turn, Atcr will no longer remain constant. Exceding Atcr for a portion of the analysis will cause a temporary instability in the integration scheme. If this coincides with energy dissipation due to plastic work the instability may be arrested. In such a situation the solution does not "blow up", but is nonetheless erroneous. This arrested instability is easily detected by performing an energy balance calculation at each time step. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.7 35 A procedure to implement the energy balance check into the solution algo-rithm is now presented. At time t + At in the analysis the internal energy of the discretized beam is given by the expression l+"Wint = lWint + 1 ADr ('P + < + A ' p ) . (2.6.16) Equation (2.6.16) represents a trapezoidal integration of the non-linear stress-strain curve for the finite element model of the beam. Similarly, the external work is given by the expression l+*lWext = lWext + iAD r (<F + <+A'F) . (2.6.17) The kinetic energy is simply given by KE=^DTMt). (2.6.18) An energy balance is then ensured at time t + At if Wint + KE- Wext < TOL, (2.6.20) where TOL is an acceptable tolerance. In problems of structural dynamics for which the duration of loading is typi-cally of the order of the fundamental period of the structure the limiting size of the critical time step of the central difference method brings into question the applicability of this integration scheme. However, for dynamical problems in which the duration of loading is small, as is the case in the blast loading of structures, the central difference method is a computationally efficient time integration algorithm. Furthermore, since the duration of a blast load is very often extremely short a small time step is advantageous in order to model all the characteristics of the loading history. 2.7 SUMMARY This chapter concludes with a flowchart , given in Table 2.2, which outlines the required steps in the solution algorithm. FORMULATION AND ANALYSIS OF THE PROBLEM / 2.7 T A B L E 2.2 Flowchart for the Solution Algorithm. 36 (1) Input the geometric and material properties of the beam. Define the finite element mesh. (2) Set the initial conditions: °Dt ° F ; set counter n = 1, t = At. (3) Assemble and then invert M . (4) Determine A * D by Equation (2.6.8). (5) Determine the internal force vector by looping through all the elements in the mesh: (a) determine strains le by Equation (2.4.5) (b) determine stresses by the material law (<7 = f(e)g(e) (c) determine the element force vector, tpe(d), by Equation (2.5.10) (d) assemble ' P ( D ) (6) Determine the external force vector ' F . (7) Check the energy balance by Equation (2.6.20). If 6 is exceded then stop. (8) Set the counter: n—*n+l\t—*t + At. (9) Determine * D by Equation (2.6.6). (10) If t = T stop; otherwise go to step (5). NUMERICAL EXAMPLES / 3.2 37 CHAPTER 3 NUMERICAL EXAMPLES 3.1 I N T R O D U C T I O N In this chapter the ability of the present numerical method to predict the non-linear transient response characteristics of dynamically loaded beams is examined through a series of examples. These examples are chosen with the intent of testing the accuracy and efficiency of the solution algorithm. As well, these examples illustrate the potential of this numerical method in the context of practical application. 3.2 E X A M P L E EXHIB IT ING M A T E R I A L NON-L INEAR B E H A V I O U R First, the dynamic response of a simply supported beam, of rectangular cross-section, under a uniformly distributed step load is studied. The material comprising the beam is modelled as being elastic perfectly-plastic. For this particular example small displacement theory is invoked so that the only type of non-linearity present during the response results from the material behaviour. The dimensions and material properties of the beam, as well as the loading history, are given in Figure 3.1 . Due to the symmetry NUMERICAL EXAMPLES / 3.2 38 of the beam configuration and the imposed loading only one half of the span is required to be modelled in the numerical solution procedure. The finite element idealization of this portion of the beam is comprised of 6 elements, each of equal length, as illustrated in Figure 3.1 . For this finite element mesh the stability criterion on the temporal integration scheme yields a critical time step size of magnitude A t c r = 7.72/isec. Since le < 4\/3rG the critical time step is determined by the speed of flexural waves propagating in the beam, and thus A f c r a / g . It is to be noted that the solution algorithm continued to give a stable solution up to a time step size of 0 .99A£ c r . The fundamental period of this beam, denoted as !Ty, is 4900 fisec for linear elastic undamped vibration. A measure of how restrictive the stability criterion is for the central difference method, as it applies to this particular problem, is obtained by observing that the critical time size can be expressed as Ty/635. The magnitude of the uniformily distributed step pressure is 0.75Po, where Pa is the static collapse load of the beam.f Figure 3.2 shows the transient response of the midspan deflection of the beam, expressed in non-dimensional terms, as predicted by the study4 The deflection ratio tw/A gives the transverse displacement of the midspan of the beam normalized with respect to A , the static elastic deflection of the midspan resulting from the load 0.75Fo. The temporal measure of the reponse has been normalized with respect to the fundamen-tal period of the beam. For comparitive purposes, the linear elastic small displacement response, as predicted by this solution procedure, is also shown in Figure 3.2 . This numer-ically generated linear response is in excellent agreement with the modal solution given by elementary dynamic beam theory. The above statement is substantiated by noting that w(L/2,t) - wFE{L/2,t) A where w{L/2,t) and wpE(L/2,t) represent, respectively, the modal solution and finite element solution for the transverse midspan deflection time history of the beam. For each j " T h e s t a t i c c o l l a p s e l o a d , Pa, f o r a s i m p l y s u p p o r t e d b e a m i s d e f i n e d a s t h e l o a d i n t e n s i t y ( f o r c e / u n i t l e n g t h ) w h i c h c a u s e s a p l a s t i c h i n g e t o o c c u r a t t h e m i d s p a n o f t h e b e a m . J T h e r e s p o n s e o f t h e b e a m , a s g e n e r a t e d b y t h e c o m p u t e r p r o g r a m , i s s p e c i f i e d b y d a t a g i v e n a t d i s c r e t e i n t e r v a l s o v e r t h e t i m e d o m a i n o f t h e p r o b l e m . I n o r d e r t o r e p r e s e n t t h e r e s p o n s e a s a c o n t i n u o u s c u r v e L a g r a n g e i n t e r p o l a t i n g p o l y n o m i a l s h a v e b e e n e m p l o y e d b y t h e p l o t t i n g r o u t i n e s . < 0.025 for 0 < t < Tf, NUMERICAL EXAMPLES / 3.2 39 6 elements in half span 0.75P„ T t V\ 1 •—?—?—*—* f t t t t t t t - I 0.7SP.-Step Load L = 30.0 in 6 =1.0 in (width) h — 2.0 in £ = 3.0x 107 psi <r0 = 5.0 x 10* psi p = 0.733 x 10"3 lb sec2/in* P„ = 2a0b{h/L)2 = 444.4 lb/in A = 0.1758 in • t, sec. F I G U R E 3.1 3.5 3-2.5-2-1.5-0.5 Simply supported beam subjected to uniformily distributed step load. / DYNAMIC ELASTIC-PLASTIC SOLUTION / w DYNAMIC ELASTIC SOLUTION \ / \ / 4-• » y STATIC ELASTIC/ * \ / — 1 — 1 I 1 I 1 I 0.2 0.4 0.6 0.8 \ SOLUTION j \J * \ ft \ 2.2 F I G U R E 3.2 1 1.2 1.4 1.6 1.8 2 TIME RATIO t/T Dynamic small displacement elastic and elastic perfectly-plastic re-sponses of a simply supported beam. NUMERICAL EX.\MPLES / 3.2 40 of the response curves given in Figure 3.2 approximately 1500 time steps were required to numerically simulate the two cycles of beam vibration, based on a time step size of 0.90Atcr. One can observe the effects that the non-linear material behaviour has on the response of the beam from Figure 3.2 . In particular, the initial period of vibration has increased to approximately 1.287y, due to the softening of the structure. As well, the peak midspan deflection has increased from 1.5A, for the linear elastic solution, to 2.2A for the elastic-plastic solution. These characteristics of the non-linear response of the beam are attributable to the significant plastification of the beam cross-section resulting from the imposed load intensity. To appreciate more fully the extent of material non-linearity present during the response of the beam the strain history of a particular point on the beam, located at (x,z) = (13.75in, 1.0in), is given in Figure 3.3 . This monitoring point is located on the outer surface of the beam at the mid-point of the 6th element, which is adjacent to the midspan of the beam. From the stress-strain relationship it is determined that incipient plastic straining of the material occurs at e0 = 0.00167. Figure 3.3 confirms that at this monitoring point along the beam a significant amount of plastic strain occurs during the first cycle of the response, with the peak strain reaching a magnitude of 6e0. In order to comment on the ability of this formulation to model the dynamic response of this beam, in the presence of non-linear material behaviour, comparison is made with the numerically generated response obtained by Bathe et ai.[35], as illustrated in Figure 3.4 . Bathe obtained a plane stress solution using the general purpose finite element program NONSAP. Six 8-node isoparametric elements were used to model one quarter of the beam. Plasticity was accounted for with the von Mises yield criterion. Time integration was performed using Newmark's method with a time step size of 50 fisec. From examination of Figure 3.4 it is seen that the present study predicts a stiffer beam response, resulting in a smaller peak midspan deflection, ( » 6% less), and a shorter initial period of vibration, ( « 4% less), than the full plane stress finite element solution obtained by Bathe. NUMERICAL EXAMPLES / 3.2 41 0.010 0.008 0.006H < on r— O 0.004-0.002-0.000 0 F I G U R E 3.3 3.5 z g p L d Q F I G U R E 3.4 / / / / / / / -A YIELD STRAIN LEVEL 0.2 0.4 0.6 0.8 1 TIME RATIO t/T Strain history of a simply supported beam. i 1.2 Legend BATHE et al. - NONSAP [35] THIS STUDY 0.4 0.6 0.8 TIME RATIO t/T 1.2 1.4 1.4 Comparison of numerically predicted response histories of an elastic plastic simply supported beam. NUMERICAL EXAMPLES / 3.3 42 It is suspected that the discrepency between the two solutions is mainly attributable to the inability of the present study to take into account the influence that the shear stress T has on the plastic yielding of the beam material. Liu and Lin [36] in their study of this problem, found that the initial peak midspan deflection of the beam was approximately 9% greater when the interaction between a and T was accounted for in the analysis (through the von Mises yield criterion and the Prandtl-Reuss flow rule) in comparison with a solution solution procedure that only took account of the axial stress, a . 3.S EXAMPLE EXHIBITING GEOMETRIC NON-LINEAR BEHAVIOUR A clamped-clamped beam, of rectangular cross-section, acted upon by a concentrated step load located at the midspan of the beam is now considered. The material comprising the beam is assumed to remain elastic throughout the analysis, so that only non-linear geometric effects are manifested during the response. The geometry and the material properties of the beam, as well as imposed loading history, are given in Figure 3.5 . Symmetry of the beam configuration and the load distribution permits only one half of the beam to be considered. A uniform finite element mesh, comprised of 5 elements, is chosen to model the spatial domain of the beam to the midspan. For this finite element mesh the resulting critical time step size on the central difference integrator is 5.77/usee. Figure 3.6 shows the non-linear transient response of the beam, as predicted by this study. The non-linear geometric effects are clearly evident in the response history of the midspan of the beam since the peak deflection is of the order of 6 beam depths. As a result of these large deflections and the slenderness of the beam cross-section the presence of membrane forces are expected to play a dominant role in the non-linear transient behaviour of the beam. The time for the beam to oscillate through the first cycle, as shown in Figure 3.6, is approximately 2300 fisec. Based on a modal solution for the linear response of the beam the fundamental period of vibration is 9056 p. sec. This confirms that the beam does stiffen considerably in the presence of these large deflections. NUMERICAL EXAMPLES / 3.3 43 h 4 5 elements in half span | P = 640/ft r i \ \ r •1/2- L/2 P = 640/6. Step Load t, see L - 20.0 in 6 =1.0 in (width) n=l /8m £ = 3.0x 107jwi p = 0.733 x 10 -3 /& see7/in* F I G U R E 3.5 Clamped beam subjected to a point step load at midspan. 1.2-1.1-1-.c 0.9-r— 1 . 1 0.8-EME 0.7-o • < 0.6-Q_ • Ul 0.5-O AN 0.4-Q_ -(/) O 0.3-2 0.2-0.1-0-Legend MONDKAR AND POWELL |37] YA.Np_ANP..s_AJGAL_ __[38j___ > THIS STUDY 1000 2000 3000 TIME (/zsec.) 4000 5000 F I G U R E 3.6 Large deflection elastic dynamic response of a clamped beam. NUMERICAL EXAMPLES / 3.4 44 The accuracy with which this solution algorithm takes account of non-linear geometric effects is checked by making comparisons with the numerically generated re-sponse given by Mondkar and Powell [37]. In their study of this problem, five 8-node plane stress elements were used to model one-half of the beam. Integration over the time do-main was performed using Newmark's method, with a time step size of 50 fisec. Very good agreement is noted between these two numerical predictions, as is illustrated in Figure 3.6 This particular beam problem has also been investigated by Yang and Saigal [38]. In their study, six 1-dimensional beam elements were used to model one-half of the beam. Time integration was performed using Newmark's method, with a time step size of 10 usee. Again, good agreement is noted between their numerical prediction and the results of this study, with the exception of two unexplained variations in the response histories occuring at approximately 2000 fisec and 4000 fisec, as illustrated in Figure 3.6 . From this example it is concluded that the solution algorithm of this study is able to accurately take account of non-linear geometric effects, with the proviso that the deformations are in accordance with the restrictions imposed in Section 2.2 . 3.4 E X A M P L E S EXH IB IT ING G E O M E T R I C A N D M A T E R I A L NON-L INEARIT IES (a) Impulsively Loaded Beam First, a clamped-clamped beam, of rectangular cross-section, subjected to an impulsive load centered over a portion of the midspan of the beam is studied. This particular example is chosen as both geometric and material non-linearities are present in the response and it affords the opportunity to compare the numerically generated results with experimental data obtained from actual test beams. Witmer et ai.[39] discuss the experiments performed on these explosively loaded test beams. The actual beams were comprised of mild steel, which exhibits both strain-hardening and strain-rate sensitivity. The static stress-strain curve for this material has been fitted by two straight lines: one NUMERICAL EXAMPLES / 3.4 45 for the elastic portion and a second approximating the plastic region up to 4% strain. This implies using two material sublayers in the mechanical sublayer model. The strain-rate model of this material, based on Equation (2.5.7), produces a stress strain-rate curve similar to the one given in Figure 2.7 . The material and geometric properties of the beam are given in Figure 3.7 . Symmetry of the beam configuration and the initial velocity distribution requires only one half of the beam span to be modelled. The finite element idealization of this portion of the beam is made up of 10 elements. This discretization of the beam leads to a critical time step size of 2.39 fisec for the central difference integrator. The critical time step can be expressed alternatively as a fraction of the fundamental period as Ty/1610. Attention now turns to the numerical modelling of the actual explosive load which was applied to the test beams. In the experiment a thin layer of high explosive covered the entire width of a portion of the midspan of the beam . Since the duration of the pressure pulse applied to the beam from the detonation of the high explosive was typically less than 5 fisec, which is extremely short compared with the overall response time of the beam, the forcing function is treated as an impulsive load. This in turn is expressible as an initial velocity distribution, as illustrated in Figure 3.8a . The finite element representation of this pure impulse, hov/ever, is not correctly given by specifying a uniform initial velocity distribution at nodes 9 through 11, as given in Figure 3.8a . Because of the way in which the beam elements are forced to deform, as dictated by the assumed displacement field, the initial energy imparted to the beam is over estimated by the finite element modelling of this uniform velocity distribution, as shown in Figure 3.86. A more accurate finite element representation of this impulsive load is offered by the velocity distribution given in Figure 3.8c . Figure 3.9 illustrates experimental-numerical comparisions of the midspan deflection time history for the mild steel beam. Discussion will initially be limited to the numerically generated response, as given by the solution procedure of this study. Two NUMERICAL EXAMPLES / 3.4 46 L h i V 10 elements in half span ^ | | | < > n = 5020 in/see * In * L, k L, = 4.035 in Z 2 = 0.965 in 6 =1.0 in (width) h = 2.0tn »- E-X 3.0 x 107 psi 3.85 x 10s psi (tangent modulus) 1.16 x 10* psi D = 40.4/sec ,6 = 5 (strain-rate parameters) P = 0.733 x 10 -316aecV»'n* F I G U R E 3.7 Impulsively loaded clamped beam. *tt = 5 0 2 0 » n / s c c ttg = °Wl0 = °WU - °W Node# 9 10 11 (o) Representation of the explosive load as a uniform velocity distribution Over estimation of the initial velocity ^ — u v e r esi: a : (6) Finite element modelling of the velocity distribution of Figure 3.8a \atdz), ' (c) A more accurate finite element modelling of the initial velocity imparted to the beam F I G U R E 3.8 Finite element modelling of the impulsive load. NUMERICAL EXAMPLES / 3.4 47 different material models for the beam are considered for comparitive purposes: namely, an elastic-plastic strain-hardening material model (EL-SH) and an elastic-plastic strain-hardening strain-rate sensitive material model (EL-SH-SR). From Figure 3.9 it is observed that for each beam non-linear geometric effects are manifested in the response history of the midspan deflection. It is to be noted, however, that the EL-SH-SR beam behaves considerably stiffer in response to the imposed loading conditions than the EL-SH beam. The peak midspan deflection for the EL-SH-SR beam is approximately 0.6 in compared with 0.75 in. for the EL-SH beam. A greater appreciation of the effect that strain-rate has on the response of the beam is afforded by examining the strain history of a particular material point, located at (x,z) — (4.76in, 0.0625in), as offered in Figure 3.10 . This monitoring point is located on the outer surface of the beam at the mid-point of the 10th element, which is adjacent to the midspan of the beam. From the material properties given for the EL-SH beam plastic straining of the material fibres occurs at ea = 0.000389. Consideration of Figure 3.10 indicates that the EL-SH beam experiences significant plastic straining under the imposed loading conditions. In fact, the peak strain registered at the monitoring point is approximately 100 ea. The peak strain experienced by the EL-SH-SR beam is considerably less at approximately 65 e0. From an examination of the strain history for the EL-SH-SR beam it is observed that the average strain-rate experienced at the monitering point over the first 100 fisec is approximately 230/sec. By consideration of the strain-rate relationship of Equation (2.5.7) this indicates that the yield stress is raised to 2.42 ua during this portion of the response. This brief analysis of the strain-rate behaviour of the material offers an indication of the important influence that this phenomenon has on the response of structures which are subjected to high intensity impulsive loads. Note, for both the EL-SH beam and the EL-SH-SR beam the strain measure during the reponse remains very much less than one in absolute value, in accordance with the kinematic restrictions imposed in Section 2.2 . NUMERICAL EXAMPLES / 3.4 48 Ul 0.9 -0.8 : 0.7 : 0.6: O < 0.5 -} Q_ 5 0-4-j z £ 0.3 Ul Q 2 0.2-1 0.1 0 EXPERIMENTAL PERMANENT DISPLACEMENT '7 I Legend • EXPERIMENT [39] EL-SH. THIS STUDY EL-SH-SR. THIS STUDY 500 1000 1500 2000 TIME (//sec.) ' i ' 2500 3000 F I G U R E 3 . 9 Comparision of predicted and experimental responses of an explosively loaded clamped beam. 0.040 < CH r -CO Legend EL-SH. THIS STUDY EL-SH-SR. THIS STUDY 1 I 1 I 1 1 1 I 1 1 —1 ' 1 ' I • 200 400 600 800 1000 1200 1400 1600 1800 2000 TIME (fisec.) F I G U R E 3 . 1 0 Strain history of an explosively loaded clamped beam. NUMERICAL EXAMPLES / 3.4 49 1 2 3 LOCATION FROM MIDSPAN F I G U R E 3.11 Deflection profiles of an explosively loaded clamped beam. Figure 3.11 shows the deflection profile of the half span of the EL-SH-SR beam at various times throughout its reponse. Moderately large rotations of the beam axis are evident in these deflection profiles. In particular, the deflection profile at 200 fisec indicates a maximum beam axis rotation of the order of 10°. Regarding the efficiency of this numerical solution procedure it is to be noted that in order to generate the non-linear transient response history of the EL-SH-SR beam for a duration of 2500usee, based on a time step size of 0.90A£ c r , only 15.4 seconds of CPU time on an Amdahl 470 V/8 computer was required. This demonstrates that this numerical method provides a relatively inexpensive means, in terms of computer time, of simulating the non-linear transient response of beams. The accuracy of these numerically generated response histories as compared with the experimental results will now be discussed. Figure 3.9 shows the response of the test beam as obtained by experiment. Clearly, the EL-SH modelling of the beam is seen to over estimate the experimentally obtained peak deflection as well as the permanent NUMERICAL EXAMPLES / 3.4 50 midspan deflection. The EL-SH-SR modelling of the beam, on the other hand, offers reasonably good results with respect to the peak deflection and the permanent midspan deflection. There remains, however, some unexplained differences between the numerically predicted and experimentally observed time histories. This may be partly attributable to an over simplification of the strain-rate behaviour of the material as accounted for by the Cowper-Symonds relationship. Also, the explosive load induced some dynamic twisting into the response of each of the test beams. This surely affects the possible correlation which can be made between the numerically predicted response history and the experimentally obtained results. This particular beam problem has also been analyzed numerically by Witmer et ai.[40], using a spatial and temporal finite difference solution procedure. The half span was modelled with 30 finite difference stations. Both E L - S H and EL-SH-SR material modelling were considered in their study. The response histories obtained by Witmer et al., for the two material models, as given in Reference 40, are essentially identical with the respective response curves obtained in this study, and shown in Figure 3.9 . (b) Step Loaded I-Beam As a final example an axially constrained pin-ended steel I-beam (10 WF 45) subjected to a uniformity distributed step load of finite duration is considered.! Besides manifesting geometric and material non-linearities in the response this example tests the ability of this solution procedure to model the response of a beam v/ith a non rectangular cross-section. The pertinent material and geometric properties of the beam, as well as the imposed loading history, are given in Figure 3.12 . Note, that for this example the duration of the step load is approximately equal to one-tenth the fundamental period of the beam and the load magnitude is 5 times the static collapse load. Symmetry of the beam configuration and the applied load allows for only one half of the beam span to be analyzed. The finite element discretization consists of a uniform mesh of 10 elements. The \ In this example the end sections of the I-beam are constrained against both axial and transverse displacements. These end sections, however, are free to rotate about the centroidal axis of the I-beam, as shown in Figure 3.12 . NUMERICAL EXAMPLES / 3.4 51 T~ 10 elements in half span - <t 1 1 t L/2 — • 4 L/2 • L = 400.0 in 6 = 8.02 in (width) d= 10.12 in = 0.618 in 5P„ 10Wf45 T Step Load i„ = 0.350 in £ = 3.0 x 107 psi aa = 3.6 x 10* psi p = 0.733 x IO - 3 lb sec2/in* P„ = 97.2 tt/tn r = 0.01 (, see F I G U R E 3.12 Axially constrained pin-ended I-beam subjected to a uniformily dis-tributed step load of finite duration. 15-14-13-/—s • • c 12-• 11-INT 10-U J 9-U l • o 8-< —J 7-to • O 6-z 5-< Q_ CO 4-O • 3-2-1-0-/ \ PERMANENT MIDSPAN DISPLACEMENT PREDICTED BY THIS STUDY. \ V / \ Legend RICID PLASTIC THEORY - SCHUBAK W\ ELASTIC-PLASTIC THEORY - THIS STUDY F I G U R E 3.13 20 30 40 50 60 70 80 90 100 110 120 130 U0 TIME (millisec.) Dynamic large displacement elastic-plastic response of an axially con-strained pin-ended I-beam. NUMERICAL EXAMPLES / 3.4 52 critical time step size, on the central difference integrator, for this finite element mesh is 65.9 fisec. The critical time step can be expressed alternatively as a fraction of the fundamental period as Iy/1760. Figure 3.13 illustrates the non-linear response history for the I-beam, as pre-dicted by this study. The peak deflection occuring during the first cycle of the response is of the order of 1.5 beam depths, indicating the presence of non-linear geometric ef-fects. Note, that this peak deflection occurs in the free vibration phase for the beam, well after the removal of the step load. Non-linear stiffening of the structure is evident from the observation that the beam oscillates with a period approximately one half that of the fundamental period (where Tj = 0.116 sec). Also, in Figure 3.13, comparison is made between the numerically generated response and an analytical solution obtained by Schubak [41]. This analytical approach assumed rigid plastic material behaviour, with a simplified linear interaction between the axial force and moment. However, the full effects resulting from non-linear geometry were included in the analysis. While the two solutions in Figure 3.13 differ significantly with respect to the overall response history, which is to be expected, fairly good agreement is achieved in regards to the estimated permanent midspan deflection. A series of deflection profiles, illustrating the modes of deformation that the beam assumes in response to the imposed loading condition, are shown in Figure 3.14 . Examination of these deflection profiles attests to the presence of moderately large beam rotations during the response. In particular, the deflection profile at 0.03 sec indicates a maximum beam axis rotation of the order of 4°. To appreciate more fully the extent of material non-linearity manifested during the response the strain history of a particular point on the beam, located at (i,z) = (190.0in, 5.08in), is given in Figure 3.15 . This monitoring point is located on the top surface of the upper flange at the mid-point of the 10th element, which is ad-jacent to the midspan of the beam. From the stress-strain relationship it is determined NUMERICAL EXAMPLES / 3.4 53 LOCATION FROM MIDSPAN (in.) F I G U R E 3.14 Deflection profiles of an axially constrained pin-ended I-beam. 0.014 0.012 0.010 0.008 0.006 0.004 0.002 0.000 -0.002 - / / -/ / / YIELD STRAIN LEVEL S 1 1 I 1 1 1 1 1 1 1 10 15 20 25 30 35 TIME (millisec.) 40 45 50 55 F I G U R E 3.15 Strain history of an axially constrained pin-ended I-beam. NUMERICAL EXAMPLES / 3.4 54 TIME (millisec.) F I G U R E 3.16 Interaction between axial force and moment over time for an axially constrained pin-ended I-beam. that incipient plastic straining of the material occurs at eQ = 0.0012. Figure 3.15 confirms that at this monitoring point along the beam a significant amount of plastic strain occurs during the first cycle of the response, with the peak strain reaching a magnitude of 12e0. An alternate approach to analyzing the extent of material non-linearity present during the response period is possible from examination of the axial force and moment at a given cross-section of the beam. Figure 3.16 gives the response histories for the resultant axial force, M/Mc, and moment, M/Mo) expressed in non-dimensional form, for the beam cross-section, located at i = 190in . As given above, M„ and M0 represent, respectively, the fully plastic axial load and moment capacity of beam. Yielding of the outer fibres of this I-beam are initiated when the following inequality is satisfied; 1.12\M/M0\ + {M/M0)>1 (3.4.1) A determination of the extent of plastification through the cross-section, over the response period, can be obtained from considering the interaction relationship, from rigid plastic the-NUMERICAL EXAMPLES / 3.4 55 ory, between axial force and moment. According to this theory a fully plastic cross-section results when the interaction relationship between axial force and moment, expressible as /(J\//JV0, M/M0) = 1, is satisfied. The explicit form of this interaction relationship for an I-beam is as follows [42]; with c = A 2/(4rwZ), el = 4bZ/A2, c2 = A/(2bd- A), Z = M0/a0, &AW = {d-2tf)t. In point of fact, it is not possible for a fully plastic moment to develop in the beam since an elastic region will always subsist in the cross-section. Therefore, in actuality Equation (3.4.2) cannot be satisfied at any given cross-section where bending moments are present. Nonetheless, Equations (3.4.1) & (3.4.2) can be used to provide an indication of the extent of plastification resulting in the I-beam. First, from examination of Figure 3.16, in the light of Equation (3.4.1), it is observed that yielding of the cross-section is initiated at approximately 0.010 sec into the response. Again with reference to Figure 3.16, the greatest plastification of this cross-section is seen to occur at approximately 0.017 sec According to Equation (3.4.2) for all intents and purposes the cross-section at this time can be assumed to have fully yielded as a result of the imposed loading condition. non-linear transient response histories of beams of non-rectangular cross-section in the literature, to the author's knowledge, does not permit more to be said regarding the accuracy of the present numerical method in this area of application. otherwise. (3.4.2) Finally, it is to be noted that the lack of numerical and/or experimental CONCLUSION / 4.1 56 CHAPTER 4 CONCLUSION 4.1 SUMMARY A numerical solution procedure for predicting the transient response of slen-der ductile beams, exhibiting geometric and/or material non-linear behaviour, has been presented in this study. The aim of this presentation has been threefold; first, to formulate this beam problem in a manner consistent with the underlying principles of continuum mechanics; secondly, to bring together the research work of others, as it applies to this problem; finally, in the light of the first two objectives, to present a numerical solution algorithm which is simple and straightforward yet yields accurate results. In the theoretical development of the problem a variational form of the non-linear governing equation of motion for a Benoulli-Euler beam, based on a Total Lagrangian Approach, was presented. In this formulation beam deformations were restricted to the simple case of planar bending and stretching. Furthermore, beam kinematics were limited to the range of small strains and moderately large displacements. Due to the assumed slenderness of the beams considered both shear deformation and rotatory inertia effects were neglected. CONCLUSION / 4.1 57 The numerical solution procedure was initiated by application of the assumed displacement version of the Finite Element Method, utilizing simple one dimensional beam elements, from which the non-linear equations of motion for the spatially discretized beam were formed. Non-linear constitutive behaviour was conveniently accounted for by appli-cation of the "mechanical sublayer model*. This method of material modelling has the ability of giving an accurate piece-wise linear approximation to the actual stress-strain behaviour of the material under cyclic loading conditions. Kinematic strain-hardening, as well as the Bauschinger effect, are automatically accounted for with this material modelling procedure. Visco-plastic effects have been acknowledged in a rather simplistic fashion by employing the Cowper-Symonds strain-rate relationship. With an objective of furnishing a simple solution algorithm the central dif-ference method was employed to integrate the equations of motion. The accuracy and efficiency of this integration scheme was found to be conveniently enhanced by employ-ing a diagonal mass matrix in the solution procedure. The major shortcoming with the central difference method is that it is conditionally stable, requiring a rather restrictive time step to ensure stability of the integration scheme. A simple formula for determining a conservative bound on the limiting time step size, for linear beam problems, was also presented. It was found that by reducing this time step 10% — 20% a stable solution to the corresponding non-linear beam problem was obtained in all the examples tested. Stability of the solution algorithm was further ensured by introducing an energy balance checking procedure. To access the accuracy and efficiency of the solution algorithm a variety of examples, which exhibited geometric and/or material non-linear behaviour, were gener-ated. Results from all the examples were compared with solutions from other sources. In general, good agreement was obtained in all cases. Where discrepancies were found in the response predictions they were identifiable with particular simplifications made during the formulation of the solution procedure. In one of the examples, it was found that the peak CONCLUSION / 4.2 58 midspan deflection, as predicted by this study, for a simply supported elastic perfectly-plastic beam, which had a span to depth ratio of 15:1, was approximately 6% less than that given by a full plane stress finite element solution. From this example, it is seen that further study is required to determine the range of applicability of this solution strategy in accurately determining the non-linear transient response of beams with various span to depth ratios. In another example, it was noted that the Cowper-Symonds strain-rate relationship, while improving the predicted response history, did not capture fully the visco-plastic behaviour of the beam. Further study is required to access the performance of this approximate constitutive law in relation to a fully visco-plastic analysis. As previously stated, the motivation behind this study has been to develop an efficient and accurate solution algorithm which could be employed in various parametric studies relating to the the non-linear transient response of slender ductile beams to air-blast pressure pulses. It is believed that this has been achieved. 4.2 F U T U R E A R E A S OF R E S E A C H There are several obvious desirable extensions of the present formulation which can be cited and which are currently being studied. The first considers the inclusion of shear deformation and rotatory inertia effects, through Timoshenko beam theory, into the solution procedure. Afterwhich, it is planned to extend the solution strategy so that the non-linear transient response of laminated fibre reinforced composite beams can be numerically predicted. 59 REFERENCES [1] Symonds, P.S. and Mentel, T.J. "Impulsive loading of Plastic Beams with Axial Constraints," Journal of the Mechanics and Physics of Solids, Vol. 6, 1958, pp. 186-202. [2] Witmer, E.A., Balmer, H.A., Leech, J.W., and Pian, H.H. "Large Dynamic Deformations of Beams, Rings, Plates and Shells," AIAA Journal, Vol. 1, No. 6, 1963, pp. 1848-1857. [3] Salus, W.P., Ip, C , and Vanderlinden, J.W. "Design Considerations of Elas-tic Plastic Structures Subjected to Dynamic Loads," AIAA/ASME 11th Structures, Structural Dynamics, and Material Conference, Denver, Col., April 1970. [4] Wu, R.W.H., and Witmer, E.A. "Finite-Element Analysis of Large Elastic-Plastic Transient Deformations of Simple Structures," AIAA Journal, Vol. 9, No. 9, 1971, pp. 1719-1724. [5] Wu, R.W., and Witmer, E.A. "Nonlinear Transient Responses of Structures by the Spatial Finite-Element Method," AIAA Journal, Vol. 11, No. 8, 1973, pp. 1110-1117. [6] Yang, T.Y, and Siagal, S. "A Simple Element for the Static and Dynamic Re-sponse of Beams with Material and Geometric Nonlinearites," International Journal for Numerical Methods in Engineering, Vol. 20, 1984, pp. 851-867. [7] Rodal, J.J.A., Steigmann, D.J., and Witmer, E.A. "Numerical Simulation of Transient Finite Deformations of Thin Beams," ASME Journal of Applied Mechanics, Vol. 50, 1983, pp. 765-769. [8] Stagliano, T.R., and Mente, L.J. "Large Deflection, Elastic-Plastic Dynamic Strutural Response of Beams and Stiffened or Unstiffened Panels- a Com-parision of Finite Element, Finite Difference and Modal Analysis, Nonlinear Finite Element Analysis and ADINA," Proceedings of ADINA Conference, 1979. [9] Bathe, K.J., Ozdemir, H., and Wilson, E.L. "Static and Dynamic Geometric and Material Nonlinear Analysis," Report No. UCSEM 74-4, 1974, Struc-tural Engineering Laboratory, University of California (Berkeley). REFERENCES 60 10] McNamara, J.F. "Solution Schemes for Problems of Nonlinear Structural Dynamics," ASME Journal of Pressure Vessel Technology, Vol. 96, 1974, pp. 96-102 11] Mikkola, M.J., Tuomala, M. , and Sinisalo, H. "Comparision of Integration Methods in the Analysis of Impulsively Loaded Elasto-Plastic and Viscoplas-tic Structures," Computers <fe Structures, Vol.14, No. 5-6, 1981, pp. 469-478. 12] Jones, N. "A Literature Review of the Dynamic Plastic Response of Struc-tures," Shock and Vibration Digest, Vol. 7, August 1975, pp.89-105. 13] Jones, N. "Recent Progress in the Dynamic Plastic Behaviour of Structures,", Part 1 and 2 , Shock and Vibration Digest, Vol. 10, September 1978, pp. 21-33 and October 1978, pp. 13-19. 14] Jones, N. "Recent Progress in the Dynamic Plastic Behaviour of Structures,", Part 3, Shock and Vibration Digest, Vol. 13, October 1981, pp. 3-16. [15] Jones, N. "Recent Progress in the Dynamic Plastic Behaviour of Structures,", Part 4, Shock and Vibration Digest, Vol. 17, February 1985, pp. 35-47. [16] Fung, Y.C. , Foundations of Solid Mechanics, Prentice-Hall, 1965. [17] Washizu, K., Variafcionai Methods in Elasticity and Plasticity, 3rd Edition, Pergamon Press, 1982. [18] Ibid, pp. 207-210. [19] Bathe, K.J. , Finite Element Procedures in Engineering Analysis, Prentice-Hall, 1982, pp. 318-334. [20] Zienkiewicz, O.C., The Finite Element Method, 3rd Edition, McGraw-Hill, 1977. [21] Cook, R.D., Concepts and Applications of Finite Element Analysis, 2nd Edition, John Wiley k Sons, 1981, pp. 68-79. [22] Ibid, pp. 27-31. [23] Zienkiewicz, O.C., Nayak, G.C., and Owen, D.R.J. "Composite and 'Over-lay' Models in Numerical Analysis of Elastic-Plastic Continua."Internationa] Symposium on Foundations in Plasticity,ed. Sawczuk, A. Warsaw, 1972. [24] Wu, R.W.H., and Witmer, E.A. "Finite-Element Analysis of Large Elastic-Plastic Transient Deformations of Simple Structures," AIAA Journal, Vol. 9, No. 9, 1971, pp. 1719-1724. [25] Bodner, S.R., and Symonds, P.S. "Experimental and Theoretical Investiga-tion of the Plastic Deformation of Cantilever Beams Subjected to Impulsive Loading," ASME Journal of Applied Mechanics, Vol. 29, 1962, pp. 719-728. REFERENCES 61 [26] Wu, R.W.H., and Witmer, E.A. "Finite-Element Analysis of Large Elastic-Plastic Transient Deformations of Simple Structures," AIAA Journal, Vol. 9, No. 9, 1971, pp. 1719-1724. [27] Ibib. [28] Stagliano, T.R., and Mente, L.J. "Large Deflection, Elastic-Plastic Dynamic Strutural Response of Beams and Stiffened or Unstiffened Panels- a Com-parision of Finite Element, Finite Difference and Modal Analysis, Nonlinear Finite Element Analysis and ADINA," Proceedings of ADINA Conference, 1979. [29] Shantram, D., Owen, D.R.J., and Zienkiewicz, O.C., "Dynamic Transient Behaviour of Two- and Three-dimensional Structures Including Plasticity, Large Deformation Effects and Fluid Interaction," Earthquake Engineering and Structural Dynamics, Vol. 4, 1976, pp. 561-578. [30] Owen, D.R.J., and Hinton, E., Finite Elements in Plasticity: Theory and Practise, Pineridge Press, 1980, pp. 431-463. [31] Krieg, R.D., and Key, S.W. "Transient Shell Response by Numerical Time Integration," Internationa] Journal for Numerical Methods in Engineering, Vol. 7, 1973, pp. 273-286. [32] Ibid. [33] Irons, B., and Ahmad, S., Techniques of Finite Elements, Elis Horwood Ltd., 1980, pp. 428-429. [34] Belytschko, T. "Explicit Time Integration of Structure-Mechanical Systems," Advanced Structural Dynamics, ed. Donea, J. , 1978. pp. 97-122. [35] Bathe, K.J. , Ozdemir, H., and Wilson, E.L. "Static and Dynamic Geometric and Material Nonlinear Analysis," Report No. UCESM 74-4, 1974, Struc-tural Engineering Laboratory, University of California (Berkeley). [36] Liu, S.C., and Lin, T.H. "Elastic Plastic Dynamic Analysis of Structures Using Known Elastic Solutions," Earthquake Engineering and Structural Dynamics, Vol. 7, 1979, pp. 147-159. [37] Mondkar, D.P., and Powell, G.H. "Finite Element Analysis of Non-Linear Static and Dynamic Response," International Journal for Numerical Meth-ods in Engineering, Vol. 11, 1977, pp. 499-520. [38] Yang, T.Y, and Siagal, S. "A Simple Element for the Static and Dynamic Re-sponse of Beams with Material and Geometric Nonlinearites," Internationa] Journal for Numerical Methods in Engineering, Vol. 20, 1984, pp. 851-867. REFERENCES 62 [39] Witmer, E.A., Balmer, H.A., Leech, J.W., and Pian, H.H. "Large Dynamic Deformations of Beams, Rings, Plates and Shells," AIAA Journal, Vol. 1, No. 6, 1963, pp. 1848-1857. [40] Ibid. [41] Schubak, R. "Finite Deflection Dynamic Response of Axially Restrained Beams," M.A.Sc. Thesis, U.B.C., Vancouver, British Columbia, (in progress). [42] Massonnent, C.E., and Save, M.A., Plastic Analysis and Design Volume 1: Beams and Frames, Blaisdel Publishing, 1965, pp. 128-135. [43] Leech, J.W., Hsu, P., and Mack, E.W. "Stability of a Finite Difference Method for Solving Matrix Equations," AIAA Journal, Vol. 3, No. 11, 1965, pp. 2172-2173. STABILITY CRITERION FOR THE CENTRAL DLFFERENCE METHOD / A.O 63 APPENDIX A STABILITY CRITERION FOR THE CENTRAL DIFFERENCE METHOD In section 2.6, it was noted that the central difference time operator is con-ditionally stable. When applied to problems of structural dynamics the critical time step, Atcr, is given by the following relationship; ^cr = 2 / w m a „ {A.l) where 0Jt7iax represents the largest natural frequency occuring in the discretized system. A proof of this fact, for a linear structural system, is given here for completeness. This proof follows closely the procedure presented by Leech et ai.[43]. The equations of motion for the discretized structure are expressible as M T J + K ' D = ' F , (A.2) where M and X are, respectively, the global mass and stiffness matrices. F is the global forcing function. Applying the central difference method to the homogeneous form of the linear equations of motion, at time t — nAt in the analysis, yields The forcing function has not been included in the above equation since it plays no role in establishing the stability criterion. A solution is sought of the form 'D = &ePnAt, where /? STABILITY CRITERION FOR THE CENTRAL DIFFERENCE METHOD / A.O 64 is an undetermined parameter and $ is an arbitrary shape vector. Substituting this trial solution into Equation (A.3) yields 2 ( e ^ ) + ( A t ' M - ' K - 2) + 1 $ = 0. {AA) For each eigenvalue-eigenvector pair in Equation (AA) their exists a separate /?. Let a given eigenvalue of M - 1 K be denoted by w2. Associated with this eigenvalue the equation for (3 can be written as f 2 - 2AZ + 1 = 0, (A.5) where, £ = e^1 and A = 1 - ^ A t V . The solution to Equation (A.5) is given by the two roots ^ A + U 2 - ! ) 1 / 2 , 6 = A - U 2 - i ) 1 / 2 . (A.6) (A.7) The stability criterion is derived from the requirement that the response must remain bounded as t —> oo, so that \e^l\ = M\<l. (A.8) This in turn, implies that — 1 < A < 1, which can be expressed alternatively as -2 < - - A t V < 0. - 2 This imposes the condition that At < 2/u2. Since the largest natural frequency, ujmax imposes the strictist requirement on At the stability criterion is given by At < 2/umax.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical simulation of the non-linear transient response...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical simulation of the non-linear transient response of slender beams Folz, Bryan 1986
pdf
Page Metadata
Item Metadata
Title | Numerical simulation of the non-linear transient response of slender beams |
Creator |
Folz, Bryan |
Publisher | University of British Columbia |
Date Issued | 1986 |
Description | A simple numerical solution strategy for predicting the transient response of slender ductile beams, exhibiting geometric and/or material non-linear behaviour, is presented in this study. In the theoretical development of the problem the governing non-linear equation of motion, in variational form, for the bending and stretching of a Bernoulli-Euler beam is established. The numerical solution procedure is then initiated by employing the assumed displacement version of the Finite Element Method with 1-dimensional 6-DOF beam elements. Elastic-plastic strain-hardening of the beam material is conveniently accounted for by means of the "mechanical sublayer model". Visco-plastic material behaviour is included in the analysis through a simple strain-rate dependent constitutive relationship. The equations of motion for the spatially discretized beam are integrated time-wise by means of the central difference method. A variety of examples are then solved and the results compared with solutions from other sources. In general, the numerical solution strategy yields an efficient and accurate modelling of the non-linear transient response of slender ductile beams. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-07-10 |
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. |
DOI | 10.14288/1.0062580 |
URI | http://hdl.handle.net/2429/26287 |
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 F66.pdf [ 4.1MB ]
- Metadata
- JSON: 831-1.0062580.json
- JSON-LD: 831-1.0062580-ld.json
- RDF/XML (Pretty): 831-1.0062580-rdf.xml
- RDF/JSON: 831-1.0062580-rdf.json
- Turtle: 831-1.0062580-turtle.txt
- N-Triples: 831-1.0062580-rdf-ntriples.txt
- Original Record: 831-1.0062580-source.json
- Full Text
- 831-1.0062580-fulltext.txt
- Citation
- 831-1.0062580.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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062580/manifest