

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


DEVELOPMENT OF A FINITE ELEMENT ANALYSIS PROGRAM FOR STRUCTURAL DYNAMICS APPLICATIONS By MUHAMMET ALI SAGLAR Bachelor of Civil Engineering Istanbul Technical University Istanbul, Turkey 2007 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE July, 2009 ii DEVELOPMENT OF A FINITE ELEMENT ANALYSIS PROGRAM FOR STRUCTURAL DYNAMICS APPLICATIONS Thesis Approved: Dr. Jonathan S. Goode Thesis Adviser Dr. Robert N. Emerson Dr. Tyler Ley Dr. A. Gordon Emslie Dean of the Graduate College iii ACKNOWLEDGMENTS I thank God for His grace, His blessings and guiding me to the right path. I would like to express profound gratitude to my advisor, Dr. Jonathan S. Goode, for his precious support, encouragement and useful suggestions throughout this research work. His moral support and continuous guidance enabled me to complete my work successfully. I am also highly thankful to Dr. Robert Emerson and Dr. Tyler Ley for their encouragements guidance rendered to me throughout my research. Their deep experience helped me a lot to be a structural engineer. I would also like to thank to Rameez Iqbal, M.S student, for his support during this research. I am really grateful to him for his friendship and understanding. My deepest appreciation goes to my parents, Remzi Saglar and Asiye Saglar for their invaluable supports during every level of my life. I am grateful also to my three brothers their unflinching courage and conviction. iv TABLE OF CONTENTS Chapter Page I. INTRODUCTION ......................................................................................................1 1.1 Background and Problem Definition .................................................................1 1.2 Objectives ..........................................................................................................2 1.3 Overview ............................................................................................................3 II. LITERATURE REVIEW ..........................................................................................6 2.1 Finite Element Method ......................................................................................6 2.2 Structural Dynamics...........................................................................................8 2.3 MATLAB Usage in Structural Dynamics........................................................10 III. FINITE ELEMENT THEORETICAL DERIVATIONS .......................................12 3.1 Introduction ......................................................................................................12 3.2 Interpolation and Shape Functions...................................................................13 3.3 Formulation of Element Stiffness Matrices .....................................................16 3.4 Element Boundary and Loading Conditions ....................................................17 3.5 Solution Techniques.........................................................................................19 3.6 Structural Dynamics.........................................................................................21 3.6.1 Equation of Motion .................................................................................21 3.6.2 Formulation of Element Mass Matrices ..................................................22 3.6.3 Dynamic Loading Conditions .................................................................23 3.6.4 Numerical Integration .............................................................................26 3.6.5 NewmarkBeta Solution Technique ........................................................27 IV. CASE STUDIES ....................................................................................................29 4.1 Introduction ......................................................................................................29 4.2 Harmonic Loading ...........................................................................................31 4.3 Impulse Loading ..............................................................................................36 4.4 Step Loading ....................................................................................................42 4.5 Multiple Dynamic Loads .................................................................................47 v Chapter Page V. CONCLUSION ......................................................................................................54 5.1 Summary and Conclusions ..............................................................................54 5.2 Recommendations ............................................................................................55 REFERENCES ............................................................................................................57 APPENDIX A – MATLAB FEA PROGRAM ............................................................59 APPENDIX B – CONFIRMATION CALCULATIONS ............................................68 vi LIST OF TABLES Table Page 41 Case Studies .....................................................................................................30 vii LIST OF FIGURES Figure Page 21 Crack Distributions at the Ultimate Load Step (Ozcan, 2008) ..........................8 31 Q4 Finite Element Nodal Degrees of Freedom ................................................15 32 Harmonic Loading ...........................................................................................24 33 Impulse Loading ..............................................................................................25 34 Step Loading ....................................................................................................25 41 Case Study 1 – Cantilever Beam .....................................................................31 42 Case Study 1 – Harmonic Loading ..................................................................32 43 Case Study 1 – Finite Element Mesh ...............................................................33 44 Case Study 1 – Displacements .........................................................................33 45 Case Study 2 – SimplySupported Beam .........................................................37 46 Case Study 2 – Impulse Loading .....................................................................37 47 Case Study 2 – Finite Element Mesh ...............................................................38 48 Case Study 2 – Displacements .........................................................................39 49 Case Study 3 – Continuous Beam ....................................................................42 410 Case Study 3 – Step Loading ...........................................................................43 411 Case Study 3 – Finite Element Mesh ...............................................................44 412 Case Study 3 – Displacements .........................................................................44 413 Case Study 4 – Continuous Beam ....................................................................48 414 Case Study 4 – Multiple Loading ....................................................................49 415 Case Study 4 – Finite Element Mesh ........................................................................... 50 416 Case Study 4 – Displacements .........................................................................50 1 CHAPTER I INTRODUCTION 1.1 Background and Problem Definition Finite element analysis (FEA) is a computational technique used to find approximate solutions of field problems that are described by differential equations or by an integral expression. The finite element solution technique helps in solving complex field problems by numerical approximations of the differential or integral expression. Finite element analysis has significant advantages when compared to the other numerical analysis methodologies. It is very powerful and applicable in many engineering problems such as displacements of a structural systems, stressstrain analysis, heat transfer and magnetic fields. In addition, for individual elements different material properties can be incorporated. One of the most important advantages of finite element analysis is that there are no limitations concerning the geometry or boundary conditions. Different types of geometry and boundary conditions can be accommodated easily. In addition to these properties, it is easy to modify the problem and increase the accuracy of the results while usually only at the expense of computational efficiency. 2 At present many commercial software packages use the finite element analysis because of the advantages mentioned above. ANSYS, ALGOR, ABAQUS, COSMOS/M and SAP2000 are some of the wellknown commercial software packages that are based on finite element analysis. That is to say that the finite element analysis has wide range of usage and is conveniently available for engineers and researchers. For research purposes, finite element analysis is very useful to allow for examination of “what if” design scenarios. As mentioned earlier, it has almost no limitation for engineering problems in terms of geometry and boundary conditions. Because of these properties it is very powerful for researchers. In this study, four types of dynamic loads are applied separately to different beams. The finite element analysis is used to define the beam and determine the problem solution for both loading cases. Furthermore, Q4 finite elements are considered during the analysis and these elements are combined to form a mesh. During the dynamic analysis calculation procedure a numerical procedure, Newmark’s method, is used to numerically integrate the equations of motion with respect to time. 1.2 Objectives The primary objective of this research is the development of a finite element analysis computer program written using MATLAB. This program will be used to determine structural displacements or response of structural systems to both static and dynamic 3 loading conditions. As a method of illustrating the use of this program, several case studies are presented that find the displacements and deformed shape of a beam under the dynamic loading with change of time by using finite element method. The purpose of this study is development of a coupled finite element analysis/structural dynamics computer program. This study does not focus on accuracy of the results. In other words, the best approximation to the exact results is not the main intention. As previously mentioned, the computer codes are written in MATLAB which is a computationally powerful programming language for research purposes. Explained in its most basic form, this computer code has the capability to calculate approximate displacements which can be calculated using design software packages. Besides these capabilities, the development of this program gives the author the understanding of what is going on behind the software packages. This program has the capacity to incorporate different scenarios for dynamic loading. In other words, any timedependent loading scenario to calculate the dynamic load can be applied to the beam to obtain displacements and the deformed shape of the structure. 1.3 Overview The next several chapters present the methodology behind this study. A summary of each chapter is provided to give a brief overview of the remaining sections of this study: • Chapter 2 – Literature Review: 4 Current and recent studies are provided from a broad viewpoint including a brief discussion about the usage of the finite element method and structural dynamics. The literature review gives several examples that use the finite element method, structural dynamics and structural dynamics combined with the finite element method. • Chapter 3 – Finite Element Theoretical Derivations: Detailed theoretical development of the finite element analysis procedure and structural dynamics is presented and discussed. This information is directly used to develop the MATLAB code for the FEA program. Theoretical equations are derived and explained in this chapter. Topics that are covered include discussions about generating element matrices, boundary and loading conditions, finite element solution techniques, and structural dynamics solution technique. • Chapter 4 – Case Studies: The capabilities of the MATLAB FEA program developed as a consequence of this study are demonstrated. Four case studies considering a dynamic loading are demonstrated on different structural beam configurations. Hand calculations (included in Appendix B of this study) are provided for confirmation of the FEA program. 5 • Chapter 5 – Summary, Conclusion and Recommendations: A summary of the results and conclusions of the study are provided. This includes a brief discussion and capability of this study. In addition, recommendations about the study and its expandability are provided. 6 CHAPTER II LITERATURE REVIEW 2.1 Finite Element Method The finite element method is an approximation method that can be used to calculate stresses, movements of loads and forces, displacements, heat transfer and other basic physical behaviors while using very large matrix arrays and mesh diagrams. In recent years, the finite element method has been used to obtain the solutions to a variety of engineering problems. A few applications of finite element modeling that have been used for engineering research are subsequently illustrated. Li et al. (2001) presented a quadratic finite element considering the principle that the local displacement fields of the elements are compatible with the global displacement field of the corresponding systems by using generalized degrees of freedom (GDOF) and a quadratic finite strip with GDOF. They developed a global displacement field based on quadratic Bspline functions and created local displacement field strips for the elements. Therefore, they formed models for the finite elements and finite strips corresponding to the GDOF by rearranging the local displacement fields of elements and strips to be compatible with their corresponding global displacement fields. In their study, several 7 numerical examples were provided for authenticity, simplicity and versatility of the element and strip in the analysis of thinwalled structures. The authors determined that the finite element and finite strip with GDOF can handle inconveniences in the analysis of beams, plates and other structures when the thickness changes. In addition to handling such problems, the authors also determined that the finite element and finite strip gave similar results with fewer degrees of freedom when compared to traditional finite element solutions. Özcan et al. (2008) presented three steel fiberadded reinforced concrete models which were analyzed experimentally and numerically with a finite element analysis. They created a finite element model using ANSYS. The authors considered eightnode solid brick elements in their finite element modeling. Experimental data and finite element analysis results were compared. Figure 21 is given as an illustration for comparison of experimental research and finite element modeling. As a result of their study, the experimental data and finite element modeling gave very close results in terms of deflections and stresses at the center line along with initial and progressive cracking, failure mechanism and load, deflections and stresses at the zero deflection point and decompression. 8 Figure 21: Crack Distributions at the Ultimate Load Step (Ozcan, 2008) 2.2 Structural Dynamics Incorporation of the dynamic response of structural systems into a finite element model is a natural extension of the numerical timestepping methods for dynamic analysis. Recently, many studies have considered this kind of “coupling” of a FEA program and structural dynamics models. A review of some of these studies is provided herein. 9 Du et al. (1992) considered a threedimensional elastic beam with an arbitrary and large moving base with six degrees of freedom by creating a general finite element structural dynamics model. In their study, six degrees of freedom of the beam base may include either a particular arbitrary motion of the base or the coupling of the beam with other structures. The beam can be pretwisted and has a mass center offset from the elasticity center. In the study, the authors derived the equations of motion using the principal of virtual work in the finite element analysis. To simplify the analysis, the authors combined the beam inertia at the end nodes of each element. At the conclusion of the study, it was determined that the model provided flexibility because of the combination of the finite element method and arbitrary base motion variables that were used in multibody dynamics and a fundamental element to solve the dynamic problems of rotating beamlike structures. Cerioni et al. (1995) presented a finite element model that was capable of analyzing the dynamic nonlinear behavior of unreinforced masonry panels in a biaxial stress field by using a nonconforming quadrilateral element. The study derived the equilibrium equations including the inertial and damping actions by a direct stepbystep integration procedure in the time domain known as Newmark’s method. The results were compared with experimental results in terms of displacement, velocity, acceleration, strain and stress. The comparison of the results indicated good performance not only for nonconforming quadrilateral elements but also timedependent and nonlinear analyses. In addition to these results, the model provided a very simple, computationally economical and convenient analysis of complex structural problems. 10 Mermertas et al. (1997) introduced a curved bridge deck to examine the dynamic interaction between a vehicle that is assumed to have fourdegrees of freedom and a curved beam. The finite element method was used to model the curved beam assumed to be simply supported neglecting any damping of the structure. The authors also applied a multipredictorcorrector procedure in conjunction with the Newmark method in the solution of the resulting equations of motion. The authors determined the midspan deflection of the bridge for different vehicle speeds and varying radii of curvature for the curved deck. 2.3 MATLAB Usage in Structural Dynamics The usage of computer programs to use the finite element method and numerical timestepping methods for the structural dynamics has become increasingly important in recent years. Advances in computer technology allow even simple computer systems to model very complex systems with ease and efficiency. A popular program in research studies is MATLAB. This program can be utilized very easily to accomplish both a finite element analysis and structural dynamics analysis due to its friendly matrix manipulations. Several recent research studies have considered MATLAB and they are further outlined here. Kiral et al. (2008) presented a fixedfixed laminated composite beam that was subjected to a concentrated force traveling at a constant velocity to determine the dynamic behavior 11 of a beam. The authors used classical lamination theory in order to create a threedimensional finite element model. In addition, the Newmark integration method to compute the dynamic response was implemented in MATLAB. The dynamic magnification, defined as the ratio between the dynamic and static displacements, was determined from their study. As a result of the study, it was determined that the load velocity and ply orientation may affect the dynamic response significantly. It was also determined that if the total traveling time is equal to the first natural period of the beam, the maximum deflection occurs at the midpoint of the beam. Bozdogan et al. (2009) demonstrated an approximate method using a continuum approach and transfer matrix method for static and dynamic analyses of multibay coupled shear walls. The authors assumed the structure system as a sandwich beam and thus wrote the differential equation whose solution gave the shape functions for each story of the sandwiched beam. The system modes and periods based on the boundary conditions and story transfer matrices, found by the shape functions, can then be calculated. Using MATLAB, a computer program was written to solve numerical examples to show the reliability of the method used. Results were compared with previous work done in the literature that was in good agreement. The authors suggested that their method is appropriate to use on a wide range of structural system applications. 12 CHAPTER III FINITE ELEMENT THEORETICAL DERIVATIONS 3.1 Introduction The primary objective of this study is to develop a finite element analysis program coupled with a structural dynamics response modeling program to determine the dynamic response of a structure due to a dynamic load. The use of the finite element method, however, influences the dynamic response of the structure. A variety of finite elements may be chosen for the analysis, each having its advantages and disadvantages. As such, one element is not always superior to another with respect to any given analysis. Often, it is the experience of the analyst that determines the appropriate finite element to be used in the analysis. The purpose of this study is not to determine the appropriateness of the finite element method chosen. In addition, there are also a variety of structural dynamic response algorithms that can be used to determine the response of a structure due to a dynamic load application. As such, it is also not the purpose of this study to determine the “best” algorithm to be used for determining the dynamic response of the structural system. Rather, this study seeks only to combine the theoretical development of the finite element method with an approach for determining the dynamic structural response into a seamless computer application that can be expanded for future use. 13 To develop the coupled finite element/structural dynamics program, an understanding of the theoretical development of both topics is necessary. Finite elements are the basis for modeling the structural system into a set of discretized pieces that can be assembled into a set of structural equations. The structural dynamics algorithm can then analyze these structural equations at discrete points in time due to a defined dynamic loading. The combination of these two procedures produces displacements of the structural system as a function of time. In the discussion to follow herein, the theoretical development will center on these topics. Included in these discussions will be formulations of element matrices, derivations of finite element method equations, dynamic analysis procedures and solution techniques which were used while developing the coupled finite element analysis/structural dynamics program. 3.2 Interpolation and Shape Functions The finite element method, at its most basic form, is a set of interpolations or approximations of a dependent variable with respect to an independent variable. Unknown degreesoffreedom (DOF) are utilized to ensure a singlevalued approximation for a set of conditions. An interpolating polynomial with dependent variable and independent variable x in terms of generalized DOF can be expressed in the form Σ or (31) where and can be written as vectors 14 1 … and … (32) For linear interpolation n can be taken 1 and for quadratic interpolation n can be taken 2. The can be written in terms of nodal values (known locations of x) of for known values of x. The relation between and the nodal values can be written as Φ (33) where each row of is calculated at the appropriate nodal location. From the equations above, substitution produces Φ where ! … (34) In matrix an individual is represented as a shape function. In this study, the bilinear rectangle (Q4) finite element is considered. Therefore, linear interpolation is considered when generating the shape functions. That is to say that between two points " , $ and " , $ on a linear line for 1 we can obtain % & ' ( where ) 1 1 * (35) Inverting and using equation 34, ! +,!+ . /1 1 0 and 1 +,!+ +,!+ +!+ +,!+ 2 (36) Equation 36 gives the shape functions of two points on a straight line which are called N1 and N2. 15 Consider the general Q4 element shown in Figure 31. To find the shape functions of the Q4 element, linear interpolation is made along the top and bottom sides to obtain side displacement u12 and u43. Thus, in equation 35, / and , so that 3 4!+ 4 3 5 46+ 4 , 378 4!+ 4 37 5 46+ 4 38 (37) Linear interpolation is then made in the y direction between u12 and u43 as 3 9!: 9 3 5 96: 9 378 (38) Substitution of equations 37 into 38 yields 3 Σ 3 which gives the shape functions of the rectangular Q4 element used in this study as 749 "1 / $"1 / ;$ 749 "1 5 $"1 / ;$ 8 749 "1 5 $"1 5 ;$ 7 749 "1 / $"1 5 ;$ (39) y x a a b b v1 v4 v3 v2 u4 u3 u1 u2 4 1 3 2 Figure 31: Q4 Finite Element Nodal Degrees of Freedom 16 A similar procedure is repeated for the DOF in the y direction. Namely, these DOF are v1, v2, v3 and v4. Similarly, the same shape functions given in equation 39 are found for the shape functions in the ydirection. Therefore, the same interpolation, or approximation, is made in both the x (or u) and y (or v) directions. 3.3 Formulation of Element Stiffness Matrices The principal of virtual work is used to obtain element stiffness matrix for the Q4 element shown in Figure 31. This is appropriate for commonly used elements, which are based on interpolation of displacements from nodal DOF. The principal of virtual work states < => ? @A < =3 B @A 5 < =3 Φ @C (310) where => , B and Φ represent the virtual strains produced by the virtual displacements, body forces and surface tractions, respectively. The displacements {u} are interpolated over an element utilizing shape functions such as those provided by equation 39 as 3 @ where 3 3 D E (311) > = 3 and > F @ where F G (312) where [B] is the straindisplacement matrix. From equations 311 and 312 =3 =@ and => =@ F (313) Substitution of equation 313 back into the statement of virtual work, equation 310, produces 17 =@ " < F H F @A / < F H I @A 5 < F ? @A / JB@A/ JΦ@C 0 (314) Equation 314 can be simplified to produce L @ M (315) As a result of the principle of virtual work, the element stiffness matrix can be determined as L < F H F @A (316) Specifically, for the Q4 element, equation 316 can be written as L < < F H F N@ @; 4 !4 9 !9 (317) For the case of twodimensional plane stress analysis (as considered in this study), the material constitutive matrix [E] is H O " !P,$ Q 1 D 0 D 1 0 0 0 " !P$ R (318) In twodimensional analyses, the thickness, t, in equation 317 is commonly taken as unity. 3.4 Element Boundary and Loading Conditions Before the solution of the structural equations, both boundary conditions and loading conditions must be prescribed for the system. Without boundary conditions the structural equations will not produce a single unique solution for the prescribed loading conditions. As such, the structural system will have rigid body motions. Without loading conditions 18 the structural equations will produce no displacements of the structural system. Thus, it is necessary that both boundary and loading conditions be prescribed for the structural system. These conditions are prescribed at particular DOF of the structural system. Boundary conditions, or support conditions, can be arranged by providing the appropriate stiffness to the related DOF to produce a prescribed displacement. For zero displacement, the prescribed stiffness can be a numerically large number (several orders of magnitude larger than the largest magnitude in the stiffness matrix) such that a relatively small displacement at that DOF is produced. The boundary conditions can be applied to any DOF on the structure no matter its direction. In the present study, only translational DOF are considered in the finite elements (i.e., there are no rotational DOF) and therefore only translational displacements are restricted with respect to particular support conditions such as a fixedsupport, pinnedsupport, or rollersupport. Loading conditions are prescribed in a fashion similar to boundary/support conditions. The load can be applied to the any DOF. The structural system may be subjected to a single loading condition or multiple loading conditions. As this study considers dynamic loads, beyond the typical static load cases considered in a typical finite element analysis, further explanation of the dynamic load is presented in Section 3.6.3. 19 3.5 Solution Techniques Equation 317 is integrated over a rectangular surface as given in Figure 31. Numerically, however, to integrate this equation is cumbersome and not efficient. Thus, integration is achieved using Gauss Quadrature. Gauss Quadrature is a numerical approximation of the integration by use of simple algebraic equations evaluated at specific points. To use Gauss Quadrature, however, the element must be formulated in the isoparametric space. The use of the isoparametric space specifies that an isoparametric element is used rather than the physical element. Thus, the integrands in the integration formulas are expressed as functions of S and T rather than x and y. For the function ø = ø (ξ, η), the Quadrature rule is given as U < < ø"ξ, η$dξ dη ! Σ ΣZY YZø"ξ, η$ ! (319) where Y and YZ are weighting factors for each Gauss point i and j. The weighting factors for twopoint Gauss Quadrature, as used in this study for the Q4 element, are taken as unity. For the Q4 element shown in Figure 31, the individual isoparametric shape functions can simply be determined from equation 39 by assigning a = 1, b = 1, x = ξ and y = η producing 7 "1 / S$"1 / T$ 7 "1 5 S$"1 / T$ 8 7 "1 5 S$"1 5 T$ 7 7 "1 / S$"1 5 T$ (320) 20 Since the isoparametric Q4 element is not in the same coordinate system as the physical Q4 element, a mapping function is needed to relate the two coordinate systems. The Jacobian matrix is used to accomplish this mapping. The Jacobian matrix is a scale factor that multiplies dξdη to produce the physical area increment dxdy and is expressed as (for the Q4 isoparametric element) [ 7 \ – "1 / T$ "1 / T$ "1 5 T$ – "1 5 T$ – "1 / S$ – "1 5 S$ "1 5 S$ "1 / S$ ^ _ ; ; 8 ;8 7 ;7 ` ) [ [ [ [ *(321) Thus, the element stiffness matrix for the Q4 isoparametric element is then L a F H F N @ @; < < F H F N [ @S @T ! ! (322) By substituting equation 322 into equation 319, the integration can be numerically calculated using Gauss Quadrature. Following the determination of the stiffness of each element in the finite element mesh, the system stiffness must be assembled. Based on the arrangement of the finite element mesh, the stiffness for each element corresponding to particular system DOF is “fed” into the system stiffness corresponding to the same system DOF. This process, called the assembly process, is a simply a mapping technique relating the DOF corresponding to each element to the DOF of the system. After assembly of the elements into the system equations, the structural equations are then in the form that can be solved. Gauss elimination is used to solve the structural equations for given boundary and loading conditions. In Gauss elimination, equations 21 [K]{D} = {R} are solved for {D}, the displacements, by reducing [K] to upper triangular form and then solving for unknowns in the reverse order by back substitution. 3.6 Structural Dynamics An integral part of this study is the incorporation of a structural dynamics response algorithm into the finite element method. These two coupled methods will produce a seamless method by which to analyze a structural system, using the finite element method, under an applied dynamic load(s) for a specified period of time. As a result, the response of the structural system as a function of time will be determined. The following sections explain the theoretical development of the structural dynamics methodology used in this study. Structural dynamics derivations and solution techniques are outlined briefly. The equations are derived for timedependent loads. This study only considers undamped structures, thus, the formulation of damping is not presented and it is taken as zero in all equations. The following sections introduce the theoretical explanations of structural dynamics that were used in the development of the coupled computer program. 3.6.1 Equation of Motion The equation of motion is the basic and fundamental part of structural dynamics. All formulations are derived based on the equation of motion in structural dynamics. The equation of motion is generally given as bü 5 d3e 5 L3 f"N$ (323) 22 The right side of the equation is the time dependent force, p(t). On the left side of equation 323, m represents mass, c represents damping, k represents stiffness and u, 3e, ü represent displacement, velocity and acceleration, respectively (the overdot represents a time derivative). The structure system, which has multiple DOFs, is of course considered in this study. The specific number of DOF depends on the finite element model of the structural system. As a result, all variables in equation 323 are matrices or column vectors of size related to the number of DOF for the system. 3.6.2 Formulation of Element Mass Matrices The formulation of the element mass matrix is based on the virtual work principal and is similar to the formulation of the element stiffness matrix as discussed in Section 3.3. The work done by externally applied loads is equal to the sum of the work absorbed by inertial, dissipative, and internal forces for any virtual displacement. For an element volume V and surface S < =3 B @A 5 < =3 Φ @C 5 Σ =3 f <" =3 g ü 5 =3 d 3e 5 => ? $@A (324) Where B and Φ are prescribed body forces and surface tractions, f and =3 are prescribed concentrated loads and their corresponding virtual displacements, g is mass density and c is the damping coefficient. Following Section 3.3, the nodal displacements, nodal velocities, nodal accelerations, and strains are approximated by 3 @ 3e h@ei ü h@ji > F @ (325) Substitution of equation 325 into equation 324 produces the virtual work expression 23 =@ k< g @Ah@ji 5 < d @Ah@ei 5 < F ? @A / JB@A/ JΦ@C/ l 1mfl 0 (326) The first integral in equation 326 provides the element mass matrix as b < g @A (327) For twodimensional analysis, equation 327 yields b < < g N@ @; 4 !4 9 !9 (328) Following the procedure for assembly of element stiffness matrix into the system matrix, the element mass matrix for each individual element is assembled in the same procedure. The DOF for the mass and stiffness matrices are identical. Thus, from a numerical perspective, the element stiffness and mass matrices can be assembled into the system stiffness and mass matrices simultaneously. 3.6.3 Dynamic Loading Conditions There are many types of dynamic loads from a realistic perspective. In this study, several common idealized dynamic loadings are considered. Namely, these include the harmonic loading, the impulse loading and the step loading. It should be noted, however, that the MATLAB program written as a consequence of this study can handle any userinput for the dynamic loading. The user only has to provide the value of the load at specified time intervals, or sampling points. These three idealized loads are only chosen based on their commonality in structural dynamic simplifications. 24 Harmonic loading is a function of the sine or cosine functions and its equation is in general f"N$ f nlmoN or f drnoN (329) where p0 is the amplitude or maximum value of the force and its frequency o is the exciting frequency or forcing frequency. Figure 32 provides an example of a harmonic loading. Examples of a harmonic loading includes wind driven loading on structures, earthquake (highly idealized) loading, and vehicular motion on bridges. Figure 32: Harmonic Loading A very large force that acts for a very short time but with a time integral that is finite is called an impulsive force. In general, an impulsive force is defined by f"N$ 1/> (330) with a time duration > starting at the time instant t = t, also called the time lag. Figure 3 3 is an illustration of an impulse loading. Examples of an impulsive load primarily include blast loadings such as those due to detonation of blast devices. 25 Figure 33: Impulse Loading Finally, another typical dynamic loading is a force that jumps suddenly from zero to magnitude p0 and stays constant at that value is called step force. In general, the step loading is defined by f"N$ f (331) Figure 34 provides an example of an impulse loading that jumps to its magnitude p0 at time 0 seconds as indicated. Examples of a step loading include the sudden application of a full load rather than being applied gradually over time. Figure 34: Step Loading 26 3.6.4 Numerical Integration There are several numerical integration methods for linear systems to numerically integrate the equation of motion previously defined as equation 323. In all methods for the purposes of this study, initial displacements and velocities are taken zero (the system is initially at rest in an underformed position) and p(t) is known at all time intervals t. The aim of numerical integration is analyzing the system over time intervals ΔN. Some of the numerical integration methods, that are called timestepping methods, are interpolation of excitation, central difference method and Newmark’s method. Interpolation of excitation uses recurrence formulas. It can be used for small ΔN and for linear systems. In addition, this method is suitable for single DOF systems but it is not appropriate for multi DOF systems. Central difference methods are based on finite difference approximations of the time derivatives of displacement, velocity and acceleration. Solution at 3 6 which represents the displacement at the step i+1 is determined from the equation of motion at time step i. Furthermore, 3 and 3 ! must be known to find the displacement at time step i+1, 3 6 . Newmark’s method, used in this study, is a family of timestepping methods based on the iteration of ü 6 . A detailed discussion is presented in the next section. 27 3.6.5 NewmarkBeta Solution Technique Newmark’s technique, also known as the more general NewmarkBeta solution technique, is based on the following equations 3e 6 3e 5 "1 / v$ΔN ü 5 "vΔN$ü 6 (332) 3 6 3 5 "ΔN$3e 5 "0.5 / z$"ΔN$ ü 5 z"ΔN$ ü 6 (333) The parameters z and v define the variation of acceleration over a time step and determine the stability and accuracy characteristics of the technique. Typical selection of v is 1/2, and z can vary between 1/4 and 1/6. For the average acceleration method, used in this study, z is taken as 1/4. All equations are matrix equations for multi DOF as is the case in this study. The NewmarkBeta solution procedure can be presented step by step for linear, multi DOF systems. Initial conditions are defined as 3 3"0$ and 3e 3e"0$ representing the initial displacement and initial velocity, respectively. Initial calculations are only calculated one time and they are presented in equations 334 thorough 337 as Solving for the initial acceleration ü ; bü f / d3e / L3 (334) After selectingΔN and the NewmarkBeta parameters z and v the effective stiffness of the system can be calculated as L{ L 5  }Δ~ d 5 }Δ~,b (335) Calculation constants, a and b, for use at each time step are calculated as 28 }Δ~b 5  } d (336) }b 5 ΔN"  } / 1$d (337) For each time step, equations 338 thorough 342 are be repeated until all time steps are done. The effective force at time step i is then calculated as Δf̂ Δf 5 3e 5 ü (338) Solving for the change in displacement during the time step Δ3 is then L{ Δ3 Δf̂ (339) During the same time step, the change in velocity and the change in acceleration are also Δ3e  }Δ~ Δ3 /  } 3e 5 ΔN 1 /  } ü (340) Δü }Δ~, Δ3 / }Δ~ 3e / } ü (341) Finally, updating the displacement, velocity and acceleration can be found relative to the previous position, velocity and acceleration as 3 6 3 5 Δ3 3e 6 3e 5 Δ3e (342) ü 6 ü 5 Δü Repeating these calculations, equations 338 through 342, the displacements, velocities and accelerations can be found for discrete time points for the structural system. 29 CHAPTER IV CASE STUDIES 4.1 Introduction Case studies are provided to illustrate the capabilities of the MATLAB FEA program. This code has one main program file, one input file and multiple subroutines to analyze the structural system. Material properties, geometrical properties, boundary conditions, loading conditions and more can be easily defined by the user in the input file of the MATLAB code. The case studies all consider twodimensional dynamically loaded beams modeled using a finite element mesh of Q4 elements. To confirm the procedure and accuracy, a basic example is solved by hand in Appendix B of this study. Additionally, the main program and input file of the MATLAB code is provided in Appendix A of this study. Table 41 shows the beams properties that are used in the four case studies. E represents modulus of elasticity; and g represent Poisson's ratio and the mass density of the beam, respectively. In the first case study, the beam is fixed at the left end. The second case study considers a beam that is simply supported. Finally, the third and fourth case studies both consider continuous beams that are supported at the left end, the right end, and at the 30 midpoint. The geometry of the beams, boundary conditions and loading conditions are provided with the case studies. Modulus of Elasticity (psi) Poisson’s Ratio Mass Density (lbs/in3) Number of Elements Number of Nodes Type of Loading Beam Length (in) Beam Depth (in) Case 1 4.35(106) 0.15 0.0868 48 65 Harmonic 144 12 Case 2 4.35(106) 0.15 0.0868 48 65 Impulse 144 12 Case 3 4.35(106) 0.15 0.0868 57 80 Step 240 12 Case 4 4.35(106) 0.15 0.0868 57 80 Multiple 240 12 Table 41: Case Studies As seen in Table 41, four types of loading cases are considered. The system is undamped for all cases. Because twodimensional analyses are considered, each node has two degrees of freedom in the x and y directions while the beam width is taken as unity for all cases. The dynamic analysis duration is 10 seconds, the time step t is taken as 0.1 seconds and initial displacements and velocities for all DOF are zero at time = 0 seconds. With the given time duration and time step, each case study has 100 time steps to solve in the dynamic response analysis. The MATLAB FEA program has the capability to produce movies of the displacements during the time duration of the dynamic loading. As it is not possible to place a movie in a text format, select “snapshots” have been provided for each case study to illustrate the dynamic response of the structural system. These figures serve only as illustrations as to the full capabilities of the MATLAB FEA program developed as a consequence of this study. 31 4.2 Harmonic Loading A harmonic loading is applied to the first beam which is cantilever beam. The illustration of the beam is given in the Figure 41. The applied load F represents the harmonic load as it is applied at the right side of the beam. The properties of the beam are provided in Table 41 as Case 1. Figure 41: Case Study 1 – Cantilever Beam The harmonic load is a sine function and is illustrated in Figure 42. The harmonic load is applied between 0 and 10 seconds with sampling points taken at intervals of 0.1 seconds. The magnitude of the load varies between 1800 lbs and 1800 lbs. 32 Figure 42: Case Study 1 – Harmonic Loading As stated in Table 41, the beam is formulated with 48 rectangular Q4 elements having 130 DOFs. Illustration of the finite element mesh of the beam is shown in Figure 43. Green circles indicate node numbering and red circles indicate element numbering in the mesh. Figure 44 illustrates the displacement of the beam considering a mesh of Q4 finite elements. Subfigures (a) – (f) provide a representation of the movie generated at discrete time intervals as indicated on each subfigure. 0 1 2 3 4 5 6 7 8 9 10 1500 1000 500 0 500 1000 1500 Time [seconds] Load [lbs] Harmonic Load 33 Figure 43: Case Study 1 – Finite Element Mesh Figure 44(a): Case Study 1 – Displacements 0 20 40 60 80 100 120 140 2 0 2 4 6 8 10 12 14 1 1 2 14 15 2 3 16 3 4 17 4 5 18 5 6 19 6 7 20 7 8 21 8 9 22 9 10 23 10 11 24 11 12 25 12 13 26 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 40 41 26 29 42 27 30 43 28 31 44 29 32 45 30 33 46 31 34 47 32 35 48 33 36 49 34 37 50 35 38 51 36 39 52 37 53 54 38 55 39 56 40 57 41 58 42 59 43 60 44 61 45 62 46 63 47 64 48 65 length [in] depth [in] Finite Element Mesh 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1 seconds 34 Figure 44(b): Case Study 1 – Displacements Figure 44(c): Case Study 1 – Displacements 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.7 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =3.5 seconds 35 Figure 44(d): Case Study 1 – Displacements Figure 44(e): Case Study 1 – Displacements 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =5 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =8.1 seconds 36 Figure 44(f): Case Study 1 – Displacements As seen from the subfigures of Figure 44, the free end of the beam moves up and down due to the harmonic loading. Comparison of hand calculations of the steadystate displacement response and computer results are presented in Appendix B. DOF 78 at node 39 is considered while comparing the results of displacements in Appendix B. 4.3 Impulse Loading An impulse loading is applied to the 12ft long simplysupported beam. The load is applied at the middle of the beam. The illustration of the beam is shown in Figure 45. The applied load F represents the impulse load in the figure. The properties of the beam are provided in Table 41 as Case 2. 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =10 seconds 37 F Length = 12' Depth = 1' 6' 6' Figure 45: Case Study 2 – SimplySupported Beam The impulse loading, as shown in Figure 46, is applied between 0 and 10 seconds with sampling points taken at intervals of 0.1 seconds. At time = 1s the impulse loading immediately jumps to a value of –2500 lbs. At all other times t, the loading is zero. Figure 46: Case Study 2 – Impulse Loading 0 1 2 3 4 5 6 7 8 9 10 2.5 2 1.5 1 0.5 0 x 10 4 Time [seconds] Load [lbs] Impulse Load 38 As stated in Table 41, the beam, which has 130 DOFs, is constructed with 48 rectangular Q4 elements. Illustration of the finite element mesh of the beam is shown in Figure 47. Figure 48 illustrates the displacement of the beam considering a mesh of Q4 finite elements. Subfigures (a) – (f) provide a representation of the movie generated at discrete time intervals as indicated on each subfigure. Figure 47: Case Study 2 – Finite Element Mesh 0 20 40 60 80 100 120 140 2 0 2 4 6 8 10 12 14 1 1 2 14 15 2 3 16 3 4 17 4 5 18 5 6 19 6 7 20 7 8 21 8 9 22 9 10 23 10 11 24 11 12 25 12 13 26 13 27 28 14 29 15 30 16 31 17 32 18 33 19 34 20 35 21 36 22 37 23 38 24 39 25 26 27 28 29 30 31 32 33 34 35 36 37 40 41 53 54 38 42 55 39 43 56 40 44 57 41 45 58 42 46 59 43 47 60 44 48 61 45 49 62 46 50 63 47 51 64 48 52 65 length [in] depth [in] Finite Element Mesh 39 Figure 48(a): Case Study 2 – Displacements Figure 48(b): Case Study 2 – Displacements 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.1 seconds 40 Figure 48 (c): Case Study 2 – Displacements Figure 48(d): Case Study 2 – Displacements 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.3 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.4 seconds 41 Figure 48(e): Case Study 2 – Displacements Figure 48(f): Case Study 2 – Displacements Until time = 1 second there is no displacement at the beam as expected. At that time, the impulse load is applied to the beam. Because the system is undamped, after 1 second the 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =5 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =10 seconds 42 beam starts oscillating up and down and goes forever with the maximum displacement at the midpoint as seen from the Figure 48 and its subfigures. 4.4 Step Loading A step loading is applied to the third case study. The 20ft long beam has a pin support at the left end and roller supports at the middle and right end. The load is located 5 feet from the right side of the beam. The illustration of the beam is shown in Figure 49. The applied load F represents the step load in the figure. The properties of the beam are provided in Table 41 as Case 3. 10' F 5' 5' Length = 20' Figure 49: Case Study 3 – Continuous Beam The step loading, as shown in Figure 410, is applied between 0 and 10 seconds with sampling points taken at intervals of 0.1 seconds. As with the impulse loading of Case 2, 43 there is no load applied before the time reaches 1 second. At this point in time, the value of the load jumps suddenly to –20000 lbs and stays constant for the duration of time. Figure 410: Case Study 3 – Step Loading As stated in Table 41, the beam is comprised of 57 rectangular Q4 elements having 160 DOFs. Illustration of the finite element mesh of the beam is shown in Figure 411. Figure 412 illustrates the displacement of the beam considering a mesh of Q4 finite elements. Subfigures (a) – (f) provide a representation of the movie generated at discrete time intervals as indicated on each subfigure. 0 1 2 3 4 5 6 7 8 9 10 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 x 10 4 Time [seconds] Load [lbs] Step Load 44 Figure 411: Case Study 3 – Finite Element Mesh Figure 412(a): Case Study 3 – Displacements 0 50 100 150 200 2 0 2 4 6 8 10 12 14 1 1 2 21 22 2 3 23 3 4 24 4 5 25 5 6 26 6 7 27 7 8 28 8 9 29 9 10 30 10 11 31 11 12 32 12 13 33 13 14 34 14 15 35 15 16 36 16 17 37 17 18 38 18 19 39 19 20 40 20 41 42 21 43 22 44 23 45 24 46 25 47 26 48 27 49 28 50 29 51 30 52 31 53 32 54 33 55 34 56 35 57 36 58 37 59 38 60 39 61 62 40 63 41 64 42 65 43 66 44 67 45 68 46 69 47 70 48 71 49 72 50 73 51 74 52 75 53 76 54 77 55 78 56 79 57 80 length [in] depth [in] Finite Element Mesh 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1 seconds 45 Figure 412(b): Case Study 3 – Displacements Figure 412 (c): Case Study 3 – Displacements 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.2 seconds 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =3.5 seconds 46 Figure 412(d): Case Study 3 – Displacements Figure 412(e): Case Study 3 – Displacements 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =5.2 seconds 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =7 seconds 47 Figure 412(f): Case Study 3 – Displacements Because of the step load, the beam displacement never becomes a positive displacement (with respect to the coordinate system used) between the roller supports and never becomes a negative displacement between the pin and roller supports. These results are as expected. 4.5 Multiple Dynamic Loads The final case study considers the same beam that was used for the third case study when the step loading was considered. However, the loading condition is different in this case. F1 represents a harmonic load and F2 represents a step load. The harmonic load is located 5ft away from the pin support and the step load is located 5ft away from the left end of the beam as shown in the Figure 413. 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =10 seconds 48 F2 5' 5' Length = 20' F1 5' 5' Figure 413: Case Study 4 – Continuous Beam The harmonic load is applied to the beam from time = 0 seconds until 10 seconds and its magnitude varies between 15000 lbs to 15000 lbs. At time = 4 seconds, suddenly another load is applied, the step load, and it stays constant for the duration of time with a magnitude of –20000 lbs. These two timedependent loads are shown the Figure 414. 49 Figure 414: Case Study 4 – Multiple Loading As stated in Table 41, the beam is comprised of 57 rectangular Q4 elements and having 160 DOFs. The finite element mesh of the beam is illustrated in the Figure 415. Figure 416 illustrates the displacement of the beam considering a mesh of Q4 finite elements. Subfigures (a) – (f) provide a representation of the movie generated at discrete time intervals as indicated on each subfigure. 0 1 2 3 4 5 6 7 8 9 10 2 1.5 1 0.5 0 0.5 1 x 10 4 Time [seconds] Load [lbs] Multiple Load harmonic load step load 50 Figure 415: Case Study 4 – Finite Element Mesh Figure 416(a): Case Study 4 – Displacements 0 50 100 150 200 2 0 2 4 6 8 10 12 14 1 1 2 21 22 2 3 23 3 4 24 4 5 25 5 6 26 6 7 27 7 8 28 8 9 29 9 10 30 10 11 31 11 12 32 12 13 33 13 14 34 14 15 35 15 16 36 16 17 37 17 18 38 18 19 39 19 20 40 20 41 42 21 43 22 44 23 45 24 46 25 47 26 48 27 49 28 50 29 51 30 52 31 53 32 54 33 55 34 56 35 57 36 58 37 59 38 60 39 61 62 40 63 41 64 42 65 43 66 44 67 45 68 46 69 47 70 48 71 49 72 50 73 51 74 52 75 53 76 54 77 55 78 56 79 57 80 length [in] depth [in] Finite Element Mesh 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1 seconds 51 Figure 416(b): Case Study 4 – Displacements Figure 416(c): Case Study 4 – Displacements 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =3 seconds 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =3.8 seconds 52 Figure 416(d): Case Study 4 – Displacements Figure 416(e): Case Study 4 – Displacements 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =4.1 seconds 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =5 seconds 53 Figure 416(f): Case Study 4 – Displacements As shown in Figure 416, the influence of the step load when applied at time = 4 seconds is readily distinguishable. Before the application of the step load, displacements are limited due to the harmonic load. However, clearly the displacement of the beam before 4 seconds is harmonic in nature as illustrated. After a period of time after 4 seconds, the displacement again begins to resemble a harmonic pattern as expected. 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =9 seconds 54 CHAPTER V CONCLUSION 5.1 Summary and Conclusions The primary objective of this study was to develop a MATLAB based computer program coupling the finite element method and structural dynamics analysis techniques. The primary reason for using the finite element method is its significant advantages in engineering problems. It is a powerful solution technique for differential and integral equations in complex engineering problems. When combined with a structural dynamic solution algorithm, the capabilities of the two become very unique. Theoretical equations are derived for both the finite element method and structural dynamics during this study. Fundamentally, formulation of element matrices, integration techniques, boundary conditions, and timestepping methods are discussed and their equations are presented. These concepts are used to develop the computer code. Several case studies are analyzed during this study with different types of dynamic loads including harmonic loading, impulse loading, and step loading. These loads are also applied to a variety of different beams including a cantilever beam, a simplysupported 55 beam and a continuous beam with three supports. Displacements and deformed shapes are found using MATLAB FEA program written as a consequence of this study. This study presents a discussion of the finite element method for structural dynamics applications. It provides a basic understanding of the behavior of structural systems under different dynamic loads. The computer code, which is written in MATLAB, allows for further flexibility not explicitly discussed in this study such as using different material properties throughout the structure, a variety of geometrical shapes of the structure, varying types of loads, and multiple finite element types beyond the Q4 element used in the finite element mesh. 5.2 Recommendations This study is open to development and further development is expected in the future. Although currently it is limited to Q4 elements for beams, it can be expanded to advanced finite elements and other types of structures (see Iqbal, 2009). The MATLAB code is written to be expanded to allow for the analysis of more complicated structural systems utilized more advanced finite elements and finite element formulations. Moreover, one can develop this study in not only finding displacements, deformed shapes, stresses, and strains but also other topics such as calculating fracture, durability or even reliability of structures. 56 Although the computer code has good properties, it is also recommended that better meshing capabilities may be implemented in the MATLAB code. The computer code covers linear programming for now. Therefore, the computer program can be incorporated with nonlinearity. If these modifications can be done, this computer code can become a more powerful finite element analysis program. Work done in conjunction with this study includes the incorporation of advanced finite element formulations beyond the Q4 element used herein (Iqbal, 2009). 57 REFERENCES [1] Bozdogan, K.B., Ozturk, D., and Nuhoglu, A. (2009), “An Approximate Method for Static and Dynamic Analyses of Multibay Coupled Shear Walls”, The Structural Design of Tall and Special Buildings, Vol. 18, pp. 112. [2] Cerioni, R., Brighenti, R., and Donida, G. (1995), “Use of Incompatible Displacement Modes in a Finite Element Model to Analyze the Dynamic Behavior of Unreinforced Masonry Panels”, Computers and Structures, Vol. 57, pp. 47  57. [3] Chopra, K.A. (2007), “Dynamics of Structures Theory and Applications to Earthquake Engineering ”, Pearson Prentice Hall. [4] Cook, R.D., Malkus, D.S., Plesha, M.E. and Witt, R.J. (2002), “Concepts and Applications of Finite Element Analysis”, John Wiley & Sons. [5] Du, H., Hitchings, D., and Davies, G.A.O. (1992), “A Finite Element Structural Dynamics Model of a Beam with an Arbitrary Moving Base – Part I: Formulations”, Finite Elements in Analysis and Design, Vol. 12, pp. 133  150. [6] Kiral, Z., and Kiral, B.G. (2008). “Dynamic Analysis of a Symmetric Laminated Composite Beam Subjected to a Moving Load with Constant Velocity”, Journal of Reinforced Plastics and Composites, Vol. 27, pp. 1932. 58 [7] Li, Q.S., Yang, L.F., and Li, G.Q. (2001), “The Quadratic Finite Element/Strip with Generalized Degrees of Freedom and Their Application”, Finite Elements in Analysis and Design, Vol. 37, pp. 325  339. [8] Mermertas, V. (1998), “Dynamic Interaction Between the Vehicle and Simply Supported Curved Bridge Deck”, Computer Methods in Applied Mechanics and Engineering, Vol. 162, pp. 125  131. [9] Ozcan, D.M., Bayraktar, A., Sahin, A., Haktanir, T., and Turker, T. (2008), “Experimental and Finite Element Analysis on the Steel FiberReinforced Concrete (SFRC) Beams Ultimate Behavior. Construction and Building Materials, Vol. 23, pp. 1064  1077. [10] Iqbal, R. (2009), “Development of a Finite Element Program Incorporating Advanced Element Types”, MS Thesis, Oklahoma State University. 59 APPENDIX A MATLAB FEA PROGRAM The Matlab FEA program that was developed during the course of this study is provided herein. First, the main FEA program (not containing the subroutines since the FEA program is still under development) is provided. Second, the general input file used to setup an analysis is provided. Both of these files are Matlab script files and can be executed in all versions of Matlab. To generate graphical output, the full version of Matlab is required. Main Program % ************************************************************************* % Finite Element Analysis (FEA) Program to Determine the Structural and Thermal Response of Structural Systems % Written By: Muhammet Saglar & Rameez Iqbal (Advised By: Dr. Jonathan S. Goode) % School of Civil and Environmental Engineering % Oklahoma State University % ************************************************************************* % ************************************************************************* % ************************************************************************* % Begin Program % Clear Variables and Settings clear all % Clear Screen clc; % Inhibit Warning Messages warning('off','all'); disp('**************************************************'); disp('Finite Element Analysis') disp('Written by: Saglar/Iqbal (Goode)'); disp('School of Civil and Environmental Engineering'); disp('Oklahoma State University'); disp('**************************************************'); 60 disp('**************************************************'); disp('Running Analysis .....'); disp('**************************************************'); % Begin Timer tic; % ************************************************************************* % PreProcessing Section % ************************************************************************* disp('PreProcessing .....'); % Call MFile for User Input Data from fea_input.m disp('..... Obtaining Input .....'); fea_input; % Call MFile for Initialization of Variables from fea_initialize.m disp('..... Initializing Program .....'); fea_initialize; % ************************************************************************* % End PreProcessing Section % ************************************************************************* % ************************************************************************* % Processing Section % ************************************************************************* disp('Processing .....'); %  % Assembly Process disp('..... Assembly Process .....'); % Determine Element Property Matrices (Constitutive, Mass, Conducitivity) [property_matrix_k,property_matrix_m,property_matrix_t] = fea_property_matrix(analysis_type_material,material_E,material_v,material_rho,material_t); % Determine Gauss Quadrature Sampling Points and Weights [gauss_points,gauss_weights] = fea_gauss_2d(element_type); % Initialize Global Index of DOF index_global = zeros(edof,nel); % Assemble System Matrices According to Analysis if analysis_type == 'S' % Static Analysis disp('..... ..... Static Analysis .....'); disp('..... ..... ..... Building Stiffness .....'); % Loop Through All Elements to Determine System Matrices (Stiffness, Mass, and Conductivity) % NOTE: Depending on the Analysis Type, Some Matrices Will Not be Generated for j = 1:nel % Call SubRoutine to Extract Nodes for the jth Element [element_nodes,element_xcoord,element_ycoord] = fea_node_coordinates(j,nnel,node_connectivity,node_coord); % Call SubRoutine to Extract System DOFs for the jth Element index = fea_index_dof(element_nodes,nnel,ndof,edof); % Save Index of DOF to Global Variable index_global(:,j) = index; % Call SubRoutine to Generate Element Stiffness Matrix 61 element_k = fea_element_matrix_stiffness(element_type,gauss_points,gauss_weights,element_xcoord,element_ycoord,nnel,edof,property_matrix_ k); % Call SubRoutine to Assemble System Stiffness Matrix system_k = fea_element_matrix_assemble(system_k,element_k,index,edof); end elseif analysis_type == 'D' % Dynamic Analysis disp('..... ..... Dynamic Analysis .....'); disp('..... ..... ..... Building Stiffness and Mass .....'); % Loop Through All Elements to Determine System Matrices (Stiffness, Mass, and Conductivity) % NOTE: Depending on the Analysis Type, Some Matrices Will Not be Generated for j = 1:nel % Call SubRoutine to Extract Nodes for the jth Element [element_nodes,element_xcoord,element_ycoord] = fea_node_coordinates(j,nnel,node_connectivity,node_coord); % Call SubRoutine to Extract System DOFs for the jth Element index = fea_index_dof(element_nodes,nnel,ndof,edof); % Save Index of DOF to Global Variable index_global(:,j) = index; % Call SubRoutine to Generate Element Stiffness Matrix element_k = fea_element_matrix_stiffness(element_type,gauss_points,gauss_weights,element_xcoord,element_ycoord,nnel,edof,property_matrix_ k); % Call SubRoutine to Assemble System Stiffness Matrix system_k = fea_element_matrix_assemble(system_k,element_k,index,edof); % Call SubRoutine to Generate Element Mass Matrix element_m = fea_element_matrix_mass(element_type,gauss_points,gauss_weights,element_xcoord,element_ycoord,nnel,edof,property_matrix_m); % Call SubRoutine to Assemble System Mass Matrix system_m = fea_element_matrix_assemble(system_m,element_m,index,edof); end end % End Assembly Process %  %  % Boundary and Loading Conditions disp('..... Applying Boundary and Loading Conditions .....'); % Apply Boundary and Loading Conditions According to Analysis if analysis_type == 'S' % Static Analysis % Call SubRoutine to Apply Boundary and Loading Conditions [system_k,system_f] = fea_system_boundary_load(system_k,system_f,boundary_dof,boundary_dof_val,load_dof,load_dof_val); elseif analysis_type == 'D' % Dynamic Analysis % Boundary and Loading Conditions Must Be Applied at Each Point in Time During the Dynamic Analysis % See the Analysis Section Below 62 end % End Boundary Conditions %  %  % Analysis disp('..... Solving .....'); % Based on Analysis Type if analysis_type == 'S' % Static Analysis disp('..... ..... Static Analysis .....'); % Solve System Equations for Static Displacements system_d = system_k \ system_f; elseif analysis_type == 'D' % Dynamic Analysis disp('..... ..... Dynamic Analysis .....'); % Call SubRoutine to Solve System Equations for Dynamic Displacements [system_d,system_d_vel,system_d_acc,system_f,dynamic_time] = fea_dynamic_analysis(system_k,system_m,boundary_dof,boundary_dof_val,load_dof,load_file,dynamic_delta_t,dynamic_beta,dyna mic_gamma,sdof); end % End Analysis %  % ************************************************************************* % End Processing Section % ************************************************************************* % ************************************************************************* % PostProcessing Section % ************************************************************************* disp('PostProcessing .....'); %  % Plot Generation % Plot Finite Element Mesh if plot_element_mesh == 'Y' disp('..... Plotting Element Mesh .....'); % Call SubRoutine to Plot Finite Element Mesh plot_handle = fea_plot_mesh(proj_name,nel,node_coord,node_connectivity,plot_element_mesh_numbers,plot_element_mesh_node_numbers); end % Based on Analysis Type if analysis_type == 'S' % Static Analysis disp('..... Static Analysis .....'); % Plot Static Displacements if plot_static_displacement == 'Y' 63 disp('..... ..... Plotting Static Displacements .....'); % Call SubRoutine to Plot Static Displacements plot_handle = fea_plot_static_displacement(proj_name,nel,node_coord,node_connectivity,system_d,index_global); if plot_static_contour_displacement == 'Y' % Call SubRoutine to Plot Displacement Contours plot_handle = fea_plot_static_displacement_contour(proj_name,node_coord,node_connectivity,system_d,sdof,nel); end end % Determine Static Stresses and Strains for Plots if plot_static_stress == 'Y'  plot_static_strain == 'Y' disp('..... ..... Determining Static Stresses and Strains .....'); % Call SubRoutine to Determine Stresses and Strains from Displacements [system_stress,system_strain,gauss_points_coord] = fea_stress_strain(nel,nnel,edof,node_connectivity,node_coord,index_global,system_d,property_matrix_k,element_type,gauss_points); end % Plot Static Stresses if plot_static_stress == 'Y' disp('..... ..... Plotting Static Stresses .....'); % Call SubRoutine to Plot Static Stresses plot_handle = fea_plot_static_stress(proj_name,nel,node_coord,node_connectivity,system_stress,gauss_points_coord); if plot_static_contour_stress % Call SubRoutine to Plot Stress Contours [plot_handle] = fea_plot_static_stress_contour(proj_name,nel,node_coord,node_connectivity,system_stress,gauss_points_coord); end end % Plot Static Strains if plot_static_strain == 'Y' disp('..... ..... Plotting Static Strains .....'); % Call SubRoutine to Plot Static Strains plot_handle = fea_plot_static_strain(proj_name,nel,node_coord,node_connectivity,system_strain,gauss_points_coord); if plot_static_contour_strain % Call SubRoutine to Plot Strain Contours [plot_handle] = fea_plot_static_strain_contour(proj_name,nel,node_coord,node_connectivity,system_strain,gauss_points_coord); end end elseif analysis_type == 'D' % Dynamic Analysis disp('..... Dynamic Analysis .....'); % Call SubRoutine to Plot Dynamic Displacements (Movie) 64 plot_handle = fea_plot_dynamic_displacement(proj_name,nel,node_coord,node_connectivity,system_d,dynamic_time,index_global); disp('..... ..... Plotting Dynamic Displacements .....'); end % ************************************************************************* % End PostProcessing Section % ************************************************************************* disp('**************************************************'); disp('Analysis Done .....'); disp('**************************************************'); % End Timer analysis_time = toc; disp(['Total Time of Analysis: ',num2str(analysis_time),' seconds']); disp('**************************************************'); % End Program % ************************************************************************* % ************************************************************************* Input File % ************************************************************************* % Finite Element Analysis (FEA)  Input File % Written By: Muhammet Saglar & Rameez Iqbal (Advised By: Dr. Jonathan S. Goode) % School of Civil and Environmental Engineering % Oklahoma State University % ************************************************************************* % ************************************************************************* % ************************************************************************* % Begin Input File %  % Project Name % Define Project Name % NOTE: Automatically Saved Figures and Data Will Be Proceeded with the Project Name proj_name = 'test'; %  % Analysis Options % Define Analysis Type % Analysis Type Options: % Input 'S' = Static Analysis % Input 'D' = Dynamic Analysis % Input 'T' = Thermal Analysis analysis_type = 'S'; %  % Input Finite Element Mesh Properties % Input Element Type % Input 'Q4' for Bilinear Rectangular Element % Input 'Q8' for Quadratic Rectangular Element % Input 'CST' for Constant Strain / Linear Triangular Element % Input 'LST' for Linear Strain / Quadratic Triangular Element 65 element_type = 'Q4'; % Input Node Coordinates (x and y coordinates) [Note Units => Length] node_coord = [ 0.0 0.0 ;... 12.0 0.0 ;... 24.0 0.0 ;... 36.0 0.0 ;... 48.0 0.0 ;... 60.0 0.0 ;... 72.0 0.0 ;... 84.0 0.0 ;... 96.0 0.0 ;... 108.0 0.0 ;... 120.0 0.0 ;... 0.0 12.0 ;... 12.0 12.0 ;... 24.0 12.0 ;... 36.0 12.0 ;... 48.0 12.0 ;... 60.0 12.0 ;... 72.0 12.0 ;... 84.0 12.0 ;... 96.0 12.0 ;... 108.0 12.0 ;... 120.0 12.0 ]; % Input Nodal Connectivity for Each Element (CCW from BottomLeft) node_connectivity = [ 1 2 13 12 ;... 2 3 14 13 ;... 3 4 15 14 ;... 4 5 16 15 ;... 5 6 17 16 ;... 6 7 18 17 ;... 7 8 19 18 ;... 8 9 20 19 ;... 9 10 21 20 ;... 10 11 22 21 ]; %  % Input Boundary Conditions % Input DOF Constrained boundary_dof = [1 2 23 24]; % Input Constrained DOF Prescribed Values [Note Units => Length] % NOTE: For Dynamic Analysis, All Prescribed Values MUST BE Zero (0) boundary_dof_val = [0 0 0 0]; %  % Input Material Properties % Material Properties % NOTE: If the property is not being used, enter 0 (zero) % Modulus of Elasticity [Note Units => Force / Length^2] material_E = 1e6; % Poisson's Ratio [Unitless] material_v = 0.3; % Mass Density [Note Units => Mass / Length^3] material_rho = 1; % Define Material Analysis Type % Material Analysis Type Options: % Input '1' = Plane Stress Analysis (2D) % Input '2' = Plane Strain Analysis (2D) % Input '3' = ThreeDimensional Anslysis (3D) analysis_type_material = 1; 66 %  % Input Load Properties % Static Loading (If Applicable) % NOTE: If Static Analysis is not being considered, NO changes are needed if strcmp(analysis_type,'S') == 1 % Input DOF Statically Loaded load_dof = [18 20]; % Input DOF Statically Loaded Prescribed Values [Note Units => Force] load_dof_val = [500 500]; % Dynamic Loading (If Applicable) % NOTE: If Dynamic Analysis is not being considered, NO changes are needed elseif strcmp(analysis_type,'D') == 1 % Input DOF Dynamically Loaded load_dof = [44]; % Input Load TimeSeries File Name [Note Units => Force] % NOTE: First Number Designates the Load at Time = 0 seconds % Column 1 > n = DOFs % Row 1 > n = Load at Time Increments load_file = 'load.dat'; % Input Time Step of Load TimeSeries [seconds] dynamic_delta_t = 0.1; % Define NewmarkBeta Parameters [Unitless] % Average Acceleration Method => Beta = 0.25 & Gamma = 0.5 (Unconditionally Stable) % Linear Acceleration Method => Beta = 0.16667 & Gamma = 0.5 (Conditionally Stable > delta_t / Tn <= 0.551) dynamic_beta = 0.25; dynamic_gamma = 0.5; end %  % Output Options % General Plots % Plot Element Mesh (Y or N) plot_element_mesh = 'Y'; % Plot Element Mesh Options % Display Element Numbers on Element Mesh Plot (Y or N) plot_element_mesh_numbers = 'Y'; % Display Element Node Numbers on Element Mesh Plot (Y or N) plot_element_mesh_node_numbers = 'Y'; % Static Analysis (If Applicable) % NOTE: If Static Analysis is not being considered, NO changes are needed % Plot Static Displacements (Y or N) plot_static_displacement = 'Y'; % Contour Plots of Displacements (Y or N) % NOTE: Must Plot Static Displacements to Plot Contours plot_static_contour_displacement = 'Y'; % Plot Static Stresses (Y or N) plot_static_stress = 'N'; % Contour Plots of Stresses (Y or N) % NOTE: Must Plot Static Stresses to Plot Contours plot_static_contour_stress = 'Y'; 67 % Plot Static Strains (Y or N) plot_static_strain = 'N'; % Contour Plots of Strains (Y or N) % NOTE: Must Plot Static Strains to Plot Contours plot_static_contour_strain = 'Y'; % Dynamic Analysis (If Applicable) % NOTE: If Dynamic Analysis is not being considered, NO changes are needed % Movie of Dynamic Displacements (Y or N) movie_dynamic_displacement = 'Y'; % End Input File % ************************************************************************* % ************************************************************************* 68 APPENDIX B CONFIRMATION CALCULATIONS Appendix B provides confirmation calculations of the element formulation for the Q4 elements as well as a few select time intervals for the structural dynamics response of the structural system. These calculations were done by hand and then compared to the results developed by the Matlab FEA program as provided in Appendix A. 69 70 71 72 73 74 75 76 77 Comparison of Hand Calculation and Computer Results  Case 1 DOF Computer Results (in) Hand Calculation Results (in) u1 7.32E48 0 u2 2.69E34 0 u3 0.003575 0.0036 u4 0.024134615 0.0241 u5 0.00715 0.0071 u6 2.69E34 0 u7 0.007 0.007 u8 0.000184615 0.00018 u9 0.003575 0.0036 u10 0.02405 0.024 u11 0.00015 0.00015 u12 0.000184615 0.00018 Comparison of the Hand Calculation and Computer Results  Case 2 DOF Computer Results (in) Hand Calculation Results (in) u1 1.40E32 0 u2 2.69E34 0 u3 0.021 0.021 u4 0.118 0.118 u5 0.028 0.028 u6 0.376 0.376 u7 1.40E32 0 u8 2.69E34 0 u9 0.021 0.021 u10 0.118 0.118 u11 0.028 0.028 u12 0.376 0.376 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 Time Step : t = 0.1 second t = 0.2 second DOF Hand Calculation Results Computer Results Hand Calculation Results Computer Results u1 0 2.29E28 0 6.24E28 u2 0 9.15E29 0 3.37E28 u3 0.00028 0.00028 0.001 0.00141 u4 0.00037 0.000374 0.001 0.00101 u5 0.00074 0.00075 0.003 0.00342 u6 0.0059 0.0059 0.029 0.02951 u7 0 2.27E28 0 6.31E28 u8 0 1.51E29 0 2.27E28 u9 0.00021 0.00021 0.001 0.00127 u10 0.00042 0.000422 0.0009 0.00093 u11 0.00087 0.000872 0.003 0.003711 u12 0.0061 0.00611 0.0298 0.02983 Time Step: t = 0.3 second t = 0.4 second DOF Hand Calculation Results Computer Results Hand Calculation Results Computer Results u1 0 8.05E29 0 2.87E27 u2 0 4.09E28 0 4.15E29 u3 0.003 0.00365 0.006 0.00672 u4 0.011 0.01184 0.036 0.03605 u5 0.007 0.00701 0.009 0.00935 u6 0.071 0.0719 0.121 0.12166 u7 0 8.65E29 0 2.87E27 u8 0 3.02E28 0 8.43E29 u9 0.003 0.003589 0.006 0.006715 u10 0.011 0.01178 0.036 0.03603 u11 0.007 0.007254 0.009 0.009488 u12 0.072 0.07217 0.121 0.12177 95 96 97 98 Note: The displacement response shown compares the steadystate solution of the handcalculation only. The computer program results include both the steadystate response and the transient response of the system. As a result, some disagreement, particularly at the peaks of the steadystate response, is expected and confirmed in the figure. 0 1 2 3 4 5 6 7 8 9 10 2.5 2 1.5 1 0.5 0 0.5 1 1.5 2 2.5 Time [seconds] Displacement [in] Comparison of Hand Calculation(SteadyState) and Computer Results Hand Calculation Computer VITA Muhammet Ali Saglar Candidate for the Degree of Master of Science Thesis: DEVELOPMENT OF A FINITE ELEMENT ANALYSIS PROGRAM FOR STRUCTURAL DYNAMICS APPLICATIONS Major Field: Civil Engineering Biographical: Personal Data: Born in Balikesir, Turkey, August 22nd 1984, son of Remzi Saglar and Asiye Saglar. Education: Graduated high school from Ataturk Anatolian, Ankara, Turkey in June 2002; received Bachelor of Engineering degree in Civil Engineering from Istanbul Technical University, Istanbul, Turkey in June 2007. Completed the requirements for the Master of Science in Civil Engineering at Oklahoma State University, Stillwater, Oklahoma in July, 2009. Experience: Employed by Oklahoma State University as a research assistant from January 2008 to July 2009. Employed by Tekfen Corporation as a civil engineering intern from May 2006 to July 2006. Employed by Yuksel Construction Corporation as a civil engineering intern from May 2005 to July 2005. Professional Memberships: American Concrete Institute, American Institute of Steel Construction. ADVISER’S APPROVAL: Dr. Jonathan S. Goode Name: Muhammet Ali Saglar Date of Degree: July, 2009 Institution: Oklahoma State University Location: Stillwater, Oklahoma Title of Study: DEVELOPMENT OF A FINITE ELEMENT ANALYSIS PROGRAM FOR STRUCTURAL DYNAMICS APPLICATIONS Pages in Study: 107 Candidate for the Degree of Master of Science Major Field: Civil Engineering Scope and Method of Study: The finite element method is a powerful method to find solutions for engineering problems. Structural dynamics is an important concept for understanding behavior of structures under different types of dynamic loads. The combination of these two concepts into a seamless, integrated computer program to analyze structural systems considering dynamic, timedependent loads provides tremendous capabilities with respect to numerical analysis. The implementation of the theoretical development of both concepts into a fullyfunctional MATLAB computer program provided the primary objective of this study. Findings and Conclusions: The development of the MATLAB computer program to analyze a structural system using the finite element method incorporating structural dynamic response due to timevarying loads was accomplished. Several case studies were presented that illustrate the capabilities of this program with respect to determining the response of the structural system. In the current study, bilinear quadrilateral (Q4) isoparametric finite elements were used to construct the finite element mesh of the structural system. The finite element method was then integrated with a numerical timestepping method, the NewmarkBeta method, to determine the response of the structural system due to a dynamic loading.
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Development of a finite element analysis program for structural dynamics applications 
Date  20090701 
Author  Saglar, Muhammet Ali 
Department  Civil & Environmental Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  The finite element method is a powerful method to find solutions for engineering problems. Structural dynamics is an important concept for understanding behavior of structures under different types of dynamic loads. The combination of these two concepts into a seamless, integrated computer program to analyze structural systems considering dynamic, timedependent loads provides tremendous capabilities with respect to numerical analysis. The implementation of the theoretical development of both concepts into a fullyfunctional MATLAB computer program provided the primary objective of this study.The development of the MATLAB computer program to analyze a structural system using the finite element method incorporating structural dynamic response due to timevarying loads was accomplished. Several case studies were presented that illustrate the capabilities of this program with respect to determining the response of the structural system. In the current study, bilinear quadrilateral (Q4) isoparametric finite elements were used to construct the finite element mesh of the structural system. The finite element method was then integrated with a numerical timestepping method, the NewmarkBeta method, to determine the response of the structural system due to a dynamic loading. 
Note  Thesis 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  DEVELOPMENT OF A FINITE ELEMENT ANALYSIS PROGRAM FOR STRUCTURAL DYNAMICS APPLICATIONS By MUHAMMET ALI SAGLAR Bachelor of Civil Engineering Istanbul Technical University Istanbul, Turkey 2007 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE July, 2009 ii DEVELOPMENT OF A FINITE ELEMENT ANALYSIS PROGRAM FOR STRUCTURAL DYNAMICS APPLICATIONS Thesis Approved: Dr. Jonathan S. Goode Thesis Adviser Dr. Robert N. Emerson Dr. Tyler Ley Dr. A. Gordon Emslie Dean of the Graduate College iii ACKNOWLEDGMENTS I thank God for His grace, His blessings and guiding me to the right path. I would like to express profound gratitude to my advisor, Dr. Jonathan S. Goode, for his precious support, encouragement and useful suggestions throughout this research work. His moral support and continuous guidance enabled me to complete my work successfully. I am also highly thankful to Dr. Robert Emerson and Dr. Tyler Ley for their encouragements guidance rendered to me throughout my research. Their deep experience helped me a lot to be a structural engineer. I would also like to thank to Rameez Iqbal, M.S student, for his support during this research. I am really grateful to him for his friendship and understanding. My deepest appreciation goes to my parents, Remzi Saglar and Asiye Saglar for their invaluable supports during every level of my life. I am grateful also to my three brothers their unflinching courage and conviction. iv TABLE OF CONTENTS Chapter Page I. INTRODUCTION ......................................................................................................1 1.1 Background and Problem Definition .................................................................1 1.2 Objectives ..........................................................................................................2 1.3 Overview ............................................................................................................3 II. LITERATURE REVIEW ..........................................................................................6 2.1 Finite Element Method ......................................................................................6 2.2 Structural Dynamics...........................................................................................8 2.3 MATLAB Usage in Structural Dynamics........................................................10 III. FINITE ELEMENT THEORETICAL DERIVATIONS .......................................12 3.1 Introduction ......................................................................................................12 3.2 Interpolation and Shape Functions...................................................................13 3.3 Formulation of Element Stiffness Matrices .....................................................16 3.4 Element Boundary and Loading Conditions ....................................................17 3.5 Solution Techniques.........................................................................................19 3.6 Structural Dynamics.........................................................................................21 3.6.1 Equation of Motion .................................................................................21 3.6.2 Formulation of Element Mass Matrices ..................................................22 3.6.3 Dynamic Loading Conditions .................................................................23 3.6.4 Numerical Integration .............................................................................26 3.6.5 NewmarkBeta Solution Technique ........................................................27 IV. CASE STUDIES ....................................................................................................29 4.1 Introduction ......................................................................................................29 4.2 Harmonic Loading ...........................................................................................31 4.3 Impulse Loading ..............................................................................................36 4.4 Step Loading ....................................................................................................42 4.5 Multiple Dynamic Loads .................................................................................47 v Chapter Page V. CONCLUSION ......................................................................................................54 5.1 Summary and Conclusions ..............................................................................54 5.2 Recommendations ............................................................................................55 REFERENCES ............................................................................................................57 APPENDIX A – MATLAB FEA PROGRAM ............................................................59 APPENDIX B – CONFIRMATION CALCULATIONS ............................................68 vi LIST OF TABLES Table Page 41 Case Studies .....................................................................................................30 vii LIST OF FIGURES Figure Page 21 Crack Distributions at the Ultimate Load Step (Ozcan, 2008) ..........................8 31 Q4 Finite Element Nodal Degrees of Freedom ................................................15 32 Harmonic Loading ...........................................................................................24 33 Impulse Loading ..............................................................................................25 34 Step Loading ....................................................................................................25 41 Case Study 1 – Cantilever Beam .....................................................................31 42 Case Study 1 – Harmonic Loading ..................................................................32 43 Case Study 1 – Finite Element Mesh ...............................................................33 44 Case Study 1 – Displacements .........................................................................33 45 Case Study 2 – SimplySupported Beam .........................................................37 46 Case Study 2 – Impulse Loading .....................................................................37 47 Case Study 2 – Finite Element Mesh ...............................................................38 48 Case Study 2 – Displacements .........................................................................39 49 Case Study 3 – Continuous Beam ....................................................................42 410 Case Study 3 – Step Loading ...........................................................................43 411 Case Study 3 – Finite Element Mesh ...............................................................44 412 Case Study 3 – Displacements .........................................................................44 413 Case Study 4 – Continuous Beam ....................................................................48 414 Case Study 4 – Multiple Loading ....................................................................49 415 Case Study 4 – Finite Element Mesh ........................................................................... 50 416 Case Study 4 – Displacements .........................................................................50 1 CHAPTER I INTRODUCTION 1.1 Background and Problem Definition Finite element analysis (FEA) is a computational technique used to find approximate solutions of field problems that are described by differential equations or by an integral expression. The finite element solution technique helps in solving complex field problems by numerical approximations of the differential or integral expression. Finite element analysis has significant advantages when compared to the other numerical analysis methodologies. It is very powerful and applicable in many engineering problems such as displacements of a structural systems, stressstrain analysis, heat transfer and magnetic fields. In addition, for individual elements different material properties can be incorporated. One of the most important advantages of finite element analysis is that there are no limitations concerning the geometry or boundary conditions. Different types of geometry and boundary conditions can be accommodated easily. In addition to these properties, it is easy to modify the problem and increase the accuracy of the results while usually only at the expense of computational efficiency. 2 At present many commercial software packages use the finite element analysis because of the advantages mentioned above. ANSYS, ALGOR, ABAQUS, COSMOS/M and SAP2000 are some of the wellknown commercial software packages that are based on finite element analysis. That is to say that the finite element analysis has wide range of usage and is conveniently available for engineers and researchers. For research purposes, finite element analysis is very useful to allow for examination of “what if” design scenarios. As mentioned earlier, it has almost no limitation for engineering problems in terms of geometry and boundary conditions. Because of these properties it is very powerful for researchers. In this study, four types of dynamic loads are applied separately to different beams. The finite element analysis is used to define the beam and determine the problem solution for both loading cases. Furthermore, Q4 finite elements are considered during the analysis and these elements are combined to form a mesh. During the dynamic analysis calculation procedure a numerical procedure, Newmark’s method, is used to numerically integrate the equations of motion with respect to time. 1.2 Objectives The primary objective of this research is the development of a finite element analysis computer program written using MATLAB. This program will be used to determine structural displacements or response of structural systems to both static and dynamic 3 loading conditions. As a method of illustrating the use of this program, several case studies are presented that find the displacements and deformed shape of a beam under the dynamic loading with change of time by using finite element method. The purpose of this study is development of a coupled finite element analysis/structural dynamics computer program. This study does not focus on accuracy of the results. In other words, the best approximation to the exact results is not the main intention. As previously mentioned, the computer codes are written in MATLAB which is a computationally powerful programming language for research purposes. Explained in its most basic form, this computer code has the capability to calculate approximate displacements which can be calculated using design software packages. Besides these capabilities, the development of this program gives the author the understanding of what is going on behind the software packages. This program has the capacity to incorporate different scenarios for dynamic loading. In other words, any timedependent loading scenario to calculate the dynamic load can be applied to the beam to obtain displacements and the deformed shape of the structure. 1.3 Overview The next several chapters present the methodology behind this study. A summary of each chapter is provided to give a brief overview of the remaining sections of this study: • Chapter 2 – Literature Review: 4 Current and recent studies are provided from a broad viewpoint including a brief discussion about the usage of the finite element method and structural dynamics. The literature review gives several examples that use the finite element method, structural dynamics and structural dynamics combined with the finite element method. • Chapter 3 – Finite Element Theoretical Derivations: Detailed theoretical development of the finite element analysis procedure and structural dynamics is presented and discussed. This information is directly used to develop the MATLAB code for the FEA program. Theoretical equations are derived and explained in this chapter. Topics that are covered include discussions about generating element matrices, boundary and loading conditions, finite element solution techniques, and structural dynamics solution technique. • Chapter 4 – Case Studies: The capabilities of the MATLAB FEA program developed as a consequence of this study are demonstrated. Four case studies considering a dynamic loading are demonstrated on different structural beam configurations. Hand calculations (included in Appendix B of this study) are provided for confirmation of the FEA program. 5 • Chapter 5 – Summary, Conclusion and Recommendations: A summary of the results and conclusions of the study are provided. This includes a brief discussion and capability of this study. In addition, recommendations about the study and its expandability are provided. 6 CHAPTER II LITERATURE REVIEW 2.1 Finite Element Method The finite element method is an approximation method that can be used to calculate stresses, movements of loads and forces, displacements, heat transfer and other basic physical behaviors while using very large matrix arrays and mesh diagrams. In recent years, the finite element method has been used to obtain the solutions to a variety of engineering problems. A few applications of finite element modeling that have been used for engineering research are subsequently illustrated. Li et al. (2001) presented a quadratic finite element considering the principle that the local displacement fields of the elements are compatible with the global displacement field of the corresponding systems by using generalized degrees of freedom (GDOF) and a quadratic finite strip with GDOF. They developed a global displacement field based on quadratic Bspline functions and created local displacement field strips for the elements. Therefore, they formed models for the finite elements and finite strips corresponding to the GDOF by rearranging the local displacement fields of elements and strips to be compatible with their corresponding global displacement fields. In their study, several 7 numerical examples were provided for authenticity, simplicity and versatility of the element and strip in the analysis of thinwalled structures. The authors determined that the finite element and finite strip with GDOF can handle inconveniences in the analysis of beams, plates and other structures when the thickness changes. In addition to handling such problems, the authors also determined that the finite element and finite strip gave similar results with fewer degrees of freedom when compared to traditional finite element solutions. Özcan et al. (2008) presented three steel fiberadded reinforced concrete models which were analyzed experimentally and numerically with a finite element analysis. They created a finite element model using ANSYS. The authors considered eightnode solid brick elements in their finite element modeling. Experimental data and finite element analysis results were compared. Figure 21 is given as an illustration for comparison of experimental research and finite element modeling. As a result of their study, the experimental data and finite element modeling gave very close results in terms of deflections and stresses at the center line along with initial and progressive cracking, failure mechanism and load, deflections and stresses at the zero deflection point and decompression. 8 Figure 21: Crack Distributions at the Ultimate Load Step (Ozcan, 2008) 2.2 Structural Dynamics Incorporation of the dynamic response of structural systems into a finite element model is a natural extension of the numerical timestepping methods for dynamic analysis. Recently, many studies have considered this kind of “coupling” of a FEA program and structural dynamics models. A review of some of these studies is provided herein. 9 Du et al. (1992) considered a threedimensional elastic beam with an arbitrary and large moving base with six degrees of freedom by creating a general finite element structural dynamics model. In their study, six degrees of freedom of the beam base may include either a particular arbitrary motion of the base or the coupling of the beam with other structures. The beam can be pretwisted and has a mass center offset from the elasticity center. In the study, the authors derived the equations of motion using the principal of virtual work in the finite element analysis. To simplify the analysis, the authors combined the beam inertia at the end nodes of each element. At the conclusion of the study, it was determined that the model provided flexibility because of the combination of the finite element method and arbitrary base motion variables that were used in multibody dynamics and a fundamental element to solve the dynamic problems of rotating beamlike structures. Cerioni et al. (1995) presented a finite element model that was capable of analyzing the dynamic nonlinear behavior of unreinforced masonry panels in a biaxial stress field by using a nonconforming quadrilateral element. The study derived the equilibrium equations including the inertial and damping actions by a direct stepbystep integration procedure in the time domain known as Newmark’s method. The results were compared with experimental results in terms of displacement, velocity, acceleration, strain and stress. The comparison of the results indicated good performance not only for nonconforming quadrilateral elements but also timedependent and nonlinear analyses. In addition to these results, the model provided a very simple, computationally economical and convenient analysis of complex structural problems. 10 Mermertas et al. (1997) introduced a curved bridge deck to examine the dynamic interaction between a vehicle that is assumed to have fourdegrees of freedom and a curved beam. The finite element method was used to model the curved beam assumed to be simply supported neglecting any damping of the structure. The authors also applied a multipredictorcorrector procedure in conjunction with the Newmark method in the solution of the resulting equations of motion. The authors determined the midspan deflection of the bridge for different vehicle speeds and varying radii of curvature for the curved deck. 2.3 MATLAB Usage in Structural Dynamics The usage of computer programs to use the finite element method and numerical timestepping methods for the structural dynamics has become increasingly important in recent years. Advances in computer technology allow even simple computer systems to model very complex systems with ease and efficiency. A popular program in research studies is MATLAB. This program can be utilized very easily to accomplish both a finite element analysis and structural dynamics analysis due to its friendly matrix manipulations. Several recent research studies have considered MATLAB and they are further outlined here. Kiral et al. (2008) presented a fixedfixed laminated composite beam that was subjected to a concentrated force traveling at a constant velocity to determine the dynamic behavior 11 of a beam. The authors used classical lamination theory in order to create a threedimensional finite element model. In addition, the Newmark integration method to compute the dynamic response was implemented in MATLAB. The dynamic magnification, defined as the ratio between the dynamic and static displacements, was determined from their study. As a result of the study, it was determined that the load velocity and ply orientation may affect the dynamic response significantly. It was also determined that if the total traveling time is equal to the first natural period of the beam, the maximum deflection occurs at the midpoint of the beam. Bozdogan et al. (2009) demonstrated an approximate method using a continuum approach and transfer matrix method for static and dynamic analyses of multibay coupled shear walls. The authors assumed the structure system as a sandwich beam and thus wrote the differential equation whose solution gave the shape functions for each story of the sandwiched beam. The system modes and periods based on the boundary conditions and story transfer matrices, found by the shape functions, can then be calculated. Using MATLAB, a computer program was written to solve numerical examples to show the reliability of the method used. Results were compared with previous work done in the literature that was in good agreement. The authors suggested that their method is appropriate to use on a wide range of structural system applications. 12 CHAPTER III FINITE ELEMENT THEORETICAL DERIVATIONS 3.1 Introduction The primary objective of this study is to develop a finite element analysis program coupled with a structural dynamics response modeling program to determine the dynamic response of a structure due to a dynamic load. The use of the finite element method, however, influences the dynamic response of the structure. A variety of finite elements may be chosen for the analysis, each having its advantages and disadvantages. As such, one element is not always superior to another with respect to any given analysis. Often, it is the experience of the analyst that determines the appropriate finite element to be used in the analysis. The purpose of this study is not to determine the appropriateness of the finite element method chosen. In addition, there are also a variety of structural dynamic response algorithms that can be used to determine the response of a structure due to a dynamic load application. As such, it is also not the purpose of this study to determine the “best” algorithm to be used for determining the dynamic response of the structural system. Rather, this study seeks only to combine the theoretical development of the finite element method with an approach for determining the dynamic structural response into a seamless computer application that can be expanded for future use. 13 To develop the coupled finite element/structural dynamics program, an understanding of the theoretical development of both topics is necessary. Finite elements are the basis for modeling the structural system into a set of discretized pieces that can be assembled into a set of structural equations. The structural dynamics algorithm can then analyze these structural equations at discrete points in time due to a defined dynamic loading. The combination of these two procedures produces displacements of the structural system as a function of time. In the discussion to follow herein, the theoretical development will center on these topics. Included in these discussions will be formulations of element matrices, derivations of finite element method equations, dynamic analysis procedures and solution techniques which were used while developing the coupled finite element analysis/structural dynamics program. 3.2 Interpolation and Shape Functions The finite element method, at its most basic form, is a set of interpolations or approximations of a dependent variable with respect to an independent variable. Unknown degreesoffreedom (DOF) are utilized to ensure a singlevalued approximation for a set of conditions. An interpolating polynomial with dependent variable and independent variable x in terms of generalized DOF can be expressed in the form Σ or (31) where and can be written as vectors 14 1 … and … (32) For linear interpolation n can be taken 1 and for quadratic interpolation n can be taken 2. The can be written in terms of nodal values (known locations of x) of for known values of x. The relation between and the nodal values can be written as Φ (33) where each row of is calculated at the appropriate nodal location. From the equations above, substitution produces Φ where ! … (34) In matrix an individual is represented as a shape function. In this study, the bilinear rectangle (Q4) finite element is considered. Therefore, linear interpolation is considered when generating the shape functions. That is to say that between two points " , $ and " , $ on a linear line for 1 we can obtain % & ' ( where ) 1 1 * (35) Inverting and using equation 34, ! +,!+ . /1 1 0 and 1 +,!+ +,!+ +!+ +,!+ 2 (36) Equation 36 gives the shape functions of two points on a straight line which are called N1 and N2. 15 Consider the general Q4 element shown in Figure 31. To find the shape functions of the Q4 element, linear interpolation is made along the top and bottom sides to obtain side displacement u12 and u43. Thus, in equation 35, / and , so that 3 4!+ 4 3 5 46+ 4 , 378 4!+ 4 37 5 46+ 4 38 (37) Linear interpolation is then made in the y direction between u12 and u43 as 3 9!: 9 3 5 96: 9 378 (38) Substitution of equations 37 into 38 yields 3 Σ 3 which gives the shape functions of the rectangular Q4 element used in this study as 749 "1 / $"1 / ;$ 749 "1 5 $"1 / ;$ 8 749 "1 5 $"1 5 ;$ 7 749 "1 / $"1 5 ;$ (39) y x a a b b v1 v4 v3 v2 u4 u3 u1 u2 4 1 3 2 Figure 31: Q4 Finite Element Nodal Degrees of Freedom 16 A similar procedure is repeated for the DOF in the y direction. Namely, these DOF are v1, v2, v3 and v4. Similarly, the same shape functions given in equation 39 are found for the shape functions in the ydirection. Therefore, the same interpolation, or approximation, is made in both the x (or u) and y (or v) directions. 3.3 Formulation of Element Stiffness Matrices The principal of virtual work is used to obtain element stiffness matrix for the Q4 element shown in Figure 31. This is appropriate for commonly used elements, which are based on interpolation of displacements from nodal DOF. The principal of virtual work states < => ? @A < =3 B @A 5 < =3 Φ @C (310) where => , B and Φ represent the virtual strains produced by the virtual displacements, body forces and surface tractions, respectively. The displacements {u} are interpolated over an element utilizing shape functions such as those provided by equation 39 as 3 @ where 3 3 D E (311) > = 3 and > F @ where F G (312) where [B] is the straindisplacement matrix. From equations 311 and 312 =3 =@ and => =@ F (313) Substitution of equation 313 back into the statement of virtual work, equation 310, produces 17 =@ " < F H F @A / < F H I @A 5 < F ? @A / JB@A/ JΦ@C 0 (314) Equation 314 can be simplified to produce L @ M (315) As a result of the principle of virtual work, the element stiffness matrix can be determined as L < F H F @A (316) Specifically, for the Q4 element, equation 316 can be written as L < < F H F N@ @; 4 !4 9 !9 (317) For the case of twodimensional plane stress analysis (as considered in this study), the material constitutive matrix [E] is H O " !P,$ Q 1 D 0 D 1 0 0 0 " !P$ R (318) In twodimensional analyses, the thickness, t, in equation 317 is commonly taken as unity. 3.4 Element Boundary and Loading Conditions Before the solution of the structural equations, both boundary conditions and loading conditions must be prescribed for the system. Without boundary conditions the structural equations will not produce a single unique solution for the prescribed loading conditions. As such, the structural system will have rigid body motions. Without loading conditions 18 the structural equations will produce no displacements of the structural system. Thus, it is necessary that both boundary and loading conditions be prescribed for the structural system. These conditions are prescribed at particular DOF of the structural system. Boundary conditions, or support conditions, can be arranged by providing the appropriate stiffness to the related DOF to produce a prescribed displacement. For zero displacement, the prescribed stiffness can be a numerically large number (several orders of magnitude larger than the largest magnitude in the stiffness matrix) such that a relatively small displacement at that DOF is produced. The boundary conditions can be applied to any DOF on the structure no matter its direction. In the present study, only translational DOF are considered in the finite elements (i.e., there are no rotational DOF) and therefore only translational displacements are restricted with respect to particular support conditions such as a fixedsupport, pinnedsupport, or rollersupport. Loading conditions are prescribed in a fashion similar to boundary/support conditions. The load can be applied to the any DOF. The structural system may be subjected to a single loading condition or multiple loading conditions. As this study considers dynamic loads, beyond the typical static load cases considered in a typical finite element analysis, further explanation of the dynamic load is presented in Section 3.6.3. 19 3.5 Solution Techniques Equation 317 is integrated over a rectangular surface as given in Figure 31. Numerically, however, to integrate this equation is cumbersome and not efficient. Thus, integration is achieved using Gauss Quadrature. Gauss Quadrature is a numerical approximation of the integration by use of simple algebraic equations evaluated at specific points. To use Gauss Quadrature, however, the element must be formulated in the isoparametric space. The use of the isoparametric space specifies that an isoparametric element is used rather than the physical element. Thus, the integrands in the integration formulas are expressed as functions of S and T rather than x and y. For the function ø = ø (ξ, η), the Quadrature rule is given as U < < ø"ξ, η$dξ dη ! Σ ΣZY YZø"ξ, η$ ! (319) where Y and YZ are weighting factors for each Gauss point i and j. The weighting factors for twopoint Gauss Quadrature, as used in this study for the Q4 element, are taken as unity. For the Q4 element shown in Figure 31, the individual isoparametric shape functions can simply be determined from equation 39 by assigning a = 1, b = 1, x = ξ and y = η producing 7 "1 / S$"1 / T$ 7 "1 5 S$"1 / T$ 8 7 "1 5 S$"1 5 T$ 7 7 "1 / S$"1 5 T$ (320) 20 Since the isoparametric Q4 element is not in the same coordinate system as the physical Q4 element, a mapping function is needed to relate the two coordinate systems. The Jacobian matrix is used to accomplish this mapping. The Jacobian matrix is a scale factor that multiplies dξdη to produce the physical area increment dxdy and is expressed as (for the Q4 isoparametric element) [ 7 \ – "1 / T$ "1 / T$ "1 5 T$ – "1 5 T$ – "1 / S$ – "1 5 S$ "1 5 S$ "1 / S$ ^ _ ; ; 8 ;8 7 ;7 ` ) [ [ [ [ *(321) Thus, the element stiffness matrix for the Q4 isoparametric element is then L a F H F N @ @; < < F H F N [ @S @T ! ! (322) By substituting equation 322 into equation 319, the integration can be numerically calculated using Gauss Quadrature. Following the determination of the stiffness of each element in the finite element mesh, the system stiffness must be assembled. Based on the arrangement of the finite element mesh, the stiffness for each element corresponding to particular system DOF is “fed” into the system stiffness corresponding to the same system DOF. This process, called the assembly process, is a simply a mapping technique relating the DOF corresponding to each element to the DOF of the system. After assembly of the elements into the system equations, the structural equations are then in the form that can be solved. Gauss elimination is used to solve the structural equations for given boundary and loading conditions. In Gauss elimination, equations 21 [K]{D} = {R} are solved for {D}, the displacements, by reducing [K] to upper triangular form and then solving for unknowns in the reverse order by back substitution. 3.6 Structural Dynamics An integral part of this study is the incorporation of a structural dynamics response algorithm into the finite element method. These two coupled methods will produce a seamless method by which to analyze a structural system, using the finite element method, under an applied dynamic load(s) for a specified period of time. As a result, the response of the structural system as a function of time will be determined. The following sections explain the theoretical development of the structural dynamics methodology used in this study. Structural dynamics derivations and solution techniques are outlined briefly. The equations are derived for timedependent loads. This study only considers undamped structures, thus, the formulation of damping is not presented and it is taken as zero in all equations. The following sections introduce the theoretical explanations of structural dynamics that were used in the development of the coupled computer program. 3.6.1 Equation of Motion The equation of motion is the basic and fundamental part of structural dynamics. All formulations are derived based on the equation of motion in structural dynamics. The equation of motion is generally given as bü 5 d3e 5 L3 f"N$ (323) 22 The right side of the equation is the time dependent force, p(t). On the left side of equation 323, m represents mass, c represents damping, k represents stiffness and u, 3e, ü represent displacement, velocity and acceleration, respectively (the overdot represents a time derivative). The structure system, which has multiple DOFs, is of course considered in this study. The specific number of DOF depends on the finite element model of the structural system. As a result, all variables in equation 323 are matrices or column vectors of size related to the number of DOF for the system. 3.6.2 Formulation of Element Mass Matrices The formulation of the element mass matrix is based on the virtual work principal and is similar to the formulation of the element stiffness matrix as discussed in Section 3.3. The work done by externally applied loads is equal to the sum of the work absorbed by inertial, dissipative, and internal forces for any virtual displacement. For an element volume V and surface S < =3 B @A 5 < =3 Φ @C 5 Σ =3 f <" =3 g ü 5 =3 d 3e 5 => ? $@A (324) Where B and Φ are prescribed body forces and surface tractions, f and =3 are prescribed concentrated loads and their corresponding virtual displacements, g is mass density and c is the damping coefficient. Following Section 3.3, the nodal displacements, nodal velocities, nodal accelerations, and strains are approximated by 3 @ 3e h@ei ü h@ji > F @ (325) Substitution of equation 325 into equation 324 produces the virtual work expression 23 =@ k< g @Ah@ji 5 < d @Ah@ei 5 < F ? @A / JB@A/ JΦ@C/ l 1mfl 0 (326) The first integral in equation 326 provides the element mass matrix as b < g @A (327) For twodimensional analysis, equation 327 yields b < < g N@ @; 4 !4 9 !9 (328) Following the procedure for assembly of element stiffness matrix into the system matrix, the element mass matrix for each individual element is assembled in the same procedure. The DOF for the mass and stiffness matrices are identical. Thus, from a numerical perspective, the element stiffness and mass matrices can be assembled into the system stiffness and mass matrices simultaneously. 3.6.3 Dynamic Loading Conditions There are many types of dynamic loads from a realistic perspective. In this study, several common idealized dynamic loadings are considered. Namely, these include the harmonic loading, the impulse loading and the step loading. It should be noted, however, that the MATLAB program written as a consequence of this study can handle any userinput for the dynamic loading. The user only has to provide the value of the load at specified time intervals, or sampling points. These three idealized loads are only chosen based on their commonality in structural dynamic simplifications. 24 Harmonic loading is a function of the sine or cosine functions and its equation is in general f"N$ f nlmoN or f drnoN (329) where p0 is the amplitude or maximum value of the force and its frequency o is the exciting frequency or forcing frequency. Figure 32 provides an example of a harmonic loading. Examples of a harmonic loading includes wind driven loading on structures, earthquake (highly idealized) loading, and vehicular motion on bridges. Figure 32: Harmonic Loading A very large force that acts for a very short time but with a time integral that is finite is called an impulsive force. In general, an impulsive force is defined by f"N$ 1/> (330) with a time duration > starting at the time instant t = t, also called the time lag. Figure 3 3 is an illustration of an impulse loading. Examples of an impulsive load primarily include blast loadings such as those due to detonation of blast devices. 25 Figure 33: Impulse Loading Finally, another typical dynamic loading is a force that jumps suddenly from zero to magnitude p0 and stays constant at that value is called step force. In general, the step loading is defined by f"N$ f (331) Figure 34 provides an example of an impulse loading that jumps to its magnitude p0 at time 0 seconds as indicated. Examples of a step loading include the sudden application of a full load rather than being applied gradually over time. Figure 34: Step Loading 26 3.6.4 Numerical Integration There are several numerical integration methods for linear systems to numerically integrate the equation of motion previously defined as equation 323. In all methods for the purposes of this study, initial displacements and velocities are taken zero (the system is initially at rest in an underformed position) and p(t) is known at all time intervals t. The aim of numerical integration is analyzing the system over time intervals ΔN. Some of the numerical integration methods, that are called timestepping methods, are interpolation of excitation, central difference method and Newmark’s method. Interpolation of excitation uses recurrence formulas. It can be used for small ΔN and for linear systems. In addition, this method is suitable for single DOF systems but it is not appropriate for multi DOF systems. Central difference methods are based on finite difference approximations of the time derivatives of displacement, velocity and acceleration. Solution at 3 6 which represents the displacement at the step i+1 is determined from the equation of motion at time step i. Furthermore, 3 and 3 ! must be known to find the displacement at time step i+1, 3 6 . Newmark’s method, used in this study, is a family of timestepping methods based on the iteration of ü 6 . A detailed discussion is presented in the next section. 27 3.6.5 NewmarkBeta Solution Technique Newmark’s technique, also known as the more general NewmarkBeta solution technique, is based on the following equations 3e 6 3e 5 "1 / v$ΔN ü 5 "vΔN$ü 6 (332) 3 6 3 5 "ΔN$3e 5 "0.5 / z$"ΔN$ ü 5 z"ΔN$ ü 6 (333) The parameters z and v define the variation of acceleration over a time step and determine the stability and accuracy characteristics of the technique. Typical selection of v is 1/2, and z can vary between 1/4 and 1/6. For the average acceleration method, used in this study, z is taken as 1/4. All equations are matrix equations for multi DOF as is the case in this study. The NewmarkBeta solution procedure can be presented step by step for linear, multi DOF systems. Initial conditions are defined as 3 3"0$ and 3e 3e"0$ representing the initial displacement and initial velocity, respectively. Initial calculations are only calculated one time and they are presented in equations 334 thorough 337 as Solving for the initial acceleration ü ; bü f / d3e / L3 (334) After selectingΔN and the NewmarkBeta parameters z and v the effective stiffness of the system can be calculated as L{ L 5  }Δ~ d 5 }Δ~,b (335) Calculation constants, a and b, for use at each time step are calculated as 28 }Δ~b 5  } d (336) }b 5 ΔN"  } / 1$d (337) For each time step, equations 338 thorough 342 are be repeated until all time steps are done. The effective force at time step i is then calculated as Δf̂ Δf 5 3e 5 ü (338) Solving for the change in displacement during the time step Δ3 is then L{ Δ3 Δf̂ (339) During the same time step, the change in velocity and the change in acceleration are also Δ3e  }Δ~ Δ3 /  } 3e 5 ΔN 1 /  } ü (340) Δü }Δ~, Δ3 / }Δ~ 3e / } ü (341) Finally, updating the displacement, velocity and acceleration can be found relative to the previous position, velocity and acceleration as 3 6 3 5 Δ3 3e 6 3e 5 Δ3e (342) ü 6 ü 5 Δü Repeating these calculations, equations 338 through 342, the displacements, velocities and accelerations can be found for discrete time points for the structural system. 29 CHAPTER IV CASE STUDIES 4.1 Introduction Case studies are provided to illustrate the capabilities of the MATLAB FEA program. This code has one main program file, one input file and multiple subroutines to analyze the structural system. Material properties, geometrical properties, boundary conditions, loading conditions and more can be easily defined by the user in the input file of the MATLAB code. The case studies all consider twodimensional dynamically loaded beams modeled using a finite element mesh of Q4 elements. To confirm the procedure and accuracy, a basic example is solved by hand in Appendix B of this study. Additionally, the main program and input file of the MATLAB code is provided in Appendix A of this study. Table 41 shows the beams properties that are used in the four case studies. E represents modulus of elasticity; and g represent Poisson's ratio and the mass density of the beam, respectively. In the first case study, the beam is fixed at the left end. The second case study considers a beam that is simply supported. Finally, the third and fourth case studies both consider continuous beams that are supported at the left end, the right end, and at the 30 midpoint. The geometry of the beams, boundary conditions and loading conditions are provided with the case studies. Modulus of Elasticity (psi) Poisson’s Ratio Mass Density (lbs/in3) Number of Elements Number of Nodes Type of Loading Beam Length (in) Beam Depth (in) Case 1 4.35(106) 0.15 0.0868 48 65 Harmonic 144 12 Case 2 4.35(106) 0.15 0.0868 48 65 Impulse 144 12 Case 3 4.35(106) 0.15 0.0868 57 80 Step 240 12 Case 4 4.35(106) 0.15 0.0868 57 80 Multiple 240 12 Table 41: Case Studies As seen in Table 41, four types of loading cases are considered. The system is undamped for all cases. Because twodimensional analyses are considered, each node has two degrees of freedom in the x and y directions while the beam width is taken as unity for all cases. The dynamic analysis duration is 10 seconds, the time step t is taken as 0.1 seconds and initial displacements and velocities for all DOF are zero at time = 0 seconds. With the given time duration and time step, each case study has 100 time steps to solve in the dynamic response analysis. The MATLAB FEA program has the capability to produce movies of the displacements during the time duration of the dynamic loading. As it is not possible to place a movie in a text format, select “snapshots” have been provided for each case study to illustrate the dynamic response of the structural system. These figures serve only as illustrations as to the full capabilities of the MATLAB FEA program developed as a consequence of this study. 31 4.2 Harmonic Loading A harmonic loading is applied to the first beam which is cantilever beam. The illustration of the beam is given in the Figure 41. The applied load F represents the harmonic load as it is applied at the right side of the beam. The properties of the beam are provided in Table 41 as Case 1. Figure 41: Case Study 1 – Cantilever Beam The harmonic load is a sine function and is illustrated in Figure 42. The harmonic load is applied between 0 and 10 seconds with sampling points taken at intervals of 0.1 seconds. The magnitude of the load varies between 1800 lbs and 1800 lbs. 32 Figure 42: Case Study 1 – Harmonic Loading As stated in Table 41, the beam is formulated with 48 rectangular Q4 elements having 130 DOFs. Illustration of the finite element mesh of the beam is shown in Figure 43. Green circles indicate node numbering and red circles indicate element numbering in the mesh. Figure 44 illustrates the displacement of the beam considering a mesh of Q4 finite elements. Subfigures (a) – (f) provide a representation of the movie generated at discrete time intervals as indicated on each subfigure. 0 1 2 3 4 5 6 7 8 9 10 1500 1000 500 0 500 1000 1500 Time [seconds] Load [lbs] Harmonic Load 33 Figure 43: Case Study 1 – Finite Element Mesh Figure 44(a): Case Study 1 – Displacements 0 20 40 60 80 100 120 140 2 0 2 4 6 8 10 12 14 1 1 2 14 15 2 3 16 3 4 17 4 5 18 5 6 19 6 7 20 7 8 21 8 9 22 9 10 23 10 11 24 11 12 25 12 13 26 13 14 15 16 17 18 19 20 21 22 23 24 25 27 28 40 41 26 29 42 27 30 43 28 31 44 29 32 45 30 33 46 31 34 47 32 35 48 33 36 49 34 37 50 35 38 51 36 39 52 37 53 54 38 55 39 56 40 57 41 58 42 59 43 60 44 61 45 62 46 63 47 64 48 65 length [in] depth [in] Finite Element Mesh 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1 seconds 34 Figure 44(b): Case Study 1 – Displacements Figure 44(c): Case Study 1 – Displacements 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.7 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =3.5 seconds 35 Figure 44(d): Case Study 1 – Displacements Figure 44(e): Case Study 1 – Displacements 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =5 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =8.1 seconds 36 Figure 44(f): Case Study 1 – Displacements As seen from the subfigures of Figure 44, the free end of the beam moves up and down due to the harmonic loading. Comparison of hand calculations of the steadystate displacement response and computer results are presented in Appendix B. DOF 78 at node 39 is considered while comparing the results of displacements in Appendix B. 4.3 Impulse Loading An impulse loading is applied to the 12ft long simplysupported beam. The load is applied at the middle of the beam. The illustration of the beam is shown in Figure 45. The applied load F represents the impulse load in the figure. The properties of the beam are provided in Table 41 as Case 2. 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =10 seconds 37 F Length = 12' Depth = 1' 6' 6' Figure 45: Case Study 2 – SimplySupported Beam The impulse loading, as shown in Figure 46, is applied between 0 and 10 seconds with sampling points taken at intervals of 0.1 seconds. At time = 1s the impulse loading immediately jumps to a value of –2500 lbs. At all other times t, the loading is zero. Figure 46: Case Study 2 – Impulse Loading 0 1 2 3 4 5 6 7 8 9 10 2.5 2 1.5 1 0.5 0 x 10 4 Time [seconds] Load [lbs] Impulse Load 38 As stated in Table 41, the beam, which has 130 DOFs, is constructed with 48 rectangular Q4 elements. Illustration of the finite element mesh of the beam is shown in Figure 47. Figure 48 illustrates the displacement of the beam considering a mesh of Q4 finite elements. Subfigures (a) – (f) provide a representation of the movie generated at discrete time intervals as indicated on each subfigure. Figure 47: Case Study 2 – Finite Element Mesh 0 20 40 60 80 100 120 140 2 0 2 4 6 8 10 12 14 1 1 2 14 15 2 3 16 3 4 17 4 5 18 5 6 19 6 7 20 7 8 21 8 9 22 9 10 23 10 11 24 11 12 25 12 13 26 13 27 28 14 29 15 30 16 31 17 32 18 33 19 34 20 35 21 36 22 37 23 38 24 39 25 26 27 28 29 30 31 32 33 34 35 36 37 40 41 53 54 38 42 55 39 43 56 40 44 57 41 45 58 42 46 59 43 47 60 44 48 61 45 49 62 46 50 63 47 51 64 48 52 65 length [in] depth [in] Finite Element Mesh 39 Figure 48(a): Case Study 2 – Displacements Figure 48(b): Case Study 2 – Displacements 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.1 seconds 40 Figure 48 (c): Case Study 2 – Displacements Figure 48(d): Case Study 2 – Displacements 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.3 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.4 seconds 41 Figure 48(e): Case Study 2 – Displacements Figure 48(f): Case Study 2 – Displacements Until time = 1 second there is no displacement at the beam as expected. At that time, the impulse load is applied to the beam. Because the system is undamped, after 1 second the 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =5 seconds 0 20 40 60 80 100 120 140 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =10 seconds 42 beam starts oscillating up and down and goes forever with the maximum displacement at the midpoint as seen from the Figure 48 and its subfigures. 4.4 Step Loading A step loading is applied to the third case study. The 20ft long beam has a pin support at the left end and roller supports at the middle and right end. The load is located 5 feet from the right side of the beam. The illustration of the beam is shown in Figure 49. The applied load F represents the step load in the figure. The properties of the beam are provided in Table 41 as Case 3. 10' F 5' 5' Length = 20' Figure 49: Case Study 3 – Continuous Beam The step loading, as shown in Figure 410, is applied between 0 and 10 seconds with sampling points taken at intervals of 0.1 seconds. As with the impulse loading of Case 2, 43 there is no load applied before the time reaches 1 second. At this point in time, the value of the load jumps suddenly to –20000 lbs and stays constant for the duration of time. Figure 410: Case Study 3 – Step Loading As stated in Table 41, the beam is comprised of 57 rectangular Q4 elements having 160 DOFs. Illustration of the finite element mesh of the beam is shown in Figure 411. Figure 412 illustrates the displacement of the beam considering a mesh of Q4 finite elements. Subfigures (a) – (f) provide a representation of the movie generated at discrete time intervals as indicated on each subfigure. 0 1 2 3 4 5 6 7 8 9 10 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 x 10 4 Time [seconds] Load [lbs] Step Load 44 Figure 411: Case Study 3 – Finite Element Mesh Figure 412(a): Case Study 3 – Displacements 0 50 100 150 200 2 0 2 4 6 8 10 12 14 1 1 2 21 22 2 3 23 3 4 24 4 5 25 5 6 26 6 7 27 7 8 28 8 9 29 9 10 30 10 11 31 11 12 32 12 13 33 13 14 34 14 15 35 15 16 36 16 17 37 17 18 38 18 19 39 19 20 40 20 41 42 21 43 22 44 23 45 24 46 25 47 26 48 27 49 28 50 29 51 30 52 31 53 32 54 33 55 34 56 35 57 36 58 37 59 38 60 39 61 62 40 63 41 64 42 65 43 66 44 67 45 68 46 69 47 70 48 71 49 72 50 73 51 74 52 75 53 76 54 77 55 78 56 79 57 80 length [in] depth [in] Finite Element Mesh 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1 seconds 45 Figure 412(b): Case Study 3 – Displacements Figure 412 (c): Case Study 3 – Displacements 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1.2 seconds 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =3.5 seconds 46 Figure 412(d): Case Study 3 – Displacements Figure 412(e): Case Study 3 – Displacements 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =5.2 seconds 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =7 seconds 47 Figure 412(f): Case Study 3 – Displacements Because of the step load, the beam displacement never becomes a positive displacement (with respect to the coordinate system used) between the roller supports and never becomes a negative displacement between the pin and roller supports. These results are as expected. 4.5 Multiple Dynamic Loads The final case study considers the same beam that was used for the third case study when the step loading was considered. However, the loading condition is different in this case. F1 represents a harmonic load and F2 represents a step load. The harmonic load is located 5ft away from the pin support and the step load is located 5ft away from the left end of the beam as shown in the Figure 413. 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =10 seconds 48 F2 5' 5' Length = 20' F1 5' 5' Figure 413: Case Study 4 – Continuous Beam The harmonic load is applied to the beam from time = 0 seconds until 10 seconds and its magnitude varies between 15000 lbs to 15000 lbs. At time = 4 seconds, suddenly another load is applied, the step load, and it stays constant for the duration of time with a magnitude of –20000 lbs. These two timedependent loads are shown the Figure 414. 49 Figure 414: Case Study 4 – Multiple Loading As stated in Table 41, the beam is comprised of 57 rectangular Q4 elements and having 160 DOFs. The finite element mesh of the beam is illustrated in the Figure 415. Figure 416 illustrates the displacement of the beam considering a mesh of Q4 finite elements. Subfigures (a) – (f) provide a representation of the movie generated at discrete time intervals as indicated on each subfigure. 0 1 2 3 4 5 6 7 8 9 10 2 1.5 1 0.5 0 0.5 1 x 10 4 Time [seconds] Load [lbs] Multiple Load harmonic load step load 50 Figure 415: Case Study 4 – Finite Element Mesh Figure 416(a): Case Study 4 – Displacements 0 50 100 150 200 2 0 2 4 6 8 10 12 14 1 1 2 21 22 2 3 23 3 4 24 4 5 25 5 6 26 6 7 27 7 8 28 8 9 29 9 10 30 10 11 31 11 12 32 12 13 33 13 14 34 14 15 35 15 16 36 16 17 37 17 18 38 18 19 39 19 20 40 20 41 42 21 43 22 44 23 45 24 46 25 47 26 48 27 49 28 50 29 51 30 52 31 53 32 54 33 55 34 56 35 57 36 58 37 59 38 60 39 61 62 40 63 41 64 42 65 43 66 44 67 45 68 46 69 47 70 48 71 49 72 50 73 51 74 52 75 53 76 54 77 55 78 56 79 57 80 length [in] depth [in] Finite Element Mesh 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =1 seconds 51 Figure 416(b): Case Study 4 – Displacements Figure 416(c): Case Study 4 – Displacements 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =3 seconds 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =3.8 seconds 52 Figure 416(d): Case Study 4 – Displacements Figure 416(e): Case Study 4 – Displacements 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =4.1 seconds 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =5 seconds 53 Figure 416(f): Case Study 4 – Displacements As shown in Figure 416, the influence of the step load when applied at time = 4 seconds is readily distinguishable. Before the application of the step load, displacements are limited due to the harmonic load. However, clearly the displacement of the beam before 4 seconds is harmonic in nature as illustrated. After a period of time after 4 seconds, the displacement again begins to resemble a harmonic pattern as expected. 0 50 100 150 200 4 2 0 2 4 6 8 10 12 14 16 length [in] depth [in] Dynamic Displacements / Time =9 seconds 54 CHAPTER V CONCLUSION 5.1 Summary and Conclusions The primary objective of this study was to develop a MATLAB based computer program coupling the finite element method and structural dynamics analysis techniques. The primary reason for using the finite element method is its significant advantages in engineering problems. It is a powerful solution technique for differential and integral equations in complex engineering problems. When combined with a structural dynamic solution algorithm, the capabilities of the two become very unique. Theoretical equations are derived for both the finite element method and structural dynamics during this study. Fundamentally, formulation of element matrices, integration techniques, boundary conditions, and timestepping methods are discussed and their equations are presented. These concepts are used to develop the computer code. Several case studies are analyzed during this study with different types of dynamic loads including harmonic loading, impulse loading, and step loading. These loads are also applied to a variety of different beams including a cantilever beam, a simplysupported 55 beam and a continuous beam with three supports. Displacements and deformed shapes are found using MATLAB FEA program written as a consequence of this study. This study presents a discussion of the finite element method for structural dynamics applications. It provides a basic understanding of the behavior of structural systems under different dynamic loads. The computer code, which is written in MATLAB, allows for further flexibility not explicitly discussed in this study such as using different material properties throughout the structure, a variety of geometrical shapes of the structure, varying types of loads, and multiple finite element types beyond the Q4 element used in the finite element mesh. 5.2 Recommendations This study is open to development and further development is expected in the future. Although currently it is limited to Q4 elements for beams, it can be expanded to advanced finite elements and other types of structures (see Iqbal, 2009). The MATLAB code is written to be expanded to allow for the analysis of more complicated structural systems utilized more advanced finite elements and finite element formulations. Moreover, one can develop this study in not only finding displacements, deformed shapes, stresses, and strains but also other topics such as calculating fracture, durability or even reliability of structures. 56 Although the computer code has good properties, it is also recommended that better meshing capabilities may be implemented in the MATLAB code. The computer code covers linear programming for now. Therefore, the computer program can be incorporated with nonlinearity. If these modifications can be done, this computer code can become a more powerful finite element analysis program. Work done in conjunction with this study includes the incorporation of advanced finite element formulations beyond the Q4 element used herein (Iqbal, 2009). 57 REFERENCES [1] Bozdogan, K.B., Ozturk, D., and Nuhoglu, A. (2009), “An Approximate Method for Static and Dynamic Analyses of Multibay Coupled Shear Walls”, The Structural Design of Tall and Special Buildings, Vol. 18, pp. 112. [2] Cerioni, R., Brighenti, R., and Donida, G. (1995), “Use of Incompatible Displacement Modes in a Finite Element Model to Analyze the Dynamic Behavior of Unreinforced Masonry Panels”, Computers and Structures, Vol. 57, pp. 47  57. [3] Chopra, K.A. (2007), “Dynamics of Structures Theory and Applications to Earthquake Engineering ”, Pearson Prentice Hall. [4] Cook, R.D., Malkus, D.S., Plesha, M.E. and Witt, R.J. (2002), “Concepts and Applications of Finite Element Analysis”, John Wiley & Sons. [5] Du, H., Hitchings, D., and Davies, G.A.O. (1992), “A Finite Element Structural Dynamics Model of a Beam with an Arbitrary Moving Base – Part I: Formulations”, Finite Elements in Analysis and Design, Vol. 12, pp. 133  150. [6] Kiral, Z., and Kiral, B.G. (2008). “Dynamic Analysis of a Symmetric Laminated Composite Beam Subjected to a Moving Load with Constant Velocity”, Journal of Reinforced Plastics and Composites, Vol. 27, pp. 1932. 58 [7] Li, Q.S., Yang, L.F., and Li, G.Q. (2001), “The Quadratic Finite Element/Strip with Generalized Degrees of Freedom and Their Application”, Finite Elements in Analysis and Design, Vol. 37, pp. 325  339. [8] Mermertas, V. (1998), “Dynamic Interaction Between the Vehicle and Simply Supported Curved Bridge Deck”, Computer Methods in Applied Mechanics and Engineering, Vol. 162, pp. 125  131. [9] Ozcan, D.M., Bayraktar, A., Sahin, A., Haktanir, T., and Turker, T. (2008), “Experimental and Finite Element Analysis on the Steel FiberReinforced Concrete (SFRC) Beams Ultimate Behavior. Construction and Building Materials, Vol. 23, pp. 1064  1077. [10] Iqbal, R. (2009), “Development of a Finite Element Program Incorporating Advanced Element Types”, MS Thesis, Oklahoma State University. 59 APPENDIX A MATLAB FEA PROGRAM The Matlab FEA program that was developed during the course of this study is provided herein. First, the main FEA program (not containing the subroutines since the FEA program is still under development) is provided. Second, the general input file used to setup an analysis is provided. Both of these files are Matlab script files and can be executed in all versions of Matlab. To generate graphical output, the full version of Matlab is required. Main Program % ************************************************************************* % Finite Element Analysis (FEA) Program to Determine the Structural and Thermal Response of Structural Systems % Written By: Muhammet Saglar & Rameez Iqbal (Advised By: Dr. Jonathan S. Goode) % School of Civil and Environmental Engineering % Oklahoma State University % ************************************************************************* % ************************************************************************* % ************************************************************************* % Begin Program % Clear Variables and Settings clear all % Clear Screen clc; % Inhibit Warning Messages warning('off','all'); disp('**************************************************'); disp('Finite Element Analysis') disp('Written by: Saglar/Iqbal (Goode)'); disp('School of Civil and Environmental Engineering'); disp('Oklahoma State University'); disp('**************************************************'); 60 disp('**************************************************'); disp('Running Analysis .....'); disp('**************************************************'); % Begin Timer tic; % ************************************************************************* % PreProcessing Section % ************************************************************************* disp('PreProcessing .....'); % Call MFile for User Input Data from fea_input.m disp('..... Obtaining Input .....'); fea_input; % Call MFile for Initialization of Variables from fea_initialize.m disp('..... Initializing Program .....'); fea_initialize; % ************************************************************************* % End PreProcessing Section % ************************************************************************* % ************************************************************************* % Processing Section % ************************************************************************* disp('Processing .....'); %  % Assembly Process disp('..... Assembly Process .....'); % Determine Element Property Matrices (Constitutive, Mass, Conducitivity) [property_matrix_k,property_matrix_m,property_matrix_t] = fea_property_matrix(analysis_type_material,material_E,material_v,material_rho,material_t); % Determine Gauss Quadrature Sampling Points and Weights [gauss_points,gauss_weights] = fea_gauss_2d(element_type); % Initialize Global Index of DOF index_global = zeros(edof,nel); % Assemble System Matrices According to Analysis if analysis_type == 'S' % Static Analysis disp('..... ..... Static Analysis .....'); disp('..... ..... ..... Building Stiffness .....'); % Loop Through All Elements to Determine System Matrices (Stiffness, Mass, and Conductivity) % NOTE: Depending on the Analysis Type, Some Matrices Will Not be Generated for j = 1:nel % Call SubRoutine to Extract Nodes for the jth Element [element_nodes,element_xcoord,element_ycoord] = fea_node_coordinates(j,nnel,node_connectivity,node_coord); % Call SubRoutine to Extract System DOFs for the jth Element index = fea_index_dof(element_nodes,nnel,ndof,edof); % Save Index of DOF to Global Variable index_global(:,j) = index; % Call SubRoutine to Generate Element Stiffness Matrix 61 element_k = fea_element_matrix_stiffness(element_type,gauss_points,gauss_weights,element_xcoord,element_ycoord,nnel,edof,property_matrix_ k); % Call SubRoutine to Assemble System Stiffness Matrix system_k = fea_element_matrix_assemble(system_k,element_k,index,edof); end elseif analysis_type == 'D' % Dynamic Analysis disp('..... ..... Dynamic Analysis .....'); disp('..... ..... ..... Building Stiffness and Mass .....'); % Loop Through All Elements to Determine System Matrices (Stiffness, Mass, and Conductivity) % NOTE: Depending on the Analysis Type, Some Matrices Will Not be Generated for j = 1:nel % Call SubRoutine to Extract Nodes for the jth Element [element_nodes,element_xcoord,element_ycoord] = fea_node_coordinates(j,nnel,node_connectivity,node_coord); % Call SubRoutine to Extract System DOFs for the jth Element index = fea_index_dof(element_nodes,nnel,ndof,edof); % Save Index of DOF to Global Variable index_global(:,j) = index; % Call SubRoutine to Generate Element Stiffness Matrix element_k = fea_element_matrix_stiffness(element_type,gauss_points,gauss_weights,element_xcoord,element_ycoord,nnel,edof,property_matrix_ k); % Call SubRoutine to Assemble System Stiffness Matrix system_k = fea_element_matrix_assemble(system_k,element_k,index,edof); % Call SubRoutine to Generate Element Mass Matrix element_m = fea_element_matrix_mass(element_type,gauss_points,gauss_weights,element_xcoord,element_ycoord,nnel,edof,property_matrix_m); % Call SubRoutine to Assemble System Mass Matrix system_m = fea_element_matrix_assemble(system_m,element_m,index,edof); end end % End Assembly Process %  %  % Boundary and Loading Conditions disp('..... Applying Boundary and Loading Conditions .....'); % Apply Boundary and Loading Conditions According to Analysis if analysis_type == 'S' % Static Analysis % Call SubRoutine to Apply Boundary and Loading Conditions [system_k,system_f] = fea_system_boundary_load(system_k,system_f,boundary_dof,boundary_dof_val,load_dof,load_dof_val); elseif analysis_type == 'D' % Dynamic Analysis % Boundary and Loading Conditions Must Be Applied at Each Point in Time During the Dynamic Analysis % See the Analysis Section Below 62 end % End Boundary Conditions %  %  % Analysis disp('..... Solving .....'); % Based on Analysis Type if analysis_type == 'S' % Static Analysis disp('..... ..... Static Analysis .....'); % Solve System Equations for Static Displacements system_d = system_k \ system_f; elseif analysis_type == 'D' % Dynamic Analysis disp('..... ..... Dynamic Analysis .....'); % Call SubRoutine to Solve System Equations for Dynamic Displacements [system_d,system_d_vel,system_d_acc,system_f,dynamic_time] = fea_dynamic_analysis(system_k,system_m,boundary_dof,boundary_dof_val,load_dof,load_file,dynamic_delta_t,dynamic_beta,dyna mic_gamma,sdof); end % End Analysis %  % ************************************************************************* % End Processing Section % ************************************************************************* % ************************************************************************* % PostProcessing Section % ************************************************************************* disp('PostProcessing .....'); %  % Plot Generation % Plot Finite Element Mesh if plot_element_mesh == 'Y' disp('..... Plotting Element Mesh .....'); % Call SubRoutine to Plot Finite Element Mesh plot_handle = fea_plot_mesh(proj_name,nel,node_coord,node_connectivity,plot_element_mesh_numbers,plot_element_mesh_node_numbers); end % Based on Analysis Type if analysis_type == 'S' % Static Analysis disp('..... Static Analysis .....'); % Plot Static Displacements if plot_static_displacement == 'Y' 63 disp('..... ..... Plotting Static Displacements .....'); % Call SubRoutine to Plot Static Displacements plot_handle = fea_plot_static_displacement(proj_name,nel,node_coord,node_connectivity,system_d,index_global); if plot_static_contour_displacement == 'Y' % Call SubRoutine to Plot Displacement Contours plot_handle = fea_plot_static_displacement_contour(proj_name,node_coord,node_connectivity,system_d,sdof,nel); end end % Determine Static Stresses and Strains for Plots if plot_static_stress == 'Y'  plot_static_strain == 'Y' disp('..... ..... Determining Static Stresses and Strains .....'); % Call SubRoutine to Determine Stresses and Strains from Displacements [system_stress,system_strain,gauss_points_coord] = fea_stress_strain(nel,nnel,edof,node_connectivity,node_coord,index_global,system_d,property_matrix_k,element_type,gauss_points); end % Plot Static Stresses if plot_static_stress == 'Y' disp('..... ..... Plotting Static Stresses .....'); % Call SubRoutine to Plot Static Stresses plot_handle = fea_plot_static_stress(proj_name,nel,node_coord,node_connectivity,system_stress,gauss_points_coord); if plot_static_contour_stress % Call SubRoutine to Plot Stress Contours [plot_handle] = fea_plot_static_stress_contour(proj_name,nel,node_coord,node_connectivity,system_stress,gauss_points_coord); end end % Plot Static Strains if plot_static_strain == 'Y' disp('..... ..... Plotting Static Strains .....'); % Call SubRoutine to Plot Static Strains plot_handle = fea_plot_static_strain(proj_name,nel,node_coord,node_connectivity,system_strain,gauss_points_coord); if plot_static_contour_strain % Call SubRoutine to Plot Strain Contours [plot_handle] = fea_plot_static_strain_contour(proj_name,nel,node_coord,node_connectivity,system_strain,gauss_points_coord); end end elseif analysis_type == 'D' % Dynamic Analysis disp('..... Dynamic Analysis .....'); % Call SubRoutine to Plot Dynamic Displacements (Movie) 64 plot_handle = fea_plot_dynamic_displacement(proj_name,nel,node_coord,node_connectivity,system_d,dynamic_time,index_global); disp('..... ..... Plotting Dynamic Displacements .....'); end % ************************************************************************* % End PostProcessing Section % ************************************************************************* disp('**************************************************'); disp('Analysis Done .....'); disp('**************************************************'); % End Timer analysis_time = toc; disp(['Total Time of Analysis: ',num2str(analysis_time),' seconds']); disp('**************************************************'); % End Program % ************************************************************************* % ************************************************************************* Input File % ************************************************************************* % Finite Element Analysis (FEA)  Input File % Written By: Muhammet Saglar & Rameez Iqbal (Advised By: Dr. Jonathan S. Goode) % School of Civil and Environmental Engineering % Oklahoma State University % ************************************************************************* % ************************************************************************* % ************************************************************************* % Begin Input File %  % Project Name % Define Project Name % NOTE: Automatically Saved Figures and Data Will Be Proceeded with the Project Name proj_name = 'test'; %  % Analysis Options % Define Analysis Type % Analysis Type Options: % Input 'S' = Static Analysis % Input 'D' = Dynamic Analysis % Input 'T' = Thermal Analysis analysis_type = 'S'; %  % Input Finite Element Mesh Properties % Input Element Type % Input 'Q4' for Bilinear Rectangular Element % Input 'Q8' for Quadratic Rectangular Element % Input 'CST' for Constant Strain / Linear Triangular Element % Input 'LST' for Linear Strain / Quadratic Triangular Element 65 element_type = 'Q4'; % Input Node Coordinates (x and y coordinates) [Note Units => Length] node_coord = [ 0.0 0.0 ;... 12.0 0.0 ;... 24.0 0.0 ;... 36.0 0.0 ;... 48.0 0.0 ;... 60.0 0.0 ;... 72.0 0.0 ;... 84.0 0.0 ;... 96.0 0.0 ;... 108.0 0.0 ;... 120.0 0.0 ;... 0.0 12.0 ;... 12.0 12.0 ;... 24.0 12.0 ;... 36.0 12.0 ;... 48.0 12.0 ;... 60.0 12.0 ;... 72.0 12.0 ;... 84.0 12.0 ;... 96.0 12.0 ;... 108.0 12.0 ;... 120.0 12.0 ]; % Input Nodal Connectivity for Each Element (CCW from BottomLeft) node_connectivity = [ 1 2 13 12 ;... 2 3 14 13 ;... 3 4 15 14 ;... 4 5 16 15 ;... 5 6 17 16 ;... 6 7 18 17 ;... 7 8 19 18 ;... 8 9 20 19 ;... 9 10 21 20 ;... 10 11 22 21 ]; %  % Input Boundary Conditions % Input DOF Constrained boundary_dof = [1 2 23 24]; % Input Constrained DOF Prescribed Values [Note Units => Length] % NOTE: For Dynamic Analysis, All Prescribed Values MUST BE Zero (0) boundary_dof_val = [0 0 0 0]; %  % Input Material Properties % Material Properties % NOTE: If the property is not being used, enter 0 (zero) % Modulus of Elasticity [Note Units => Force / Length^2] material_E = 1e6; % Poisson's Ratio [Unitless] material_v = 0.3; % Mass Density [Note Units => Mass / Length^3] material_rho = 1; % Define Material Analysis Type % Material Analysis Type Options: % Input '1' = Plane Stress Analysis (2D) % Input '2' = Plane Strain Analysis (2D) % Input '3' = ThreeDimensional Anslysis (3D) analysis_type_material = 1; 66 %  % Input Load Properties % Static Loading (If Applicable) % NOTE: If Static Analysis is not being considered, NO changes are needed if strcmp(analysis_type,'S') == 1 % Input DOF Statically Loaded load_dof = [18 20]; % Input DOF Statically Loaded Prescribed Values [Note Units => Force] load_dof_val = [500 500]; % Dynamic Loading (If Applicable) % NOTE: If Dynamic Analysis is not being considered, NO changes are needed elseif strcmp(analysis_type,'D') == 1 % Input DOF Dynamically Loaded load_dof = [44]; % Input Load TimeSeries File Name [Note Units => Force] % NOTE: First Number Designates the Load at Time = 0 seconds % Column 1 > n = DOFs % Row 1 > n = Load at Time Increments load_file = 'load.dat'; % Input Time Step of Load TimeSeries [seconds] dynamic_delta_t = 0.1; % Define NewmarkBeta Parameters [Unitless] % Average Acceleration Method => Beta = 0.25 & Gamma = 0.5 (Unconditionally Stable) % Linear Acceleration Method => Beta = 0.16667 & Gamma = 0.5 (Conditionally Stable > delta_t / Tn <= 0.551) dynamic_beta = 0.25; dynamic_gamma = 0.5; end %  % Output Options % General Plots % Plot Element Mesh (Y or N) plot_element_mesh = 'Y'; % Plot Element Mesh Options % Display Element Numbers on Element Mesh Plot (Y or N) plot_element_mesh_numbers = 'Y'; % Display Element Node Numbers on Element Mesh Plot (Y or N) plot_element_mesh_node_numbers = 'Y'; % Static Analysis (If Applicable) % NOTE: If Static Analysis is not being considered, NO changes are needed % Plot Static Displacements (Y or N) plot_static_displacement = 'Y'; % Contour Plots of Displacements (Y or N) % NOTE: Must Plot Static Displacements to Plot Contours plot_static_contour_displacement = 'Y'; % Plot Static Stresses (Y or N) plot_static_stress = 'N'; % Contour Plots of Stresses (Y or N) % NOTE: Must Plot Static Stresses to Plot Contours plot_static_contour_stress = 'Y'; 67 % Plot Static Strains (Y or N) plot_static_strain = 'N'; % Contour Plots of Strains (Y or N) % NOTE: Must Plot Static Strains to Plot Contours plot_static_contour_strain = 'Y'; % Dynamic Analysis (If Applicable) % NOTE: If Dynamic Analysis is not being considered, NO changes are needed % Movie of Dynamic Displacements (Y or N) movie_dynamic_displacement = 'Y'; % End Input File % ************************************************************************* % ************************************************************************* 68 APPENDIX B CONFIRMATION CALCULATIONS Appendix B provides confirmation calculations of the element formulation for the Q4 elements as well as a few select time intervals for the structural dynamics response of the structural system. These calculations were done by hand and then compared to the results developed by the Matlab FEA program as provided in Appendix A. 69 70 71 72 73 74 75 76 77 Comparison of Hand Calculation and Computer Results  Case 1 DOF Computer Results (in) Hand Calculation Results (in) u1 7.32E48 0 u2 2.69E34 0 u3 0.003575 0.0036 u4 0.024134615 0.0241 u5 0.00715 0.0071 u6 2.69E34 0 u7 0.007 0.007 u8 0.000184615 0.00018 u9 0.003575 0.0036 u10 0.02405 0.024 u11 0.00015 0.00015 u12 0.000184615 0.00018 Comparison of the Hand Calculation and Computer Results  Case 2 DOF Computer Results (in) Hand Calculation Results (in) u1 1.40E32 0 u2 2.69E34 0 u3 0.021 0.021 u4 0.118 0.118 u5 0.028 0.028 u6 0.376 0.376 u7 1.40E32 0 u8 2.69E34 0 u9 0.021 0.021 u10 0.118 0.118 u11 0.028 0.028 u12 0.376 0.376 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 Time Step : t = 0.1 second t = 0.2 second DOF Hand Calculation Results Computer Results Hand Calculation Results Computer Results u1 0 2.29E28 0 6.24E28 u2 0 9.15E29 0 3.37E28 u3 0.00028 0.00028 0.001 0.00141 u4 0.00037 0.000374 0.001 0.00101 u5 0.00074 0.00075 0.003 0.00342 u6 0.0059 0.0059 0.029 0.02951 u7 0 2.27E28 0 6.31E28 u8 0 1.51E29 0 2.27E28 u9 0.00021 0.00021 0.001 0.00127 u10 0.00042 0.000422 0.0009 0.00093 u11 0.00087 0.000872 0.003 0.003711 u12 0.0061 0.00611 0.0298 0.02983 Time Step: t = 0.3 second t = 0.4 second DOF Hand Calculation Results Computer Results Hand Calculation Results Computer Results u1 0 8.05E29 0 2.87E27 u2 0 4.09E28 0 4.15E29 u3 0.003 0.00365 0.006 0.00672 u4 0.011 0.01184 0.036 0.03605 u5 0.007 0.00701 0.009 0.00935 u6 0.071 0.0719 0.121 0.12166 u7 0 8.65E29 0 2.87E27 u8 0 3.02E28 0 8.43E29 u9 0.003 0.003589 0.006 0.006715 u10 0.011 0.01178 0.036 0.03603 u11 0.007 0.007254 0.009 0.009488 u12 0.072 0.07217 0.121 0.12177 95 96 97 98 Note: The displacement response shown compares the steadystate solution of the handcalculation only. The computer program results include both the steadystate response and the transient response of the system. As a result, some disagreement, particularly at the peaks of the steadystate response, is expected and confirmed in the figure. 0 1 2 3 4 5 6 7 8 9 10 2.5 2 1.5 1 0.5 0 0.5 1 1.5 2 2.5 Time [seconds] Displacement [in] Comparison of Hand Calculation(SteadyState) and Computer Results Hand Calculation Computer VITA Muhammet Ali Saglar Candidate for the Degree of Master of Science Thesis: DEVELOPMENT OF A FINITE ELEMENT ANALYSIS PROGRAM FOR STRUCTURAL DYNAMICS APPLICATIONS Major Field: Civil Engineering Biographical: Personal Data: Born in Balikesir, Turkey, August 22nd 1984, son of Remzi Saglar and Asiye Saglar. Education: Graduated high school from Ataturk Anatolian, Ankara, Turkey in June 2002; received Bachelor of Engineering degree in Civil Engineering from Istanbul Technical University, Istanbul, Turkey in June 2007. Completed the requirements for the Master of Science in Civil Engineering at Oklahoma State University, Stillwater, Oklahoma in July, 2009. Experience: Employed by Oklahoma State University as a research assistant from January 2008 to July 2009. Employed by Tekfen Corporation as a civil engineering intern from May 2006 to July 2006. Employed by Yuksel Construction Corporation as a civil engineering intern from May 2005 to July 2005. Professional Memberships: American Concrete Institute, American Institute of Steel Construction. ADVISER’S APPROVAL: Dr. Jonathan S. Goode Name: Muhammet Ali Saglar Date of Degree: July, 2009 Institution: Oklahoma State University Location: Stillwater, Oklahoma Title of Study: DEVELOPMENT OF A FINITE ELEMENT ANALYSIS PROGRAM FOR STRUCTURAL DYNAMICS APPLICATIONS Pages in Study: 107 Candidate for the Degree of Master of Science Major Field: Civil Engineering Scope and Method of Study: The finite element method is a powerful method to find solutions for engineering problems. Structural dynamics is an important concept for understanding behavior of structures under different types of dynamic loads. The combination of these two concepts into a seamless, integrated computer program to analyze structural systems considering dynamic, timedependent loads provides tremendous capabilities with respect to numerical analysis. The implementation of the theoretical development of both concepts into a fullyfunctional MATLAB computer program provided the primary objective of this study. Findings and Conclusions: The development of the MATLAB computer program to analyze a structural system using the finite element method incorporating structural dynamic response due to timevarying loads was accomplished. Several case studies were presented that illustrate the capabilities of this program with respect to determining the response of the structural system. In the current study, bilinear quadrilateral (Q4) isoparametric finite elements were used to construct the finite element mesh of the structural system. The finite element method was then integrated with a numerical timestepping method, the NewmarkBeta method, to determine the response of the structural system due to a dynamic loading. 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


