Dynamic Analysis of Multi-Degree-Of-Freedom Systems Using a Pole-Residue Method

This dissertation is concerned with the development of a robust pole-residue method for the dynamic analysis of multiple-degree-of-freedom structures. The utilization of this new pole-residue method requires additional analysis and a determination of its effectiveness in calculating a structure’s response to an external excitation. The components of the computational framework are a deterministic model of the dynamic loading and a well defined structure. The finite element method will be used to describe the system, simplify the problem, and develop the mathematical expressions to describe the system’s behavior. A commercial finite element software program was used to verify the natural frequencies of the structures and ensure that the governing equations were defined correctly. A poleresidue method was implemented and compared to the frequency domain method and a time-domain method. The system’s total response was found using the timedomain and pole-residue methods, while the frequency-domain method found the steady-state response.


Introduction
This thesis is concerned with the development of a robust pole-residue method for the dynamic analysis of multiple-degree-of-freedom structures. The strategy developed involves utilizing finite element analysis to determine the system response resulting from time-varying loads. A commercial finite element analysis software package is used to validate the obtained numerical results. [1] The study of structures is performed in the initial design phase to ensure that public safety, finished product reliability, and determination of the required construction material is well understood. In order to accomplish this feat, the dynamic loading needs to be characterized. A dynamic load is one whose magnitude, direction, or point of application varies with time. The resulting time-varying displacements and stresses constitute the dynamic response of the structure. [2] The essential characteristics of a dynamic problem differ from a static problem not only in the time-varying nature of the load, but also in the complexity of the solution for related system characteristics. A static problem has a single solution since there is no time-varying change. However, since a dynamic problem is timevarying there is a single solution per snapshot of time. [3] A static problem when subjected to a static load can be solved through using established principles of force equilibrium equations. The internal moments and shears of the deflected shape rely solely upon the static load. In a dynamic problem the resulting displacements of the system not only depend upon the load but also upon the opposing inertial forces. Thus, the resulting internal forces must be solved in equilibrium with the internal moments and shears, and not rely solely upon the dynamic loading. [4] The two basic approaches for evaluating dynamic structural response involve deterministic and nondeterministic methods. The selection of which approach to use depends on how the loading is defined. If the loading is fully known, the time-varying magnitude is understood even if the loading is highly oscillatory or irregular. The loading is then referenced as being prescribed dynamic loading.
The analysis of any system's response to prescribed dynamic loading is performed using deterministic methods. However, if the time-varying nature of the load is not known but can be described through the use of statistics, it is referenced as being random dynamic loading. The analysis of a system's response to random dynamic loading is performed using nondeterministic methods. [5] The structural response to any dynamic loading is expressed through the structure's displacement. Deterministic analysis allows for the calculation of the time-varying displacements, as well as other useful related time-varying response quantities. A nondeterministic analysis provides only statistical information about the displacements resulting from the statistically defined loading, and thus related statistical response quantities.
There are two categories of prescribed or deterministic loading, periodic and nonperiodic. Periodic loading can be described as a simple continuous width sinusoidal pulse, also known as a simple harmonic. Complex periodic loads can be described by a summation of a series of simple harmonic components through the incorporation of the fourier transform. This allows the same general procedure to be used in the analysis of any periodic loading. [6] Nonperiodic loadings are created from impulsive loadings or long-duration general loads. Impulsive responses are more commonly observed in Shock testing, typically the result of the test apparatus being hit by a large hammer or weight.
Long-duration general loads are observed during vibration testing, resulting in the shaking of the structure.
The analysis of dynamic loading encompasses three major steps: defining the analytical model, creating and applying the corresponding mathematical model, and solving for the dynamical response. The most demanding step in any structural analysis is the creation of the mathematical model used to represent the structure.
A mathematical model is developed from the idealized model of the structural system to be studied similar to the real system. This process reduces the dynamical problem to an analytical model, thus reducing the complexity of the problem and making it easier to analyze mathematically. Once the analytical model has been developed, the application of physical laws (e.g. Newton's laws, stress-strain relationships) are used to obtain the differential equations of motions representing the corresponding mathematical model. The differential equations are then solved to obtain the predicted dynamical response. [7] There are three methods for solving dynamic loading problems: lumped-mass procedure, generalized displacements, and finite element method. Lumped-mass procedure utilizes Newton's laws dealing with force and acceleration in vector quantities, to derive the equation of motion for a single particle or rigid body. A simple spring-mass model is a common example of a lumped-mass model. When using the lumped-mass model to analyze a structure, the physical components for that structure are related to displacement, velocity, and acceleration. The generalized displacements method utilizes a generalized parameter model to approximate the system's deformation. This model is used for a continuous system, whose deformation is described through one or more functions of space and time. In order to be solvable, the model requires that a geometric boundary condition be specified, constraining the kinematic displacement on a portion of the structure. The finite element method utilizes key components of both lumped-mass and generalized dis-placements methods. This method discretizes the structure in multiple segments or elements and solves for multiple highly coupled interpolation functions describing the structure's geometry. [8] This thesis will implement a pole-residue method to determine the system response. In order to understand the effectiveness of any new method, the results must be compared against a well known and understood method. This comparison will be performed using the assumed-modes method which combines the lumpedmass procedure with Newton's laws for the analysis of a six-degree-of-freedom mass-spring-dashpot system, and the finite element method to investigate a simple five story frame structure. In order to depict conceptual understanding of the topic at hand, this method will first be implemented in Matlab and later verified in a commercially available finite element analysis software application, such as Abaqus.

CHAPTER 2 Structural Dynamics
In order to perform dynamic analysis, several concepts need to be understood.
There are a number of factors that influence which model is implemented. The complexity of the system and the number of influencing parameters need to be determined to appropriately decide on which model to use. This section will cover the assumed-modes method, the finite element method, complex analysis using the Fourier transform, as well as the commercial finite element analysis software Abaqus.

Single Degree of Freedom Systems
The fundamental equation in structural dynamics and linear vibration theory is expressed as second-order differential equation which relates force to displacement, velocity, and acceleration This equation of motion represents a single lumped-parameter model and is derived using Newton's laws, (Figure 1).

Figure 1. Single-Degree-of-Freedom Lumped Parameter Model
The elements which comprise this model include the mass, the spring, and the damping component. According to the model, an object is placed into motion through the displacement of its mass at a specified distance from a reference. The spring element then goes into tension due to being stretched, and it will have the desire to compress and restore the system to equilibrium. The damping component, which occurs when impacted by the spring's restoring motion, will dissipate system energy decreasing the time required to reach equilibrium.

Newton's Laws to Lumped-Parameter Models
By connecting several lumped-parameter models together, a simple singledegree-of-freedom model can be constructed into a more complex multiple-degreeof-freedom model. Depending on the complexity of this model, the correct mathematical models can then be developed through the application of Newton's laws and/or the assumed-modes method. Newton's laws are best used on simpler models comprised of a few lumped parameter models described by a mass, spring, and damping values. For more complex systems, the assumed-modes method is preferred due to the influence of additional parameters on the system.

Newton's Laws
In this section, the mathematical model for a multiple-degree-of-freedom model will be developed using Newton's laws, (Figure 2). The notation m i represents mass, c i represents a dashpot, and k i represents a spring. The force related to the dashpots are represented by F ci and the force related to the springs are represented by F si .  After each free-body diagram has been drawn and annotated, the equation of motion for each particle must be defined per Newton's second law The mass element can easily be plugged into the equation of motion. However the linear elastic spring and viscous damping forces still need to be related to the displacements u i for each mass element. This relationship is illustrated for the three springs by This relationship is described for the three dashpots through Incorporating the defined spring and dashpot equations of motions with respect to the displacement u i with the previously defined governing equations of motions, the mathematical model for the multiple-degree-of-freedom system is defined through the use of matrix notation. This will allow for the evaluation and determination of the exact solution. These matrices define the mass M, stiffness K, damping C. The corresponding excitation vector is defined as p.
The mass matrix M is The damping matrix C is The stiffness matrix K is The excitation vector p is

Assumed-Modes Method
The assumed-modes method incorporates an extension of the virtual displacement method to produce a generalized-parameter model of a continuous system in order to approximate the deformation of that system. In order to complete this transformation, several definitions are required to define the generalized-parameter single-degree-of-freedom model. These definitions are: 1. A continuous system is a system whose deformation is described by multiple spatial variables and time To create the generalized-parameter single-degree-of-freedom model of a continuous system, a single assumed mode is used The shape function should represent a single deformation experienced by the structure. The generalized displacement coordinate is defined as q v (t) with the subscript v denoting the generalized coordinate in relation to the physical displacement v(x, t). The shape function ψ(x) can be represented by any function, however a shape representing the deforming structure should be chosen. To generate a multiple-degree-of-freedom model of a continuous system, Eq 15 is expanded to include N shape functions, thus allowing the continuous displacement to be ap- The assumed-modes method consists of substituting this equation into the expressions for kinetic energy, T , and strain energy, V, through representation of the material by its density ρ, area A, length L, and modulus of elasticity E, as shown below The analyst defines the N -degree-of-freedom assumed-modes through the selection of the shape function. The shape functions must form a linear independent set and must possess derivatives up to the order appearing in the strain energy V.
The shape function must also satisfy all the prescribed boundary conditions.

Finite Element Modeling of Structures
The finite element method is a far more powerful approximation method than the assumed-modes method. Through the approximation of the displacement function by using a continuous function, the deflected shape of the entire structure is then described by each shape function ψ The N -degree-of-freedom mathematical model of the structure is obtained, using the generalized stiffness, mass coefficients, and generalized forces obtained from integrals involving the ψ i 's and their derivatives. While using the computer to solve difficult engineering problems has facilitated the evaluation of these integrals, the draw backs in the assumed modes method: 1. It is difficult for an engineer to choose a set of ψ i 's for a structure that has complex geometry 2. The equations that result from the assumed mode method are coupled 3. There is little reuse from one problem to the next The finite element method overcomes these difficulties and has become the main computational tool for structural dynamics. The finite element method utilizes a set of ψ i 's to represent a portion of the structure through the assembling of elements to form the structure. The elements are linked together at nodes and the displacement compatibility is enforced at these linkages.

Axial Motion
Axial motion is defined through the study of a uniform bar element of length L, mass density ρ, modulus of elasticity E, and a cross-sectional area A. The approximation of the axial displacement u(x, t) within the element employs linear interpolation between the displacements u 1 (t) and u 2 (t) at the two ends and shape functions ψ 1 (x) and ψ 2 (x) The shape functions ψ 1 (x) and ψ 2 (x) must satisfy the boundary conditions, which are derived by considering axial deformation under static loads and are used to approximate axial displacement along an element of length L.
When the axial motion of the plane cross section at x is represented by u(x, t) where the generalized coordinates are labeled u i (t), Eq. 16 can be written as When Eq. 21 is subsituted in Eq. 18 the resulting expression for strain energy Where the stiffness coefficients for axial deformation is given by In a similiar fashion, the combination of Eq. 21 with Eq. 17 produces the following quadratic expression for the kinetic energy of a bar Where the mass coefficient is When the mass matrix and stiffness matrix are formed through an equation like Eqs. 23 and 25, they are defined as being consistent, as they were formed through a similiar manner using the same shape functions.
If the bar experiences an external excitation force, then the generalized forces are determined by employing virtual work δW Thus, when the virtual work equation, Eq. 26, is integrated with the displacement The generalized force for each shape function is obtained from Expressions for the stiffness coefficients k ij , the mass coefficients m ij , and the generalized forces p i (t) for axial motion were defined earlier in Eqs. 23, 25, and 28. Through substituting ψ 1 and ψ 2 into those equations, the stiffness and mass matrices for a uniform element undergoing axial deformation is defined by

Transverse Motion
Transverse motion is defined through the study of a uniform beam element of length L, mass density ρ, elastic modulus E, cross sectional-area A, and moment of inertia I. The approximation of the transverse displacement v(x, t) with the element utilizes linear interpolation between the displacements v 1 (t) and v 2 (t) at the two ends and shape functions ψ 1 (x) and The shape functions ψ 1 (x) and ψ 2 (x) must satisfy the boundary conditions.
The shape functions are derived by considering the beam element to be loaded statically by end shears and bending moments to produce the various static deflection shapes. Expressions for the stiffness coefficients k ij , the mass coefficients m ij , and the generalized forces p i (t) for axial motion was defined earlier in Eqs. 23, 25, and 28. Through the substitution the boundary conditions for ψ 1 through ψ 4 into those equations, the stiffness and mass matrices for a uniform element undergoing transverse deformation is

Torsion
Torsional deformation is defined through the study of the rotation along the element of a uniform rod of length L, mass density ρ, polar moment of inertia I p , and torsional stiffness GJ. The strain V and kinetic energies T for pure torsion The rotation θ(x, t) along the element is given by the assumed-modes form where the torsion members are represented by θ 1 (t) and θ 2 (t) and shape functions ψ 1 (x) and The shape functions ψ 1 and ψ 2 must satisfy the boundary conditions. The shape functions are derived by considering a torsion member in equilibrium but under static end torques. Combining Eqs. 34, 35, and 36 leads to the development of the mass, stiffness, and force expressions shown below Through the incorporation of shape functions, the stiffness matrix and mass matrix for uniform torsion can be determined for an element by

Three Dimensional Frame Element
The stiffness and mass matrices for a three-dimensional frame element can be obtained from the axial, bending, and torsion elements. Bernoulli-Euler beam theory was used in deriving the bending elements. A three dimensional frame element, (Figure 4), shows the local references and the displacement coordinates u i . The x axis follows along the centroids of the beam, while the y and z axis are principal to the cross section. Displacements u 2 , u 6 , u 8 , and u 12 are based on bending in the x-y plane, Eqs 32 and 33. While displacements u 9 , u 11 , u 3 , and u 5 are for bending in the x-z plane.
Finally, the coefficients associated with torsional degrees of freedom, u 4 and u 10 are obtained from Eqs. 40 and 41. Combining axial, transverse, and torsion, a three dimensional stiffness matrix and mass matrix for a uniform three-dimensional element can be determined

Vibration Properties of MDOF Systems
Mass matrix M and Stiffness matrix K are related to strain energy and kinetic energy, respectively. M and K are positive definite matrices for most structures.
This implies that any arbitrary displacement of a system with a positive definite K strain energy will be positive and that any arbitrary velocity displacement of a system with positive definite M will result in a positive kinetic energy.

Eigensolution of Undamped Free-Vibration
The free vibration equation of motion for a multiple-degree-of-freedom system is written in symbolic matrix form M and K are N × N matrices and u (t) is the corresponding N × 1 vector of physical or generalized displacement coordinates. When harmonic motion defined by displacement U and natural frequency ω When harmonic motion is then incorporated into Eq. 44 the N th order eigen value problem is formed For the solution to be nontrivial, the determinant must equal zero This is known as a characteristic equation. When the determinant is expanded, the corresponding values for the squared natural frequencies, ω 2 n , are determined from the roots of the polynomial. Corresponding to each eigenvalue, ω 2 n , there will be an eigenvector, or natural mode U n .

Eigensolution of Damped MDOF Systems
When a damping component is incorporated into the undamped N -degree-offreedom matrix equation of motion, Eq. 44, through the inclusion of the viscous dampening component C then the equation of motion becomes A procedure for determining the natural modes, ω 2 n , requires the transformation of the N second-order equations of motion into a 2N first-order equations.
The generalized state-space form is developed when the displacement vector is transformed into the state vector Which is comprised of the displacement vector and an effective velocity vector The effective velocity vector v could also be defined asu. Where A and B are the two 2N × 2N constant coefficient matrices and F is the 2N state forcing vector shown below Since this is a homogeneous set of ordinary differential equations with constant coefficients, the solution is comprised of the eigen values, λ, and the eigen vectors, Through the division of Eq. 55 by e λt the generalized algebraic eigenvalue equation is obtained by The solution for both the eigenvectors and eigenvalues consist of i = 1, 2, ..., 2N and must satisfy the characteristic equation

Response of Systems to Periodic Excitation 2.5.1 Complex Fourier Series
A common and reliable method used in structural analysis to separate the harmonic components found in these excitations is the Fourier series The fundamental frequency Ω 1 (rad/s) is related to the period of the function by Ω 1 T 1 = 2π. The bar over a parameter symbolizes that the coefficients may be complex. When the excitation is not periodic, it can be represented by a Fourier integral. In developing expressions for the Fourier integral transform pair, it is convenient to use complex notation by allowing the period, T 1 , approach infinity.
Determining the structural response requires the multiplication of the excitation and the impulse response in the frequency domain and then performing the inverse Fourier transform to return to the time domain ( Figure 5).

Figure 5. Fourier Transform Process
The complex Fourier expansion of p(t) is achieved by multiplying Eq. 58 by e i(mΩ 1 t) and integrating over one period The complex conjugate pair ofP n is represented byP −n =P * n and that the average value of p(t) can be found through

Complex Frequency Response Function
The steady state response,ū(t) shown below is in complex form The frequency responseH is defined by Ω the forcing frequency and ω n the undamped natural frequency of the system is When the excitation is periodic as shown in Eq. 58, the steady state response can be written in complex notation as The corresponding expression for periodic response becomes The process for determining the response through the use of the Fourier series is shown in Figure 5. The method utilizes Eq. 59 to convert the forcing excitation such that it relates the magnitudes to their respective frequencies. Eq. 64 is used to couple the harmonic excitationP n and the frequency responseH n .

Fourier Integral Transform
The Fourier coefficientP n in relation to p(t) was shown in Eq. 59, providing that the integral exists. If you allow T 1 → ∞ then it is convenient to allow the following notation: Ω 1 = ∆Ω, Ω n = nΩ 1 ,P (Ω n ) = T 1Pn = 2π ∆ΩP n then Eq. 58 becomes Through the same notation the inverse is The limits of integration will include the entire time history of p(t) when T 1 → ∞. As T 1 → ∞, Ω n can be represented as a continuous variable Ω and ∆Ω becomes the differential notation dΩ. Thus allowing Eqs. 65 and 66 to be represented respectively by These two equations are the Fourier transform pairs.P (Ω) is known as the Fourier transform of p(t) and p(t) is called the inverse Fourier transform ofP (Ω).
Finally, the Fourier Transform pair can be written in terms of frequency, when the relationship f = Ω 2π . Then Eqs. 67 and 68 can be written to illustrate the transform to and from the frequency domain

Frequency Response Function
The periodic excitation can be expressed as Eq. 63 and the relationship between the excitation, the periodic response, and the system response can be expressed as Eq. 64. However, if the relationships used to convert the Fourier transform for a periodic excitation function are incorporated into the response, then we obtain the Fourier transform pair for u(t) andŪ (f ) The relationship Eq. 64 thus becomes Which is the product of the system frequency response functionH(f ) and the Fourier transform of the excitationP (f ).
Therefore, the response can be expressed as the following inverse Fourier transform

System Response Function
The system response function,H(f ) can be obtained from the Fourier transforms of the measured time histories p(t) and response u(t). The system parameters for natural frequency ω n and the damping factor ζ can be extracted from the system frequency response function. However, other analytical methods are available for estimating the system response if either p(t) or u(t) is unknown, such as the mode-superposition method.

Mode-Superposition for Undamped Systems
The mode-superposition method is used to obtain the response u(t) through the use of the natural frequencies and natural modes of the system. The natural frequencies and modes must satisfy the algebraic eigen problem, shown earlier in Eq. 57. Given the natural frequency, ω 2 r , and the modes, φ r , for r = 1, 2, ..., N , the modal masses, M r , and modal stiffness, K r are calculated through Any natural frequencies that are redundant or repeated have been orthogonalized by This is satisfied when r = s. The modes are then gathered to create the modal Using the modal matrix we can rewrite Eqs. 75 and 76 for the modal mass and modal stiffness to be The key step of the mode-superposition method is the coordinate transform shown by The coordinates η r (t) is referred to as a principal coordinates or modal coordinates and can be defined by where: The generalized modal mass matrix is The generalized modal damping matrix is The generalized modal stiffness matrix is The generalized modal force vector is

Mode-Superposition for Damped MDOF Systems
When a viscous damped system satisfies the orthoganality requirements for Eq. 81, then the uncoupled modal equation of motion can be written as If the system is under harmonic excitation, p(t) = P cos(Ωt), then the modesuperposition solution for the steady-state response becomes Where F r = φ T r P. Using the complex-frequency response techniques shown earlier, Eq. 87 can be defined asη The steady-state solution forη r becomes Where r r is the modal frequency ratio, r r = Ω ωr . The magnitude and phase of Eq. 89 can be found through The complex frequency response in generalized coordinates can be obtained by writingū(t) in the complex form Combining the above equations a solution forū(t) is shown in to bē The complex frequency response function in physical coordinates,H ij (Ω) gives the response at coordinate u i due to harmonic excitation at p j Where φ ir is the element in row i of the r-th mode φ r , that is the element in row i and column r of the modal matrix Φ. The steady-state response u(t) can then be obtained from

Frequency Response of MDOF Systems based on Complex Modes
The mode-superposition method also applies when a state space form which is used to transform N second order differential equations to 2N first order differential equations. The complex modes from the state-space problem have the following λ r is the complex scalar eigenvalue and θ r is the corresponding 2N state eigenvector. The subscripts u and v reference the displacement and velocity, respectively. The eigenvalue, λ r can be written in the following forms The natural frequency, ω r , and the dampening factor, ζ r are given by The modal matrix containing the state eigenvectors is The orthogonality equations can be combined with definitions of diagonal terms modal a r and modal b r Such that λ r = −br ar . The harmonic excitation at the forcing frequency, ω, will have the steady state response z = Ze iωt . Where Z can be expressed in the modesuperposition form The displacement partition can be extracted by The complex receptance frequency response function is (θ ui ) r is the i-th element in the displacement partition of the r-th eigenvector. Since the eigenvectors and eigenvalues occur in complex conjugate pairs, Eq. 105 can be written as

Abaqus
Abaqus is a commercial finite element analysis software application that can be purchased from Dassault Systemes. Abaqus provides product simulation to engineers in order to vary and simulate various design attributes. The Abaqus Unified FEA product suite offers powerful and complete solutions for both routine and sophisticated engineering problems covering a vast spectrum of industrial applications. Abaqus allows engineering working groups to consider full loads, dynamic vibration, multibody stems, impact/crash, nonlinear static, thermal coupling, and acoustic-structural coupling using a common model data structure and integrated solver technology. Abaqus provides four analaysis packages to the user. Abaqus utilizes three steps to develop, perform, and present the results for each finite element model. The first step is the pre-processing or modeling stage.
This involves the creation of an input file which contains the design of the system, the system parameters, and the analysis to be performed. The second step is the processing of the model which is performed behind the scenes. The final step is post-processing where a report, image, animation, data table, etc is provided from an output file. Abaqus CAE provides the capability for pre-processing, postprocessing, and monitoring of the processing stage.

Abaqus Spring-Dashpot Commands
Abaqus provides a graphical user interface for the design and evaluation of mechanical systems. These systems are typically defined by the component(s) being drawn or integrated and their specific material properties. However, Abaqus does provide the capability to evaluate mass-spring-dashpot systems through the use of specific definitions and commands. In order to incorporate spring and dashpots to a system these specific functions must be called. In the example provided in this paper, these commands are implemented behind the scenes in the Abaqus input file, (Figure 6).  The specific commands used to generate this model are: 1. * M ass allows for the introduction of a concentrated mass at a point and the ability to associate with the three translational degrees of freedom at a node.
2. * Spring allows for a spring element to be defined. This spring can couple a force with a relative displacement.
3. * Dashpot allows for a dashpot element to be defined. This dashpot can couple a force with a relative velocity.
4. * RotaryInertia allows for rotary inertia to be included at a node. This inertia can be associated with the three rotational degrees of freedom at a node and can be paired with a mass element.

Three Dimensional Frame Element
Abaqus provides a complex simulation environment for the design of structures. Through varying specific attributes of the structure an engineer can determine the most feasible option for their client through material selection and the method used for joining linkages. In the example provided in this paper a 3-Dimensional Framed Structure is evaluated, (Figure 8). The specific commands used to generate this model are: 1. * N ode allows for the user to define the nodes by directly specifying its coordinates. Nodal coordinates are given in a local system.
2. * Element allows for the definition of an element by specifying its nodes.
3. * F rameSection allows for the definition of the cross-section for frame elements.

Analysis
Abaqus provides an extensive suite of analysis tools to the user to evaluate the system being considered.
For the examples utilized in this paper the analysis focused on finding the natural frequencies. The command used to perform this task is the * f requency. This command tasks Abaqus to perform eigenvalue extraction to calculate the natural frequencies and the corresponding mode shapes of a system.
The analysis will include initial stress and load stiffness effects due to preloads and initial conditions if geometric nonlinearity is accounted for in the base state.

Output
Abaqus has several reporting features which can be used to extract solutions for the model, system, or at the declared nodes. The * Output command is used to write contact, element, energy, nodal, or diagnostic output. There are various commands that can be used in parallel to efficiently extract the desired parameters.
The specific output commands that are available are: 1. * Contact Output allows for the extraction of variables associated with surfaces in contact, coupled temperature-displacement, coupled thermalelectrical, and crack propagation problems.
2. * Element Output allows the request of element variables such as stresses, strains, section forces, element energies, etc.
3. * Energy Output allows for the extraction of the total energy of the model or of a specific element.

* Modal
Output allows for the output in generalized coordinates the modal amplitude and phase values.
5. * Node Output allows for the displacements, reaction forces, etc. that describe the nodal variables.
6. * Radiation Output allows for the extraction of cavity radiation variables.

Matlab State-Space Methods
MatLab provides several built in functions which allow for the modeling of complex systems through the use of the appropriate mathematical models derived from physical laws. Dynamic systems change in time per the set mathematical model. The functions to be reviewed in this section are the ss and lsim functions.
When utilizing these MatLab functions for a mass-spring-dashpot system it is best to construct a free body diagram and follow the procedures shown earlier.

MatLab Function ss
The MatLab function ss is used to create state-space models with real or complex matrices. This function takes the following inputs: 1. a: a state matrix that is square real or complex with as many rows as states such that it is X by X 2. b: an input state matrix that is real or complex with as many rows as states and as many columns as inputs and its size is X by U 3. c: a real or complex state to output matrix that has as many rows as outputs and as many columns as states and its size is Y by X 4. d: a real or complex feedthrough matrix that has as many rows as outputs and as many columns as inputs and its size is Y by U Where X is equal to the number of states, U is equal to the number of inputs, and Y is equal to the number of outputs. For dynamic systems, governing equations are determined through the summation of forces following Newtons second law.
A state-space representation of the system is then developed in order to reduce the N -second order differential equations to 2N -first order differential equations.
Once the state space representation is gained, the information for a, b, and c is available as an input. Since feedback is not being studied, d is set to zero.

MatLab Function lsim
The MatLab function lsim is used to simulate the time response of continuous or discrete linear systems to arbitrary inputs. This function returns the system response, sampled at the same rate as the time history. The output will have the same number of elements as the time history and the same number of columns as the system outputs. Lsim uses a linear interpolation (first-order hold) method based on the smoothness of the signal u. This function takes the following inputs: The numerical methods used to solve the system response of multiple-degree of freedom systems primarily operate in the time or frequency domain. While more complex problems used the Fourier transform or numerical integration procedures in the time domain, the Laplace transform has been used for some simple problems which require the use of look-up tables. This section will detail how to implement a semi-numerical pole residue method to solve more complex problems.

System Response by the Laplace Transform
Transforms are used to simplify the solution of certain dynamics problems.
Their purpose is to shift a complex problem from one domain to another, allowing a solution to be found more readily. Once the solution is obtained, these results can be transferred back to the original domain. The Fourier transform allows the problem to be shifted from the time domain to the frequency domain and back.
The Laplace transform shifts the problem from the time-domain to the s-domain.
The Laplace transform is particularly useful if the forward and inverse transform can be identified in a lookup table. However, this is typically not the case with the vast majority of engineering problems. In this section a new pole and residue method is to be discussed that falls within the category of Laplace domain method.
Thus far, this method has proven to be applicable to arbitrary input functions, but it has yet to be utilized in the determination of the system response of a structure.
This method focuses on obtaining the poles and residues of the response function.
Knowledge of these parameters will allow one to express the response function in both the time and s domain. The following steps are required for this method: 1. Preparing the poles and residues of input and system transfer functions 2. Conducting algebraic computation to obtain the poles and residues of the response function based on those of the input and system transfer functions 3. Providing a solution in the time domain from the poles and residues of the response

Utilizing a Prony's Method for Developing the Poles and Residues
Structural excitations can be represented as an arbritrary signal y(t). This signal can be decomposed into a finite number of exponential components In this equation, N is the number of terms, α n and λ n are constants. When excitation y(t) is a real-valued signal, λ n must be real or form complex conjugate pair. When λ n is equivalent to −δ n + iΩ n , then δ n is the damping factor in seconds −1 and Ω n represents damped frequency in radians per seconds. The coefficients α n corresponds to the complex exponents, thus λ n must appear in complex conjugate pairs. If α n is equivalent to A n e iθn , then A n represents the amplitude, and θ n is the sinusoidal initial phase in radians in association with e λnt .
When digital signal analysis is performed the continuous y(t) its components are not known and the captured signal is usually represented by equally spaced samples. The Prony series corresponding to this discrete signal is A Prony series has its components being complex exponentials, which could be purely harmonic, damped harmonic, or purely damped exponential components.
A Prony series is comprised of a smaller subset of terms, p, in comparison to the sampled signal length, N . Most often p is much smaller than N . Finally, z n must be estimated first from the discrete signal y k before computing the complex coefficients γ n .

Differential/Difference Equation for Prony Series
The arbitrary signal y(t) contains the exponenetial functions e λnt , for n = 1, 2, ...N which forms a basis on the open internval 0 ≤ t < T . Thus allowing the Prony series, Eq. 107, to be viewed as the general p-th order homogenous linear ordinary differential equation with constant real valued coefficents a n p n=0 a n y n (t) = 0, for 0 ≤ t < T In Eq. 109, y n (t) = d n y(t)/dt, the exponents λ n for n = 1, ..., N are the roots of the characteristic equation

Implemented Prony's Method
A realization of A is possible through the substitution of the discrete signal y k into the Hankel matrix Upon the substitution of y k into H(k), it can be shown that H(k) is the product of three matrices Where and Applying singular value decomposition of H(0) yields The model order of the dynamic system, p, is equal to the number of non-zero singular values. The singular values should go to zero when the rank of the matrix is exceeded. However, when performed on measured data, the singular values will not become zero due to random errors and inconsistencies and will instead become very small. When choosing the size of H(0) both ξ and η must be greater than p.
A direct comparison between Eq. 119 and Eq. 121 suggests that P ξ is related to U 1 and Q η to V T 1 . Substituting these relationships into Eq. 120 leads to Rearranging Eq. 122 and premultiplying by S 1 2 1 U T 1 and post multiplying by The computed eigenvalues of A are z n , for n = 1, ..., p corresponding to λ n .
Finally the complex coefficient γ n , for n = 1, ..., p can be computed in the same way as that used in the original third step of Prony's method.

Prony's Method to Laplace
It is evident that the exponential function αe λt and the pole-residue form α s−λ are a Laplace transform pair. Thus, after signal decomposition has been performed using the modified Prony's method, the Laplace transform yields the following equation in the s-domainỹ

Integration of Pole and Residues with Structural Systems
In the previous sections, background information on the laplace transform and the pole-residue method was presented. Through the coupling of this information with material already presented in Chapter 2, another method capable of solving complex structural problems will be presented.

Equation of Motion
The governing equations for a N -degree-of-freedom system is This second order differential equation is comprised of matrices for M, C, and K representing the mass, damping, and stiffness. These matrices are N × N .
Through the annotation of T(s) with (j, k), one can determine the appropriate transfer function for a response at coordinate, j, caused by the impact of an input force at coordinate, k. There are two scenarios to be considered when deriving the pole-residue form of T jk (s). The first is when the problem is composed of a symmetric M, C, and K matrices and the second is the asymmetric case.

Symmetric System Response
The majority of structural dynamic problems are composed of symmetric M, C, and K. This is especially true when the matrices are developed through a finite element procedure. In order to determine the system response, the N -second order differential equations will need to be simplified into 2N -first order equations. This can be performed through the procedure mentioned earlier in Section 2.4. Since M, C, and K are symmetric, the state space model maintains this property through coefficient matrices A and B. From Eq. 49, the generalized eigenvalue equation is The solution consists of 2N eigenvalues φ n , which are real or occurs in complex conjugate pairs The first N modes can be extracted from the first N elements of the eigenvector θ n . Since A and B are symmetric the orthogonality property, θ T r Aθ s = 0 and θ T r Bθ s = 0 must hold true for µ r = µ s . Through the utilization of the modal The physical-modal coordinates transformation relationship is One can then convert, Eq. 49 into 2N uncoupled modal equations of motion where a n = θ T n Aθ n , b n = θ T n Bθ n , and g n = θ T n g a nẏn + b n y n = g n (134) From Eq. 134 the transfer function of the n-th modal equation is found to be T n (s) = 1 a n s + b n = 1 a n (s − µ n ) The pole, µ n is defined as being equal to −bn an . One can express the transfer function T n (s) as T jk (s) through the use of mode super-position as Where the residue β jk,n is β jk,n = θ k,n θ j,n a n (137) While poles are global parameters and are independent of the input and output coordinates, j and k, the residues, β jk,n are dependent on j and k.

Asymmetric System Response
When M, C, and K, are not symmetric, another approach is required to determine the Transfer Function, T jk without the use of coordinate transformation.
An intermediate step is required to gain the zero-pole form, prior to obtaining the pole-residue form T jk (s). Since the Transform, T jk (s), contains the inverse of Z, the explicit formula for the inverse of a nonsingular square matrix, Z The cofactor, C jk is equal to (−1) j+k det(Z jk (s)), where Z jk denotes the submatrix obtained from Z by deleting the j-th row and k-th column from Z. From Eq. 138 the (j, k) elements of T can be written as Note the subscript for T and Z are reversed. Since Eq. 128, is an N x N matrix, for a non-singular matrix Z (s) the determinate must be a 2N -order polynomial in the s-domain. The determinant of the non-singular matrix for Z jk (s) must be a 2(N − 1) order polynmial in the s-domain. The zero-pole form for T jk is The remaining task is to compute poles, µ n , the jk-dependent zeroes τ jk,n and leading coefficient κ jk from M, C, and K. The poles, µ n are computed by carrying out eigen analysis of a state space model. Since M, C, and K are real-valued matrices, the solution to the generalized eigenvalue problem yields 2N eigenvalues.
Computing the zeros, τ jk,n follows the same procedure as computing the poles by replacing M, C, and K by M jk , C jk , and K jk , respectively. The coefficient for κ jk can be computed using Eq. 141, through the multiplicative rule of the determinants Eqs. 142 and 143.
Once the zero-pole form has been obtained through Eq. 140, the corresponding pole-residue form can be found to be The jk dependent residues of β jk,n is

Displacement Response
The displacement response at coordinate j to a force at k is given by the convolution integral of the system and force The total response at coordinate j produced by a force involving all components of the load vector is obtained through the summation all load components and thus capturing their individual contributions The pole-residue analysis is similar to the frequency and analytical procedures in that superposition is performed. The resulting s-domainx jk (s) response to the displacement at the j-th coordinate due to a force at the k-th coordinate is expressed asx The total response at coordinate j produced by a general loading at all k can be obtained to be Eq. 149 forms the general solution to the coupled equation of motion, assuming zero initial conditions. In order to successfully perform Laplace analysis, one must be able to efficiently compute the pole-residue form of the response,x jk (s).
In order to implement Eq. 148 and compute the pole-residue of the response,f k (s) and T jk (s) can be expressed asf When Eqs. 150 and 151 is substituted into Eq. 148 the following form for When Eq. 152 is transferred into the pole and residue form, it is evident that there are a total of 2N + L poles resulting from the L excitation poles and 2N system polesx The first L poles correspond to the excitation poles such that ν m = λ m , where m = 1, 2, ..., L. The remaining 2N poles corresponding to the system poles such that, ν m+L = µ m , where m = 1, 2, ..., 2N . The corresponding residue can be found for each response pole, ν m by Through Eq. 154, the residues corresponding to the first L response poles can be found using The residues corresponding to the last 2N response poles can be found through The paper analyzed a real world six degree of freedom problem that is comprised of two closely spaced modes, two weak modes, and a wide range of modal damping.
The authors determined that the sensitivity of implemented methods was a direct result of the model order, design parameters, and user expertise.

Problem Overview
The structure presented is a simplified model of one-half of a railway vehicle, ( Figure 9). The railway car is represented as a rigid body with two degrees of freedom, vertical displacement u 3 and pitch angle u 4 . The railway car is connected to two nodes located at the front and rear end of the car representing the secondary suspension. These nodes are modeled as a rigid body with a single degree of freedom, comprised of a vertical displacement u 2 or u 4 . The secondary suspension nodes are connected through the primary suspension to a wheelset modelled as an unsprung mass with a single degree of freedom, comprised of a vertical displacement u 1 or u 6 . The interaction between the wheel set and the track is modelled by linearized springs. The track is assumed to be fixed and rigid. As one-half of a railway vehicle the values for the system parameters are nearly symmetrical with respect to the vertical axis through the car's center of mass.   Through the combination of the related forces to their respective displacement and velocities, a mathematical model can be written in matrix notation. The mode shapes, modal frequencies, and modal dampings, can be determined after transitioning the problem to a state space model.

Free Body Diagram
A free body diagram is developed utilizing the process depicted in Section 2.2.
Through the completion of this activity the equation of motion for each particle can then defined, thus allowing for the creation of the mathematical model representing this problem, (Figure 10).

Equation of Motions
Now that the free body diagram for each element has been drawn and annotated, we can now define the equation of motions as shown below.

Relation to Displacements and Velocities
As mentioned earlier in Section 2.2 the mass element can be easily plugged into the equations above. However, the linear elastic spring and viscous damping forces still need to be related to the displacements u i for each mass element. This relationship is illustrated in the equations below for the six elastic spring forces.
This relationship is illustrated in the equations below for the six viscous damping forces.
Now that we have defined the relationship between the linear elastic spring and viscous damping forces to the displacements, we can complete the equations of motion and complete the mathematical model.

State Space Model
Through the use of the eigensolution of damped multiple-degree-of-freedom systems, the modal frequencies and modal dampings were found through the use of a generalized state space model.

Abaqus
In order to ensure that the modal frequencies that were calculated in Matlab are correct, this six degree-of-freedom model was implemented in Abaqus and analyzed through the command prompt. An input file representing the model was developed and analysis was performed through the command prompt was launched.

Impact and Measurement Locations
In order to analyze a multiple degree of freedom problem and determine the response, the two key parameters are the input, p, and the output, q. In this example the input location was chosen to be node 2 and the output location was node 4. The displacement response is to be represented as x 42 .

Excitation
An excitation, p(t) was developed utilizing the parameters for the frequencies components, amplitudes, phase angle, and decay factor, (Table 6). The excitation was generated using A n e δnt cos(2πf n + θ n ) The three component signals are shown in red and the excitation is shown in black, ( Figure 11).  The poles and residues components of the excitation were found using the modified Prony's method , (Table 3). The poles and residues for the system transfer function were found following the symmetric process. The poles and residues were found by constructing the state space model from M, C, and K and then solving the generalized eigenvalue problem, (Table 4).

Steady-State Response through the Poles and Residues
Once the poles and residues of the excitation and the transfer function are known, then determining the poles and residues of the time history forũ 4 can be accomplished, (Table 5).
From these corresponding poles and residues the time history u 42 (t) can be obtained, (Figure 14).

Comparison of Results
The displacement response x 4 due to the excitation at node 2 was found through the use of a frequency-domain, time-domain, and pole-residues methods.
The displacement response from these three methods do agree, (Figure 15). The time-domain and pole-residue methods agree strongly since both method depict the steady-state portion of the response, but also the transient component.
The frequency-domain result is composed of only the steady-state solution.

Function Completion Times
The time to complete the execution of the above methods for the one-half of a railway cart was determined through the use of the profile function in MatLab, ( Figure 16). The problem being examined in this chapter is a three dimensional frame structure. This problem represents a real world multiple-degree-of-freedom system depicting a five story building with four fixed points representing the earthstructure interface. There are four fixed nodes, twenty nodes that are not fixed, and forty linkages. This structure will be analyzed through the use of Bernoulli-Euler beam theory.

Problem Overview
The structure is a simplified five story building, (Figure 17). The five story building is represented as a three dimensional frame element with six degrees of freedom. The building will suffer from axial motion, transverse motion, and torsion. The axial motion will be defined by the structures cross sectional area A, Elastic Modulus E, and mass density ρ. The transverse motion is governed by the structure's mass density ρ, elastic modulus E, cross sectional area A, and moment of inertia I. This structure will undergo torsion which is governed by the shear modulus of elasticity G, elasticity J, and polar moment of inertia I p . The linkages are modelled as a three dimensional frame element with 12 degree of freedoms, six displacements, and six moments.

Matlab Finite Element Analysis Solution
One solution to this problem is to implement the finite element method.
Through the representation of the structure as an assembly of elements the displacement can then be found. In order to perform this task, the structure needs to be defined in such a way that the following tasks could be accomplished: 1. The structure can be modeled 2. The parameters for each beam and the corresponding beam elements can be defined 3. This finite element problem was implemented in Matlab utilizing the three dimensional 12x12 matrices for mass and stiffness 4. The eigenvalues were found using the undamped free-vibration eigensolution

State Space Model
The M and K matrices were generated utilizing the finite element method.
The undamped eigensolution was then used to determine the natural frequencies

Impact and Measurement Locations
In order to analyze a multiple-degree-of-freedom determine the response, the two key parameters are the input, p, and the output, q. In this example the input was chosen to be node 20 and the output location was node 14. These two nodes were selected randomly. The displacement response is to be represeneted as x 14 20 .

Excitation Signal
An excitation, p(t) was developed utilizing the parameters for the frequencies components, amplitudes, phase angle, and decay factor, (Table 6). The excitation was generated by A n e δnt cos(2πf n + θ n ) The three component signals are shown in red and the encompassing excitation signal is shown in black, (Figure 18).  The poles and residues components of the excitation were found using the modified prony's method , (Table 7). The poles and residues for the system transfer function were found following the symmetric process. The poles and residues were found by constructing the state space model from M, C, and K and then solving the generalized eigenvalue problem, (Table 8). Once the poles and residues of the excitation and the transfer function are known, then determining the poles and residues of the time history forũ 14 can be accomplished, (Table 9). From these corresponding poles and residues the time history u 14 (t) can be obtained, (Figure 21).

Comparison of Results
The displacement response x 14 due to an excitation at node 20 was found through the use of a frequency-domain, time-domain, and pole-residues methods.
The displacement response from these three methods do agree, (Figure 22). The time-domain and pole-residue methods agree strongly as they not only depict the steady-state portion of the response, but also the transient component.
The frequency-domain results depict just the steady-state portion.

Function Completion Times
The time to complete the execution of the above methods for the five-story building was determined through the use of the profile function in MatLab. The program was implemented twenty times and the execution time was captured for each method, (Figure 23). The time-domain method was the fastest with completion times just above the 0.05 seconds mark. The frequency-domain method was second and had completion times between 0.15 and 0.20 seconds. The pole-residue method was last with completion times between the 0.30 and 0.40 seconds.

CHAPTER 6 Conclusions
The displacement response was calculated through the use of a pole-residue method, frequency-domain method, and a time-domain method. There were two structures studied in this thesis. The first structure was a six degree-of-freedom system and the second was a one hundred and twenty degree-of-freedom system.
The calculated displacement response from these three methods were compared.
The first structure studied represented one half of a railway vehicle. The railway car was represented as a rigid body with two degrees of freedom, vertical displacement, and pitch angle. The railway car was connected to two nodes located at the front and rear end of the car representing the secondary suspension. These nodes were modelled as a rigid body with a single degree of freedom comprised of a vertical displacement. The secondary suspension nodes were connected through the primary suspension to a wheelset modelled as an unsprung mass with a single degree of freedom, comprised of a vertical displacement. The interaction between the wheel set and the track was modelled by linearized springs. The track was assumed to be fixed and rigid.
The second structure studied was a five story building. The five story building was represented as a three dimensional frame element with six degrees of freedom.
The building suffered from axial motion, transverse motion, and torsion. The axial motion was defined by the structure's cross sectional area, elastic modulus, and mass density. The transverse motion was governed by the structure's mass density, elastic modulus, cross sectional area, and moment of inertia. This structure underwent torsion which was governed by the shear modulus of elasticity, elasticity, and polar moment of inertia. The linkages were modelled as a three dimensional frame element with twelve degree of freedoms, six displacements, and six moments.
A frequency-domain method was implemented through the use of Fourier transform. The excitation was transferred from the time domain to the frequency domain through the fast Fourier transform. The complex frequency response was determined from the system's modes which were found by implementing a statespace model to simplify the problem. The complex frequency response was multiplied with the complex periodic excitation, and then shifted back to the timedomain through the inverse fourier transform. The resulting displacement response was found to be steady-state and lacks the transient compnents.
A time-domain method was implemented utilizing the functions ss and lsim found within MatLab. The ss command transferred the system into the required state-space model. The lsim command simulated the displacement response of continuous or discrete linear systems to arbitrary inputs through linear interpolation.
The displacement response was calculated through the implementation of a pole-residue method. The poles and residues associated with the excitation were determined from the implementation of a modified Prony's method. The poles and residues for the system were found using eigen analysis of the state-space model. These excitation and system poles and residues were then used to construct the displacement response.
The different displacements from the frequency-domain, time-domain, and pole-residue methods were compared. The displacement responses from these three methods are shown in Figures 15 and 22. From the initial inspection, one can say that three methods are in agreement. The frequency method provides just the steady-state response, and does not inlcude the transient portion. The timedomain and pole-residue method provided the complete displacement response containing the steady-state and transient portion.
The time required to execute the respective methods was accumulated over a large number of executions in order to determine execution speed. This was measured from the first step of rearranging the M, C, and K matrices to the required state-space form to the calculation of the displacement response. The execution time for each method to complete twenty iterations is shown in Figures   16 and 23. From the initial inspection, one can see that the frequency-domain method is the quickest when operating within a simple model. The analytical method was the fastest when operating on the more complex five-story building.
The pole-residue method was the slowest in both situations.
From the results provided in this paper, the built in analytical functions provided by MatLab performed the best. The pole-residue method provided the complete response, however it was the slowest. The frequency-domain method was the fastest for the simple six degree of freedom structure, but it provided just the steady-state response.
There are still several opportunities available to demonstrate the capabilities of this pole-residue method, especially since the work presented has been based on theoretical structures. The pole-residue method is an input output function and requires knowledge of the input and system. The utilization of this method in signal processing may allow for gains in signal detection and isolation. The inclusion of this algorithm in the detection of signals in the undersea environment may allow for gains in the study of anthropogenic noise and the study of marine mammals. The accuracy of the method may allow for the replacement of timedomain methods and allow for the quick iteration through material choices.