Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Dynamic analysis and control of one-dimensional distributed parameter systems
(USC Thesis Other)
Dynamic analysis and control of one-dimensional distributed parameter systems
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
DYNAMIC ANALYSIS AND CONTROL OF ONE-DIMENSIONAL
DISTRIBUTED PARAMETER SYSTEMS
by
Kyoungrae Noh
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)
August 2012
Copyright 2012 Kyoungrae Noh
ii
Acknowledgements
I would like to express my deepest thanks to my academic advisor, Dr. Bingen Yang,
for his guidance and encouragement throughout my graduate studies. I am also grateful to
the members of my Ph. D dissertation committee: Dr. H. Flashner, Dr. G. Shiflett and Dr.
L. C. Wellford for their constructive suggestions and advice. I would like to thank the
staffs, Sam and April, of Dept. of Aerospace and Mechanical engineering in USC for
providing a congenial research environment.
I am forever indebted my parents for giving me the opportunity and support to
accomplish my study. My work is theirs. Special thanks to my wife, Woojung, who was
there with love, trust and help.
iii
Table of Contents
Acknowledgements …………..………………………..…………………………….. ii
List of Figures ….…………..………………………..…………………………….. vii
List of Tables .........…………..………………………..…………………………….. x
Abstract …..……..…………..………………………..…………………………….. xi
Chapter 1 Introduction …....……………………………………………………… 1
1.1 Motivation and Purpose ..…………………………………..………. 1
1.2 Literature and Background ………..…..………………………….... 2
1.3 Summary …………..……..…………………..……………………. 10
Chapter 2 Fundamentals …................................................................... ……..…… 13
2.1 Introduction ……….…………..………………………..…………. 13
2.2 Single-Body Distributed Parameter Systems ……..…………..……14
2.2.1 Distributed Transfer Function Formulation …..………..…… 14
2.2.2 The s-domain Solution ……..……….……..…………..…… 15
2.2.3 Eigenvalue Problem ………..……..……..……………..…… 17
2.2.4 Frequency Response ……………..……..……..………..…… 18
2.2.5 Transient Response ……………..……..……..……. …..…… 18
2.3 Determination of Green’s Function …………..………..……..…… 19
2.3.1 Evaluation of State Transition Matrix …………….……..…… 19
2.3.2 Determination of Transfer Function Residues …………..…… 24
2.4 Multi-Body Distributed Parameter Systems ……..……..……..…… 28
2.4.1 Distributed Transfer Function Formulation …………....…… 28
2.4.2 Eigenvalue Problems ………..…..……..……..……..…..…… 32
2.4.3 Applications …………………..……..……..……..……..…… 33
iv
2.5 Mid- and High-Frequency Analyses ………..…..……..……..…… 33
2.5.1 Eigensolutions ………..…..……..……..……..……..…..…… 34
2.5.2 Transient Response …………..……..……..……..……..…… 35
Chapter 3 Vibration Analysis and Control of Bars, Shafts and Strings …..…… 43
3.1 Introduction …………..……..……..……..……..……..……..…… 43
3.2 System Description ………..……....……..……..……..……..…… 44
3.3 Distributed Transfer Function Formulation ………....………..…… 46
3.4 Dynamic Analysis ………..……...……..……..……..………..…… 47
3.4.1 The s-domain Solution ………..……..……..…………..…… 47
3.4.2 Eigenvalue Problem ……..…..……..……..……..……..…… 48
3.4.3 Frequency Response …….…..……..……..……..……..…… 49
3.4.4 Transient Response ………..……..……..……..………..…… 49
3.5 Specific Distributed Systems …………..……..……..………..…… 51
3.6 Boundary Control ………..…..…..……..……..……..………..…… 56
3.6.1 Boundary Controller ………………………………………….. 56
3.6.2 Stability of Root-Locus Analysis …………………………..… 57
3.7 Numerical examples ………....……..……..……..……..……..…… 60
Chapter 4 Analysis and Control of Acoustic Fields in a Duct ……........……..… 72
4.1 Introduction ……..………..……..……..……..……..………..…… 72
4.2 System Description ……..…..…..……..……..……..………..…… 72
4.3 Distributed Transfer Function Formulation ……....…………..…… 74
4.4 Dynamic Analysis ……..……………..……..……..…………..…… 75
4.5 Specific Duct Systems ……..………..……..……..…………..…… 78
4.6 Acoustic Control of Duct Systems ……...……..…….………..…… 81
4.6.1 Closed-Loop System …..…..…..……..……..…………..…… 81
4.6.2 Time-Delay Relations and Stability …..…….…………..…… 82
4.6.3 Perfect Noise Cancelling Controller …..…….…………..…… 85
4.7 Examples …..……..……..……..……..……..…….…………..…… 89
v
Chapter 5 Dynamic Analysis of Transversely Vibrating Beams …....…....…… 101
5.1 Introduction …..……..……...……..……..……..……………..…… 101
5.2 System Descriptions …..……....……..……..…….…………..…… 101
5.2.1 Euler-Bernoulli Beam Model …..……..…….…………..…… 101
5.2.2 Timoshenko Beam Model …..……..……..……………..…… 102
5.3 Distributed Transfer Function Formulation …...……………..…… 105
5.3.1 Euler-Bernoulli Beam ……...……..……..……………..…… 105
5.3.2 Timoshenko Beam …..……..……..……..……………..…… 106
5.4 Dynamic Analysis …..……..……..……..……..……………..…… 108
5.5 Spinning Timoshenko Beam …..………..……..……………..…… 117
5.5.1 System Description .…..…………………………………… 117
5.5.2 Dynamic Analysis .……..………………………………… 118
5.6 Timoshenko Beam Column ….………...……..……………..…… 125
5.6.1 System Description ..…..…………………………………… 125
5.6.2 Dynamic Analysis ..…..…………………………………… 126
5.7 Axially Moving and Loaded Timoshenko Beam ……………..…… 131
5.7.1 System Description .….…………………………………… 131
5.7.2 Dynamic Analysis .…...…………………………………… 132
Chapter 6 Dynamic Analysis of Stepped Bars, Shafts and Strings ………..…… 135
6.1 Introduction …..……..……..………..……..…….…………..…… 135
6.2 System Description …..….…..……..……..…….…………..…… 135
6.3 Distributed Transfer Function Formulation …..……………..…… 137
6.4 Dynamic Analysis …………....……..……..……..…………..…… 140
6.4.1 Eigensolutions ………..…. …..…..……..……..……...…… 140
6.4.2 Computation of Eigenvalues …...…………. …..……...…… 142
6.4.3 Transient Response ……………. …..……..……..……..…… 144
6.5 Example ………………...……..……..……..……..……..….…… 146
vi
Chapter 7 Multi-Body Structures …..…… ….……..……..…..……...…….…… 154
7.1 Introduction …………….…...……..……..……..……..…….…… 154
7.2 Stepped Beams ……………..…..……..……..……..……...…… 155
7.2.1 System Description ………………………………………… 155
7.2.2 Distributed Transfer Function Formulation ………………… 158
7.2.3 Applications . ……………………………………………… 160
7.3 Pointwise Connected Systems …..………….....……..……. .…… 167
7.4 Distributed Connected Systems ………………..……..…….…… 171
7.5 Combined Systems ……………….……..……..……..…….…… 177
7.6 Truss/Frame Structures ………………....……..……..……...…… 186
Chapter 8 Conclusions ……………..….…..……..……...….…..……..…….…… 194
References ...…………………..……..……..……..……..……..……..………..…… 198
Appendix A: Derivation of Solutions in Section 3.5 ………...……..……….…… 207
Appendix B: Proof of Eq. (4.36) …..……..…..……..……..…….…………..…… 209
Appendix C: Proof of Eq. (4.37) …..……..…..……..……..…….…………..…… 211
Appendix D: Proof of Eq. (4.39) ……....……..……..……..…….…………..…… 213
Appendix E: Proof of Eq. (4.40) …..…..……..……..……..…….…………..…… 214
Appendix F: Derivation of Eq. (4.62) ……..……..……..……..……………..…… 215
vii
List of Figures
Figure 2.1 Longitudinal bars connected with a spring .………………………….... 30
Figure 2.2 A stepped beam …………………………………...……………………... 39
Figure 2.3 A non-uniform beam with tip mass, damper and spring .………………... 40
Figure 2.4 The frequency response of the beam at x = 0.5 ………………………... 41
Figure 2.5 The transient displacement of the beam at x = 0.5 ..…………………….. 42
Figure 3.1 Schematic of a non-uniformly distributed system ....…………………… 45
Figure 3.2 Schematic of a controlled non-uniform system …....…………………… 59
Figure 3.3 A shaft with end disk in torsional vibration ……....…………………… 63
Figure 3.4 Spatial profiles of the forced vibration of the shaft ..…………………… 64
Figure 3.5 The forced vibration of the shaft at x = L ………….…………………… 65
Figure 3.6 The transient responses ( ,0.1) x of the shafts under the external
excitation
01.
( ) sin3355.4 ( )
t
t e t x x
with / 0.5 xL
...…………. 66
Figure 3.7 The eigenvalue loci of the bar with a controller (0 ≤ μ ≤ 2) ..…...……… 69
Figure 3.8 The frequency response of the bar at the right end (x / L=0.8) ……. …… 70
Figure 3.9 Spatial profiles of the forced vibration of the bar ..……………...……… 71
Figure 4.1 Schematic of noncollocated control of ducts ..…………………...……… 88
Figure 4.2 Relative acoustic pressure amplitude at the sixth eigenfrequency …….… 97
Figure 4.3 Transient waves in the duct due to a boundary excitation at the left end ... 98
Figure 4.4 Spatial distributions of the acoustic pressure due to
internal sound generation ………….……....………………….………… 99
viii
Figure 4.5 Time history for internal sound generation at the closed end (x = 0) ……. 99
Figure 4.6 Root loci of the duct under the noncollocated controller ..……………… 100
Figure 4.7 The frequency response of a duct with linear temperature profile ……. 100
Figure 5.1 Schematic of a beam model ………………..…………………...……... 105
Figure 5.2 The first natural frequency curves of the two models of a beam …….… 113
Figure 5.3 The first four mode shapes of beams ..………………………………… 114
Figure 5.4 Spatial profiles of the forced vibration of the Timoshenko beam ….…… 115
Figure 5.5 The forced vibration of the beams at x = L ..……………………...…… 116
Figure 5.6 A spinning beam and its inertial frame oxyz ..…………………….…… 117
Figure 5.7 Loci of the first 3 eigenvalues of the spinning beam with a varying
spin rate Ω ……………………………………………..…..…….…… 122
Figure 5.8 Spatial profiles of forced vibration of the spinning Timoshenko beam … 123
Figure 5.9 The time history (0 1) t of displacement (at x = L) of the spinning
Timoshenko beam with a (forward) spin rate Ω ………………….…… 124
Figure 5.10 Schematic of a beam column ………………....…………………..…… 125
Figure 5.11 Loci of the first 3 natural frequencies of the Timoshenko beam column
with a varying axial force ………………………..……………….…… 129
Figure 5.12 The time history (0 2) t of displacement (at x / L = 0.8) of
the Timoshenko beam column for an axial force N
x
…………...……… 129
Figure 5.13 Spatial profiles of the forced vibration of the Timoshenko beam column 130
Figure 5.14 The time history (0 2) t of displacement (at x / L = 0.8) of the axially
moving Timoshenko beam with a speed, c …………..……………….. 134
Figure 6.1 Schematic of an n-component stepped distributed dynamic system ……. 137
Figure 6.2 (a) Characteristic function of the stepped shaft from ref. (Yang, 2010)
(b) its derivative of the characteristic function with respect to ω …….... 144
ix
Figure 6.3 A three-segment bar in longitudinal vibration ………………....……… 150
Figure 6.4 Mode shapes of the stepped bar ………………..………………..…… 152
Figure 6.5 Spatial profiles of the transient displacement of the bar subject to
a sinusoidal boundary excitation
0
( ) sin
b
z t Z t ….…………….…… 153
Figure 7.1 Specification of constraints for stepped beams (Yang, 2005) ………….. 164
Figure 7.2 A stepped beam with hinged constraints …………………………….. 165
Figure 7.3 Mode shapes of the stepped beam …….…………………………..….. 165
Figure 7.4 The frequency response of the stepped beam at x = L …….……………. 166
Figure 7.5 A coupled beam structure …..…….………………………………….. 170
Figure 7.6 A distributed connected system …..…….…………………………….. 171
Figure 7.7 An elastically connected beam …..…….…………………………….. 175
Figure 7.8 Loci of the natural frequencies versus the change of the beam length l 176
Figure 7.9 Loci of the natural frequencies (the extended view in boxes in Fig. 7.8). 176
Figure 7.10 Connection of a distributed subsystem and a lumped parameter system ...177
Figure 7.11 A beam connected with an oscillator …………………..…….………. 181
Figure 7.12 Loci of the eigenvalues versus the coefficient of a spring k …………… 182
Figure 7.13 A beam with elastically mounted rigid body …………..…….…….…. 185
Figure 7.14 Schematic of a frame member …..…….…………………………….. 186
Figure 7.15 A frame structure …..……..……..………………………………..….. 192
Figure 7.16 Mode shapes of the frame structure …....…….…………………….…. 193
x
List of Tables
Table 2.1 The natural frequencies
k
of the non-uniform Euler beam ….………….. 27
Table 2.2 The natural frequencies
k
of the non-uniform Timoshenko beam ….….. 27
Table 2.3 The natural frequencies
k
of the stepped beam ..………..…………… 39
Table 2.4 The first ten eigenvalues of the beam …………………………………… 40
Table 3.1 The natural frequencies
k
of the non-uniform shaft ……….…………… 63
Table 3.2 The eigenvalues s
k
of the non-uniform bar ………....…………………… 69
Table 4.1 The eigenfrequencies
k
of the duct with three temperature distributions .. 97
Table 5.1 The natural frequencies
k
of the Timoshenko beam …..……………… 113
Table 5.2 The natural frequencies
kk
sj of the spinning Timoshenko beam ..… 122
Table 5.3 The eigenvalues
kk
sj of the moving Timoshenko beam ..………… 134
Table 6.1 The natural frequencies
k
of the stepped bar …….…………………… 151
Table 7.1 The natural frequencies of the stepped beam ..………………………… 166
Table 7.2 The natural frequencies of the coupled beam ..………………………… 170
Table 7.3 The natural frequencies of the elastically connected beam ..…………… 175
Table 7.4 The natural frequencies of the beam with an oscillator ..……………… 182
Table 7.5 The natural frequencies of the beam mounted a rigid body ..…………… 185
Table 7.6 The natural frequencies of the frame structure …...…………………… 192
xi
Abstract
An analytical method, based on the distributed transfer function method (DTFM), is
generalized and developed for modeling, analysis and control of one-dimensional
distributed parameter systems. Distributed parameter systems encompass a variety of
applications in wide engineering fields. Aerospace structures, building, bridges, turbines,
rotating machines and active sound suppression are examples of the application of
structural dynamics and control.
Distributed parameter systems studied in this work are classified as a single-body
distributed systems and a multi-body (complex) distributed system. A single-body
distributed system contains one continuum with boundary constraints; otherwise, a multi-
body distributed system is composed of multiple continua combined with lumped
parameter systems. The distributed transfer function formulation is generalized for both
types of the systems. Through use of state transition matrix, the formulation covers the
utility of non-uniformly distributed systems. By the formulation and the residue formula
for inverse Laplace transform, this method gives the exact transient response of the
system subject to arbitrary external, boundary and initial excitations. For case where the
exact solutions of non-uniformly distributed components are not available, a numerical
procedure is utilized to evaluate the state transition matrix with enough precision.
Eigenvalue distribution of stepped systems is investigated, and then a new algorithm is
introduced for determination of system eigenvalues.
xii
Active noise control of a duct with axial mean temperature gradients is considered in
this study. Through the introduction of specific time delays, the destabilizing effects of
non-collocation of sensors and actuators can be eliminated. By the DTFM, the exact and
closed-form solutions of the closed-loop duct system under delayed feedback controllers
is provided. The efficiency of delayed feedback control method is verified in numerical
examples.
The proposed method is engaged to analyze the natural frequency, mode shape,
frequency response, transient response and design of a controller of a distributed system.
This method does not require the system eigenfunctions. In the DTFM-based analysis,
system configurations such as the type of the continuous system, boundary conditions,
constraints and combined subsystems are treated in a systematic and unified manner. As a
result, the method is efficient in numerical computations and convenient in computer
coding.
1
Chapter 1. Introduction
1.1 Motivation and Purpose
Distributed parameter systems describe many important physical processes. Such
systems encompass a variety of applications in wide engineering fields such as aerospace
structures, building, bridges, turbines, rotating machines and active sound suppression.
The modeling of distributed parameter systems is described by partial differential
equations to be satisfied over the domain of the system, is subject to boundary conditions
at the end points of the domain.
In both time and frequency domain, complete solutions of a general distributed
parameter system has not been well addressed yet. Although modal analysis in theory can
give exact vibration solutions for simple distributed systems under external forces, it
requires orthogonal eigenfunctions, and has to deal with spatial integrals of system
eigenfunctions, which can be tedious in computation. In eigenfunction expansion,
handling general boundary excitations (i.e., inhomogeneous terms in boundary
conditions) is quite difficult, if not impossible. Furthermore, conventional analytical
methods require different derivations for different boundary conditions. Because of these
issues, analyses of distributed systems usually rely on numerical methods.
Although numerical methods are capable of delivering accurate results for complex
structures, exact analytical solutions are desirable for providing deep physical sight into
the problem and for validating the accuracy of numerical algorithms in development.
2
Additionally, from the viewpoint of vibration theory, exact analytical solutions for
flexible dynamic systems are always welcome.
The purpose of this study is to develop an analytical solution method for analysis of
distributed parameter systems, and design of a control system. Special attention is given
to development of the transient solution method for vibration problem of distributed
systems of which orthogonal eigenfunctions are difficult to achieve. The goals of the
investigation are of the following:
(1) To predict exact solutions for transient vibration of a distributed parameter system
with/without various non-self-adjoint effects, eigenvalue-dependent boundaries, and
under arbitrary external, initial and boundary disturbances
(2) To formulate the distributed transfer functions for non-uniformly distributed systems.
(3) To design stabilizing controllers for non-uniform systems
(4) To establish a technique to alleviate the risk of extensive computation rates and
missing eigenvalues in determining eigenvalues of a stepped system.
(5) To formulate the distributed transfer functions for complex distributed systems.
1.2 Literature and Background
Dynamic response prediction is important in modeling and analysis of flexible
vibrating systems, such as buildings, bridges, turbines, machines and aerospace structures
(Rayleigh, 1945; Timoshenko, 1981; Morse, 1948; Caughey and O’Kelly, 1965;
Meirovitch, 1967). Due to operational requirements and performance specifications,
flexible vibrating systems in many applications have spatially non-uniform distributions
3
of geometric and material properties and they are combined with lumped masses and
energy dissipative devices such as viscous dampers. These design considerations
mathematically lead to non-uniformly distributed systems that are non-proportionally
damped (Caughey and O’Kelly, 1965; Meirovitch, 1967). In vibration analysis of this
type of a distributed system, closed-form solutions in general are difficult to obtain;
approximate methods are commonly used. These flexible systems are the basic element
of structures and machines in a variety of engineering applications. Previous
investigations on non-uniform bars and shafts were mainly focused on free vibrations
relating to self-adjoint eigenvalue problems. Pouyet and Lataillade (1981) examined the
eigensolutions of a shaft with a power distribution of polar moment and inertia by
introducing a new variable through integral. Eisenberger (1991) obtained exact natural
frequencies of bars with polynomial variation in the cross-sectional area and mass
distribution. Bapat (1995) obtained exact natural frequencies and mode shapes of bars
with cross section of exponential and catenoidal distributions. Abrate (1995) determined
exact eigensolutions of bars with cross section of conical variation. Li, et al. (1996)
obtained free vibration solutions of undamped bars of cross-sections with polynomial and
exponential distributions. Kumar and Sujith (1997) obtained eigensolutions of bars with
cross section being a power function or squared sine function. Through use of coordinate
transformation, Li (2000, 2003), devised a formula for exact eigensolutions of certain
types of onestep, multi-step bars whose linear density and longitudinal rigidity are not
proportional to each other. Li, et al. (1998) presented exact eigensolutions for certain
4
types of proportionally damped bars. These results have helped better understand
dynamic behaviors of non-uniform systems.
While eigenvalue problems of non-uniform bars and shafts have been intensively
studied, closed-form solutions of forced vibration of these distributed systems with non-
proportional damping have not been well addressed. Even for uniform vibrating systems,
such as beams and bars, classical methods face difficulties in determination of closed-
form solutions in time domain, as discussed in References (Wu, 1998; Yang, 2008). Modal
expansion and Laplace transform method are two commonly-used methods in structural
dynamics. Yet these methods have not been shown to deliver closed-form solutions of
transient and frequency responses for general non-uniform and non-proportionally
damped vibrating systems. In the literature, no work has been reported on closed-form
solutions of transient and frequency responses of non-uniform and non-proportionally
damped systems by Laplace transform.
The modeling and analysis of acoustic ducts with temperature gradients are problems
of considerable interest in many engineering applications and thus has been extensively
studied. Due to design considerations, these duct systems often have non-uniform
temperature distributions, which are known to be of importance for the sound
propagation in duct systems. Examples are diverse, including automotive mufflers,
propulsion and power generating systems. Understanding of acoustic behaviors in such
duct systems subject to a variety of loading conditions is essentially important in a variety
of applications. Previous researches on analytical solutions have been mainly focused on
ducts without sound generation. Sujith, et al. (1995) and Kumar and Sujith (1997, 1998)
5
obtained exact solutions for ducts with several types of mean temperature profiles by
using the suitable transformation. Sujith (1996, 2001) also proposed a transfer matrix
formulation and obtained exact solutions for sound propagation with oscillating heat
release. Raspet, et al. (1993) derived analytical solutions for a duct with a linear
temperature gradient. Karthik, et al. (2000) obtained exact solutions for sound
propagation in ducts with mean temperature gradients and mean flow. Munjal and Prasad
(1986) and Peat (1988) developed exact solutions for ducts with small temperature
gradients. Peat (1997) presented a transfer matrix formulation for a duct with a linear
temperature gradient. Subrahmanyam, et al. (2001, 2003) produced exact analytical
solutions for wave propagation in ducts with mean temperature and area variations.
Besides analytical methods, approximate methods have also been used in analysis of duct
systems with or without sources. Dokumaci (1998) formulated an exact transfer matrix
and obtained approximate analytical solution for transmission of plane sound waves with
temperature gradients and area variations. Through use of Runge-Kutta method, Kapur, et
al. (1972)
obtained numerical solutions for sound propagation in ducts with axial
temperature gradients. Cummings (1977, 1978) developed an approximate method for
sound generation in ducts with temperature gradients. These studies have greatly
enhanced the understanding of acoustic behavior of flow duct systems with various
effects such as mean flows, area variations and temperature gradients. While much
research results have been reported, determination of analytical transient response of
ducts in the presence of sound generations has not been well addressed in the literature. It
is well known that analysis of ducts with temperature gradients has to deal with an
6
acoustic wave equation with spatially varying coefficients, whose exact solutions for
arbitrary temperature profiles in general are not available (Kumar and Sujith, 1997).
Conventional methods have difficulties in determination of transient solutions. Although
numerical methods are capable of delivering accurate results, exact analytical solutions
are desirable for better understanding the acoustic behaviors of one-dimensional wave
propagation in ducts with axial temperature gradients (Sujith, et al., 1995; Kumar and
Sujith, 1997; Kumar and Sujith, 1998).
The active noise control has been received interests since Lueg (1936) first proposed
the concept, and extensively studied. Recent advances in electronics and microprocessor
have made this concept be practical and economical. An active noise control system in a
duct usually has one or more control speakers and microphones at various measurement
locations. A variety of investigations for active sound attenuation have been presented in
several decades (Swinbank, 1973; Ross, 1981; Trindler and Nelson, 1983; LaFontaine
and Shepherd, 1983; Eriksson, et al, 1987). The concept of acoustic feedback is raised
with the attempt to isolate microphones’ measurement from control speakers (LaFontaine
and Shepherd, 1983; Hu, 1995). Hull, et al. (1990) developed a feedback control
technique for a discretized acoustic model of a one-dimensional duct. Hu (1995)
investigated the problem of active sound cancellation in a duct. Most research of active
noise control has been focused on a duct system with a uniform temperature distribution.
In design of a controller, non-collocation of sensors and actuators generally results in loss
of controllability and stability of the system due to spillover effects. Yang and Mote
(1991) introduced that such destabilizing effect can be alleviated through the introduction
7
of a specific time delay block(s) in the control system and implemented their approach
for active vibration control of an axially moving string and showed that all the vibrating
modes can be stabilized and the spillover instability can be avoided.
Han, et al. (1999) represented solutions of frequency equations and forced vibration of
four models of a beam using the method of eigenfunction expansion. He indicated that
there is no paper which presents the complete solution from the formulation of the
governing differential equation to its solution for all four beam models. Katz, el al.
(1988) employed a finite integral transform to calculate transient response for spinning
Rayleigh and Timoshenko beams with the simply supported boundary condition and
moving load. Zu and Han (1992, 1994) obtained closed-form solutions for free vibration
analysis of spinning Timoshenko beams for general boundary conditions and subjected to
a moving load using modal analysis technique. Tan and Kuang (1995) presented for
analyzing the vibration of a stepped rotating shaft with general boundary conditions. In
1980, Kounadis derived equations of motion for vibrating Timoshenko column.
Aristizabal-Ochoa (2004) investigated free vibration problems of Timoshenko beam
column with general boundary conditions. Based on Haringx theory, Aristizabal-Ochoa
(2008) studied the behavior of a compressed member of symmetrical cross-section with
semi-rigid connections. Blaauwendraad (2008) indicated that Ochoa’s theory appears to
be reliable for helical springs, but is inadequate for shear-weak-beam-columns. Lee, et al.
(2004) developed the spectral element method for vibration analysis of an axially moving
Timoshenko beam. These studies have enhanced the understanding of dynamic behaviors
of transversely vibrating beams with/without various effects. Determination of complete
8
solutions of all beam models subject arbitrary external, initial, and boundary conditions
has not been addressed.
Stepped second-order distributed parameter system is assemblages of one-dimensional
distributed component (like elastic bars, torsional shafts and taut strings), and lumped
masses. Vibration of stepped second-order distributed dynamic systems is important in
many engineering applications, such as buildings, bridges, rotating machines, turbines,
helicopters and aerospace structures, and hence has been extensively studied (Snowdon,
1968; Gorman, 1968; Dimarognas and Paipeties, 1983). Most investigations on this type
of vibrating systems have been focused on free vibration, forced response to sinusoidal
excitations and wave propagation (Munjal, et al., 1973; Mioduchowski, 1986; Bapat and
Bhutani, 1994; Yang, 1995; Hsueh, 1999; Li, 2003). Determination of transient vibration
of stepped distributed systems, however, has mainly relied on numerical methods,
although certain analytical solutions can be obtained by eigenfunction expansion (Yang,
1996; Fedotov, et al., 2006, 2008), and by a time-domain receptance method that is valid
for stepped systems with lumped parameter components (Li and Stone, 1999). The utility
of the existing analytical techniques in transient analysis of stepped systems is restricted
(Yang, 2010). Recently, Yang (2010) has presented compact distributed transfer function
formulation of stepped systems composed of uniform distributed components.
In structural dynamics analysis, the eigenvalue (natural frequency) contains important
information of a system. The natural frequencies of vibration of an infinite-dimensional
system can be determined by finding the roots of characteristic equation. By use of
dynamic stiffness matrix, the characteristic equation has a transcendental form depending
9
on the square of (circular) frequency ω (Wittrick and William, 1971; Lee, et al., 2002).
Wittrick and William (1971) developed an algorithm called ‘W-W algorithm’ to calculate
the number of natural frequencies that lie below any chosen frequency. Lee and his co-
workers (2002) developed W-W algorithm to determine natural frequencies of elastic-
piezoelectric two-layer beams by the spectral element model. The algorithm is applicable
and useful for determining all required natural frequencies from the characteristic
equation which may be obtained by dynamic stiffness matrix. However, the DTFM
provides the condensed characteristic equation for which W-W algorithm may not be
valid.
The dynamics and vibration of multi-body structures or complex distributed parameter
systems that composed of multiple distributed and lumped parameter subsystems have
been extensively studied by many researchers. Variety of methods have been developed
such as the transfer matrix method (Prestel and Leckie, 1963; Lin, 1969), wave
propagation analysis (Mead, 1971; Mead and Yaman, 1990), modal analysis (Chen, Y.,
1963; Bergman and Nicholson, 1985), the Green’s function method (Nicholson and
Bergman, 1986; Bergman and McFarland, 1988), and the distributed transfer function
formulation (Butkovisky, 1983; Yang and Tan, 1992). Instead of the rich literature on
complex distributed parameter systems, there still remain some research problems of
interest. Most works rely on discretization or truncated eigenfunction expansions. In
addition, many existing techniques require different derivation in different cases of
system type, constraints, boundary conditions and external disturbances. As such,
10
computer coding for general complex distributed parameter systems is difficult. The
distributed transfer function synthesis has been presented by Yang (1994).
1.3 Summary
The distributed transfer function method is briefly introduced in chapter 2. The general
one-dimensional distributed parameter systems in consideration are classified as single-
body distributed systems and multi-body distributed systems. Through the use of state
transition matrix, the method is generalized for non-uniformly distributed systems. The
DTFM method is engaged to analyze eigensolutions, frequency response, and transient
response of distributed systems subject to arbitrary initial, boundary and external
disturbances. By an augmented spatial state form and Jacobi’s formula, evaluation of
state transition matrix and transfer function residues is studied. For arbitrary non-
uniformly distributed systems (components) whose exact solutions are not available, a
numerical procedure can be utilized with enough precision for determination of state
transition matrix. For multi-body distributed systems, a new distributed transfer function
formulation are presented and investigated. The DTFM method is implemented for
analysis of the system in mid- and high-frequency regimes. In transient vibration analysis,
the method is utilized to obtain asymptotic transient solutions by the selection of modes
in Green’s functions.
In chapter 3, vibration analysis and boundary control of non-uniform systems under
initial, external and boundary excitations is discussed. The DTFM-based method gives
closed-form transient and frequency response for a class of non-uniformly distributed
11
systems, namely, bars in longitudinal vibration, shafts in torsional vibration, and taut
strings in transverse vibration. Boundary controllers which may contain lumped masses
and energy dissipative devices mathematically lead to non-uniformly distributed systems
that are non-proportionally damped. In vibration analysis of this type of a distributed
system, closed-form solutions in general are difficult to obtain; approximate methods are
commonly used. By the DTFM-based method, time and frequency-domain solutions are
determined in exact and closed-form.
Chapter 4 presents analytical solutions of transient response and active noise control of
one-dimensional duct systems with axial mean temperature gradients that are under
arbitrary sound generations, boundary disturbances and initial conditions. Design of
stabilizing controllers is carried out for sound attenuation of ducts with axial temperature
gradients. Two different types of boundary conditions of a duct are considered for design
of non-collocated controllers. Noncollocation of sensors and actuator can cause spill-over
instability. However, through introduction of specific time delays in the feedback
controller, the destabilizing effect caused by noncollocation of sensors and actuators is
eliminated. The root locus method is employed for the stability analysis. In addition, the
frequency response of the closed-loop system is presented.
In chapter 5, time-domain analysis of two models for beams with non-self-adjoint
effects is studied. Two beam models are the Euler-Bernoulli and Timoshenko. For the
Timoshenko beam model, the proposed method deals with transient analysis of a beam
with the effects such as a spin, an axial loading and an axial movement, and the method
delivers an exact and closed-form solution in time domain.
12
In chapter 6, vibration problem of stepped distributed dynamic systems which are
assemblages of non-uniformly distributed components (like non-uniform bars, shafts and
strings), and lumped masses are studied. An analytical method is developed for exact
solution of eigenvalue problem, frequency response, and transient vibration of a stepped
distributed system. One of the key steps of the proposed method is to evaluate the
Eigensolutions. In order to guarantee no eigenvalues to be missed, eigenvalue distribution
is analyzed and the new algorithm is introduced into the distributed transfer function
method as shall be seen in section 6.4.2.
In chapter 7, the new formulation is introduced for vibration problem of multi-body
structures based on the distributed transfer function method. Multi-body structures are
composed of distributed components such as bars, beams, trusses and frames, and some
lumped parameter subsystems. The formulation is explicitly justified in various complex
distributed parameter systems, i.e. point-wise connected beams, distributed connected
beams, combined beams, and truss/frame structures. The conclusion is represented in
chapter 8.
13
Chapter 2. Fundamentals
2.1 Introduction
This chapter is concerned with analytical method for modeling and analysis of general
one-dimensional distributed parameter systems. The distributed parameter system can be
considered as an assembly of distributed and (or) lumped parameter elements. For those
cases usually involve local discontinuity and non-uniformity such as a joint and a varying
cross-section. In such cases, it is difficult to obtain an exact and closed-form solution of
the system by conventional methods. In this research, an analytical method is developed
based on the Distributed Transfer Function Method (DTFM), which provides a closed-
form analytical solution for analysis modeling of solids and structures (Yang and Tan,
1992; Yang, 2005). In the development, a spatial state equation in the Laplace transform
domain is first devised. Through determination of a state transition matrix, the distributed
transfer function of the system is precisely evaluated; the s-domain solution of the state
equation is then obtained through use of the distributed transfer functions. With a residue
formula, inverse Laplace transform of the s-domain solution leads to exact and closed-
form transient solutions. The proposed method, which does not depend on system
eigenfunctions, is convenient in computer coding and numerically efficient in that the
method can avoid numerical instability conventional analytical methods generally
encounter in mid- and high-frequency ranges.
14
2.2 Single-Body Distributed Parameter Systems
2.2.1 Distributed transfer function formulation
Consider a one-dimensional continuum of length L, where ( , ) w x t is the generalized
response of the continuum. The analysis of the system in consideration leads to a set of n-
order partial differential equation (Yang, 1992).
2
2
0
( ) ( ) ( ) ( , ) ( , ), (0, )
k n
k k k k
k
a x b x c x w x t f x t x L
t t x
(2.1a)
The system is also subject to boundary and initial conditions
0
( , ) ( , ) ( ), 1,2, ,
j j j
x x L
M w x t N w x t t j n
(2.1b)
00
( ,0) ( ), ( ,0) ( ) w x u x w x v x
t
(2.1c)
where ( , ) f x t is the external load and ()
j
t are the boundary disturbances;
j
M and
j
N
are temporal-spatial, linear differential operators which describe the boundary conditions;
()
k
ax , ()
k
bx and ()
k
cx describe the spatial distributions of system parameters, such as
inertia, damping, gyroscopic forces, stiffness and axial loading;
0
() ux and
0
() vx are
given functions describing the initial states of the system.
Equations (2.1), non-self-adjoint in general, describe many vibration continua such as
bars, shafts, strings, beams, beam columns, flexible rotating shafts, and axially moving
materials [17-18].
15
First, the s-domain solution of the system is devised by the distributed transfer function
method (Yang and Tan, 1992; Yang, 2005). By taking Laplace transform to Eqs. (2.1)
with respect to time and defining the spatial state vector as
T
1
1
( , ) ( , )
ˆ ( , ) ( , )
n
n
w x s w x s
x s w x s
xx
(2.3)
Equations (2.1) are rewritten in the spatial state form
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ) x s x s x s x s
x
Fp (2.4a)
ˆ ˆ (0, ) ( , ) ( )
b b b
s L s s MN (2.4b)
where the overbar stands for Laplace transform, T is the transpose, and s is the Laplace
transform parameter; The matrix ( , ) xs F is composed of the coefficients of the differential
operators; The matrices
b
M
and
b
N
contain those of the boundary operators; ˆ( , ) xs p and
()
b
s are vectors to contain the external and the boundary disturbances with the initial
states, respectively.
2.2.2 The s-domain Solution
In this study, the state equation (2.4) is solved through use of state transition matrix.
The state transition matrix of the system is the unique solution of (Chen, 1999)
( , , ) ( , ) ( , , ), , (0, ) x s x s x s x L
x
F (2.5)
subject to ( , , ) x x s I , where I is the identity matrix. The state transition matrix has
the properties
16
1
( , , ) ( , , )
( , , ) ( , , ) ( , , )
x s x s
x z s x y s y z s
(2.6)
Furthermore, the state transition matrix can be written as
1
( , , ) ( , ) ( , ) x s x s s
UU (2.7)
where ( , ) xs U is any fundamental matrix which is a non-singular solution of
( , ) ( , ) ( , ), (0, ) x s x s x s x L
x
U F U (2.8)
With the state transition matrix, the s-domain solution of the distributed system is taken
as
0
ˆ ˆ
ˆ ˆ ( , ) ( , , ) ( , ) ( , ) ( ), , (0, )
L
b
x s x s x s d x s s x L
G p H (2.9)
where
ˆ
G and
ˆ
H are the distributed transfer functions of the system, and they are given
by
1
ˆ
( , ) (0, , ),
ˆ
( , , )
ˆ
( , ) ( , , ),
ˆ
( , ) ( ,0, ) ( )
b
b
x s s x
xs
x s L s x
x s x s s
HM
G
HN
HZ
(2.10)
with
( ) ( ,0, )
bb
s L s Z M N (2.11)
where ( , , ) xs is the state transition matrix of the system.
17
2.2.3 Eigenvalue problems
The proposed transient solution method requires the knowledge of system eigenvalues.
Define the eigenvalue problem of the system by vanishing the excitation terms in Eq.
(2.4):
()
( , ) ( ), (0, )
x
x s x x L
x
F
(2.12)
subject to the conditions
(0) ( ) 0
bb
L MN (2.13)
where s is an eigenvalue, and () x the eigenfunction associated with s. Assume that
( ) ( , , ) x x s a (2.14)
with a being a constant vector to be determined. The () x , which satisfies Eq. (2.12), is
plugged in the boundary condition (2.13), yields
( ,0, ) 0
bb
Ls M N a = (2.15)
From Eq. (2.15), the characteristic equation of the system is
det ( ) det ( ,0, ) 0
bb
s L s Z M N (2.16)
whose roots,
k
s , 1,2, , k
are the eigenvalues of the system. Although eigenfunctions
are not needed in the proposed transient analysis, the eigenfunction (mode shape)
associated with the kth eigenvalue
k
s can be obtained by
( ) ( , , )
kk
x x s a (2.17)
where
k
a is a nonzero solution of the homogenous equation ( ) 0
k
s Za .
18
2.2.4 Frequency response
Consider the system subject to harmonic (sinusoidal) excitations, ( , ) ( )
jt
x t x e
pp
and
0
()
jt
b
te
. By Eq. (2.9), the steady-state solution of the system is
0
0
ˆ ˆ
ˆ ( , ) ( , , ) ( ) ( , ) , (0, )
L
j t j t
x t x j d e x j e x L
G p H (2.18)
where 1 j , and is the excitation frequency. The
ˆ
( , , ) xj G and
ˆ
( , ) xj H are
known as frequency response function of the system. By the DTFM-based analysis (2.10),
the frequency response of the system is given in the exact and closed form.
2.2.5 Transient response
Inverse Laplace transform of Eq. (2.9) yields the Green’s function formula
0 0 0
( , ) ( , , ) ( , ) ( , ) ( )
t L t
b
x t x t d d x t d
G p H (2.19)
where ( , ) xt is the inverse Laplace transform of ˆ ( , ) xs ; and the matrix Green’s
function G
and H
are the inverse Laplace transform of the transfer functions
ˆ
G
and
ˆ
H ,
respectively.
By the theorem of residues (Morse, 1953; Yang, 2008, 2010), the Green’s functions
are determined as
1
1
( , , ) 2 Re ( ,0, ) (0, , )
( , ) 2 Re ( ,0, )
k
k
st
k k b k
k
st
kk
k
x t x s s e
x t x s e
G R M
HR
(2.20)
where Re{} is the real part of {} ;
k
s is the kth eigenvalue of the system, and
19
1
adj ( )
Res ( )
()
k
k
k
k
ss
ss
s
s
d
s
ds
Z
RZ
Z
(2.21)
where adj ( ) s Z and () s Z are the adjoint and determinant of Z(s), respectively. R
k
shall be
called the residues of the distributed transfer functions or transfer function residues. The
transient solution given by Eqs. (2.19) and (2.20) is in exact if the integrals can be
evaluated by the exact quadrature. In this solution process, neither approximation nor
system eigenfunctions has been used.
2.3 Determination of the Green ’s Function
The transient response of the distributed parameter systems, Eq. (2.19), is determined
by the Green’s function formulations. By the theorem of residue, the Green’s functions of
the system are given as Eqs. (2.20) and (2.21). The state transition matrix and transfer
function residues are required to determine the distributed transfer functions and the
Green’s functions. In the DTFM-based analysis, the key is to determine exact state
transition matrix and transfer function residues. Once the state transition matrix and
transfer function residues are known, the distributed transfer functions and the Green’s
functions are calculated, and various dynamic problems can be conveniently investigated.
2.3.1 Evaluation of State Transition Matrix
This section is concerned with evaluating the state transition matrix of non-uniformly
distributed parameter systems. Although it is always desirable to determine exact state
20
transition matrix, there are not always exact solutions yet. In other words, the state
transition matrix of a non-uniformly distributed parameter system, in general, cannot be
expressed in exact and closed-form and has to be evaluated numerically. Yang and Fang
(1994) presented the methods for the evaluation of state transition matrix: step-function
approximation and truncated Taylor series. Under the existence of exact solutions of
Equation (2.1a), a fundamental matrix is exactly obtained with non-singular solutions of
Eq. (2.8), and thus state transition matrix can be exactly evaluated by Eq. (2.7).
Consider the state transition matrix ( , , ) xs in Eq. (2.5). The state transition matrix
has the property ( , , ) ( , , ) ( , , ) x z s x y s y z s . Define the commutative condition as
00
( , ) ( , ) ( , ) ( , )
xx
xx
x s s d s d x s
F F F F for all and x (2.22)
If ( , ) xs F satisfies the commutative property ( 0) , then the solution of Eq. (2.5) is
written as (Chen, 1999)
0
( , )
0
( , , )
x
x
sd
x x s e
F
(2.23)
When ( , ) xs F is constant, i.e. ( , ) ( ) x s s FF , the state transition matrix is determined by
exponential matrix
( )( )
( , , )
sx
x s e
F
. Furthermore, the approximate state transition
matrix can be obtained by Eq. (2.23) when ε is relatively very small. Assume that the
system is composed of the sequence of subsystems. Define that
k
is the state transition
matrix of the kth subsystem. By the property, Eq. (2.6), the state transition matrix of the
system is given as
0 1 2 1
( , , )
n n n n
x x s
(2.24)
21
where
11
( , , ),
k k k k k
x x s x x
and 1,2, , kn .
Thereby, the state transition matrix is computed by
0 0 1 1
ˆ
( , , ) ( , , )
n n n n
x x s x x s
(2.25)
where
1
( , )
x
k
x
k
sd
j
e
F
. The accuracy of the matrix
ˆ
increases as the number of the
interval.
The distributed transfer functions can be computed with the known state transition
matrix. With the distributed transfer function formulation, various dynamic problems of
non-uniformly distributed parameter systems can be conveniently investigated. To prove
the efficiency, the method developed is illustrated on eigenvalue problems of two non-
uniform beams.
Example 2.1: Eigensolutions of a Nonuniform Euler-Bernoulli Beam
As Euler-Bernoulli beam model, the displacement ( , ) w x t of the system is governed by
2 2 2
2 2 2
( , ) ( , )
( ) ( ) 0, (0, )
w x t w x t
x EI x x L
t x x
(2.26)
where ρ(x) is the linear density, EI(x) the bending stiffness as functions of the spatial
coordinate x.
Let the parameters of the beam be
22
00
( ) (1 ) , ( ) (1 ) x ax EI x EI ax (2.27)
Assume that the beam is fixed at the left end (x = 0) and free at the right end (x = L). The
boundary conditions are assigned as
22
22
22
(0, )
Left end : (0, ) 0, 0
( , ) ( , )
Right end : ( ) 0, ( ) 0
xL
wt
wt
x
w L t w x t
EI L EI x
x x x
(2.28)
Eqs. (2.26) and (2.28) can be expressed by Eq. (2.4) and the state matrix is given by
2
2 0
2 2 2
2
0
0 1 0 0
0 0 1 0
0 0 0 1
( , )
(1 ) 2
0
(1 ) (1 )
(1 )
xs
ax a ax
s
ax ax
EI ax
F (2.29)
Choose the values of the system parameters be
00
0.5, 2, 0.2, 1 EI a L . The
system considered here is non-uniformly distributed. It is known that the exact and
closed-form of state transition matrix of the system is not available. The approximate
state transition matrix
ˆ
can be evaluated by Eq. (2.25). With the approximate state
transition matrix, the characteristics equation is expressed as
ˆ ˆ
det ( ) det ( ,0, ) 0
k b b k
j L j Z M N (2.30)
The natural frequencies , 1,2,...
k
k , satisfying Eq.(2.30), are listed in Table 2.1. The
numerical results by the proposed approximate method have good agreements with the
FEM results. As the number of elements increases, the accuracy of the proposed method
and FEM goes higher.
23
Example 2.2: Eigensolutions of a Nonuniform Timoshenko Beam
As a Timoshenko beam model, the displacement ( , ) w x t and the rotation ( , ) xt of the
system are governed by
2
2
2
2
( , ) ( , )
( ) ( ) ( , ) 0, (0, )
( , ) ( , ) ( , )
( ) ( ) ( ) ( , ) 0
s
s
w x t w x t
x k x x t x L
t x x
x t x t w x t
I x EI x k x x t
t x x x
(2.31)
where ρ(x) is the linear density, EI(x) the bending stiffness, I
ρ
(x) the rotary inertia and
k
s
(x) the shear stiffness as functions of the spatial coordinate x.
Assume that the beam is fixed at the left end (x = 0) and free at the right end (x = L).
The distributions of the system parameters are of exponential form:
2 / 2 / 2 / 2 /
0 0 0 0
( ) , ( ) , ( ) , ( )
x L x L x L x L
s
x e EI x EI e I x I e k x k e
(2.32)
where
0 0 0 0
, , and EI I k are constants.
Eq. (2.31) can be expressed by Eq. (2.4) and the state matrix is given by
2
0
0
2
0 0 0
00
0 0 1 0
0 0 0 1
22
1 ( , )
2
0
s
xs
k L L
I s k k
EI EI L
F (2.33)
In numerical simulation, the following non-dimensional values for the system
parameters are assigned as
4
0 0 0 0
1, 0.02, 10, 10 , I EI k L=1 and α = 0.2. The state
transition matrix of the system can be obtained in exact and closed-form because the state
matrix is constant, i.e. ( , ) ( ) x s s FF . Table 2.2 lists the natural frequencies of a non-
24
uniform Timoshenko beam which are computed by det ( ) 0
k
j Z . As the number of
elements increases, the FEM results approach the results by the proposed method.
2.3.2 Determination of Transfer Function Residues
This section represents evaluation of the transfer function residues by Eq. (2.21). The
denominator term of Eq. (2.21), by Jacobi’s formula, can be expressed as
()
( ) tr adj ( )
k
k
k
ss
ss
d d s
ss
ds ds
Z
ZZ (2.34)
where tr ( ( )) s Z is the trace of () s Z .
With Eq. (2.11), the derivative of () s Z with respect to s is given as
( ) ( , , )
( , , )
bb
b
dd d s d L s
Ls
ds ds ds ds
MN Z
N
(2.35)
With known the derivative of state transition matrix with respect to s, the derivative of
() s Z and () s Z can be computed by Eqs. (2.35) and (2.34), respectively. Then, the
transfer function residues are obtained and the Green’s functions can be evaluated by the
residue theorem. In the proposed method, it is required to calculate the derivative of state
transition matrix with respect to s.
First, devise an augmented state form. Differentiating Eq. (2.8) with respect to s,
( , ) ( , ) ( , )
( , ) ( , ) (0, )
d x s d x s d x s
x s x s x L
x ds ds ds
U F U
UF (2.36)
Eqs. (2.8) and (2.28) can be rewritten as the augmented state form
25
ˆ ˆ ( , ) ( , ) ( , )
a a a
x s x s x s
x
F (2.37a)
where
( , ) ( , )
ˆ ( , ) , ( , )
( , ) ( , )
( , )
aa
x s x s
x s x s
d x s d x s
xs
ds ds
U F 0
F
UF
F
(2.37b)
Here, 0 is a zero matrix. With the state transition matrix ( , , )
a
xs , the solution of Eq.
(2.37) can be expressed as
ˆ ˆ ( , ) ( , , ) ( , )
a a a
x s x s s . By applying the same procedure
to Eq. (2.7), it can be rewritten as
( , ) ( , , ) ( , )
( , ) ( , , ) ( , )
( , , )
x s x s s
d x s d x s d s
xs
ds ds ds
U 0 U
UU
(2.38)
By comparing it with
ˆ ˆ ( , ) ( , , ) ( , )
a a a
x s x s s , it is shown that the state transition
matrix ( , , )
a
xs is composed of the state transition matrix ( , , ) xs and its derivative
with respect to s.
( , , )
( , , )
( , , )
( , , )
a
xs
xs
d x s
xs
ds
0
(2.39)
The state transition matrix of the augmented form, Eq. (2.39), is the unique solution of
( , , ) ( , ) ( , , )
a a a
x s x s x s
x
F (2.40)
As explained in section 2.2.1, the solution of Eq. (2.40) is determined, and then the state
transition matrix ( , , ) xs and its derivative are simultaneously obtained as submatrices
26
in Eq. (2.39). With the given derivative of state transition matrix, the transfer function
residues are determined by
adj ( )
()
tr adj ( )
k
k
k
k
sj
j
ds
j
ds
Z
R
Z
Z
In the proposed transient analysis, exact solutions of the system can be obtained by the
transfer function residues.
27
Table 2.1 The natural frequencies
k
of the non-uniform Euler beam (in rad/s)
k
The propose method FEM
N
*
= 1 4 8 N
*
= 16 48
1
2
3
4
5
6
7
8
7.873
48.754
135.337
264.566
436.894
652.295
910.773
1212.327
8.309
49.188
135.179
265.190
437.507
652.302
913.053
1214.290
8.348
49.443
136.092
265.469
438.001
653.665
912.462
1214.323
8.349
49.453
136.125
265.564
438.267
654.377
914.210
1218.278
8.352
49.471
136.171
265.622
438.249
654.029
912.968
1215.075
N
*
: number of elements
Table 2.2 The natural frequencies
k
of the non-uniform Timoshenko beam (in rad/s)
k
The proposed
method
Finite element method
N
*
= 20 40 100
1
2
3
4
5
6
7
8
11.932
55.857
123.921
197.502
272.964
347.772
421.900
494.876
11.932
56.095
125.220
201.208
280.718
361.575
444.093
528.057
11.932
55.916
124.245
198.422
274.887
351.189
427.392
503.098
11.932
55.866
123.973
197.649
273.271
348.317
422.775
496.186
N
*
: number of elements
28
2.4 Multi-Body Distributed Parameter Systems
2.4.1 Distributed transfer function formulation
Multi-body distributed parameter systems are composed of multiple continua and
lumped parameter subsystems. Each continuum may be under various constraints,
combined with lumped parameter systems, and subject to external and boundary
disturbances. The analysis and modeling of the multi-body distributed parameter system
in consideration leads to a set of governing equation along with boundary and initial
conditions. In addition, matching conditions are required at the nodes where the
subsystems are interconnected. The multi-body distributed system in consideration is an
assembly of N one-dimensional subsystems. Two nodes,
, si
x and
,
,
ei
x are representing the
ends of the sub-system and
,, i e i s i
l x x is the length of the ith subsystem. Let the sub-
domain be
i
,,
( , )
s i e i
xx . The response ( , )
i
w x t of the ith subsystem is governed by the
nth order, linear partial differential equation (Yang, 1994)
2
, , , 2
0
( ) ( ) ( ) ( , ) ( , ),
k n
k i k i k i i i i k
k
a x b x c x w x t f x t x
t t x
(2.41a)
0, 0,
( ,0) ( ), ( ,0) ( ), 1,2, ,
i i i i
w x u x w x v x i N
t
(2.41b)
where ( , )
i
f x t is the external disturbance;
,,
( ), ( )
k i k i
a x b x and
,
()
ki
cx describe the spatial
distributions of subsystem parameters; and
0,
()
i
ux and
0,
()
i
vx are given profiles of initial
states of a subsystem.
The system is subject to the boundary conditions and the matching conditions
29
,,
( , ) ( , ) ( )
s i e i
ij i ij i ij
x x x x
M w x t N w x t t
(2.42)
where
ij
M
and
ij
N
are linear temporal-spatial differential operators of proper order and
()
ij
t is related to a boundary disturbance and an external load applied at the nodes where
the subsystems are interconnected.
By the distributed transfer function formulation (Yang, 2005), define the state vector
1
1
ˆ ( , ) , 1,2, ,
T
n
ii
ii n
ww
x s w i N
xx
(2.43)
The subsystem equations can be expressed in state forms
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ),
i i i i i
x s x s x s x s x
x
Fp (2.44)
where
ˆ ( , )
i
xs p is related with the external disturbance and the initial conditions of
subsystems.
For convenience of analysis, define a global domain
12 N
and
convert to Eq. (2.44) to a global state equation
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ), x s x s x s x s x
x
Fp (2.45a)
where
12
12
12
ˆ ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
ˆ ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
( , ) ( , ) ( , ) ( , )
T
T T T
N
T
T T T
N
N
x s x s x s x s
x s x s x s x s
x s diag x s x s x s
p p p p
F F F F
(2.45b)
The boundary and matching conditions, Eq. (2.42), are rewritten in the augmented form
ˆ ˆ ˆ ( , ) ( , ) ( )
G s G e G
x s x s s MN (2.46a)
30
where
1 ,1 2 ,2 ,
1 ,1 2 ,2 ,
ˆ ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
ˆ ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
T
T T T
s s s N s N
T
T T T
e e e N e N
x s x s x s x s
x s x s x s x s
(2.46b)
Here, the matrices
G
M and
G
N are composed of the coefficients of the operators
ij
M and
ij
N , and
ˆ ()
G
s is described with the initial conditions, the boundary disturbances and
external load at the internodes. Note that boundary and matching conditions are classified
clearly to formulate a global state form, Eq. (2.46).
x
s,1
x
e,1
x
s,2
x
e,2
ρ
1
, EA
1
ρ
2
, EA
2
k
l
1
l
2
Figure 2.1 Longitudinal bars connected with a spring
As an example, consider longitudinal bars connected with a spring; see Fig. 2.1. The
displacement of the bars is governed by the differential equation
22
22
( , ) ( , ) ( , ), 1,2
i i i i i
u x t EA u x t f x t i
tx
(2.47)
The system is subject to the boundary conditions
,1
,2
1
1 1 1 ,1 2 ,2 1 ,1
2
2 2 2 ,2 1 ,1 2 ,2
0, ( , ) ( , ) ( , )
0, ( , ) ( , ) ( , )
s
e
e s e
xx
s e s
xx
w
EA EA w x t kw x t kw x t
xx
w
EA EA w x t kw x t kw x t
xx
(2.48)
31
Assume that the initial conditions are zero. Taking Laplace transform to Eqs. (2.47)
and (2.48), the state equations of the system are assigned as
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
ˆ ˆ ˆ ( , ) ( , ) ( )
G s G e G
x s x s x s x s
x
x s x s s
Fp
MN
(2.49a)
where
2
ˆ ( , ) { / } , 1,2
01
( , )
/0
T
i i i
i
ii
x s w w x i
xs
s EA
F
(2.49b)
12
12
ˆ ˆ ˆ ( , ) ( , ) ( , )
( , ) ( , ) ( , )
T
TT
x s x s x s
x s diag x s x s
F F F
(2.49c)
1
22
1
2
22
2
00 0 0 0
0 0 0
,
00 0
0 0 0 00
GG
EA
k EA k
k EA k
EA
0
MN
0
(2.49d)
Here,
22
0 is a 2 by 2 zero matrix.
As seen in section 2.1, the state equations of a multi-body distributed system have the
same forms as Eq. (2.4) of a single-body distributed system. The s-domain solution of the
distributed system, Eq. (2.45), can be obtained through use of state transition matrix. The
state equations and the s-domain solution lay out a foundation for determination of the
eigensolutions and applications of the multi-body systems, as shall be seen in the
subsequent sections.
The vibration analysis of stepped systems that are assemblages of one-dimensional
distributed components (like elastic bars, torsional shafts and taut strings) and lumped
32
masses can be easily addressed based on the DTFM (Yang, 2010). For stepped second-
order distributed systems, the distributed transfer function formulation can be expressed
in simpler forms because matching matrices are not singular. Such stepped systems may
be considered as a special case of multi-body distributed parameter systems and the
details are explained in chapter 6.
2.4.2 Eigenvalue problems
Define the eigenvalue problem of the multi-body distributed system by vanishing the
excitation terms in Eq. (2.45).
()
( , ) ( ),
x
x s x x
x
F
(2.50)
subject to the conditions
( ) ( ) 0
G s G e
xx MN (2.51)
where and
se
xx are representing the end nodes of the subsystems, s is an eigenvalue, and
() x
the eigenfunction associated with s. Assume that
( ) ( , , )
s
x x x s a (2.52)
with a being a constant vector to be determined and ( , , ) xs the state transition matrix
which is the unique solution of (Chen, 1999)
( , , ) ( , ) ( , , ), , x s x s x s x
x
F (2.53)
The () x , which satisfies Eq. (2.53), is plugged in the condition (2.51), yields
( ) ( , , ) 0
G G G e s
s x x s Z a M N a = (2.54)
33
The characteristic equation of the system then is det ( ) 0
G
s Z , from which the system
eigenvalues , 1,2,
k
sk , are determined. The corresponding eigenfunction is obtained
by ( ) ( ,0, )
k k k
x x s a where
k
a is a nonzero solution of the homogeneous equation
( ) 0
Gk
s Za .
2.4.3 Applications
For the complex distributed system under harmonic excitations of the frequency ,
the frequency response or steady-state response of the system is determined by the s-
domain solutions with replacing s by j where 1 j . The solution is exact and of
closed form. To obtain the transient (time-domain) response of the multi-body distributed
system to initial and external disturbances, inverse Laplace transform of the s-domain
solution is required by the residue theorem.
2.5 Mid- and High-Frequency Analyses
Many engineering structures such as aircraft, satellite and rocket are often exposed to
mid- and high-frequency excitations. Prediction of mid- and high-frequency structural
response is important in such structures. A structural dynamic analysis, both at steady
state and at transient state, has been generally performed using the finite element methods
(FEM). It is well known that the FEM is not suitable for analysis in mid- and high-
frequency regions. Developing accurate and numerically stable methods for the analysis
of high-frequency vibration is an extremely challenging task owing to the numerical
34
instability in existing methods (Langley, 1991, 1995; Langley and Bardell, 1998; Wei, et
al., 2002). Wu (1997) investigated avoiding the numerical instability and matrix
singularity of high-frequency vibration analysis of a uniform beam using the DTFM. This
section is concerned with mid- and high-frequency analyses of distributed parameter
systems.
2.5.1 Eigensolutions
The determination of system eigensolutions is important in the proposed method and
given as Eq. (6.13). The evaluation of system eigensolutions often suffers numerical
instability at high frequencies. For instance, the numerical issues at high frequency ranges
can be observed in the model of an Euler-Bernoulli beam. For conventional analytical
methods, the system eigenvalues subject to certain boundary conditions are computed by
the roots of the characteristic equation which contains the hyperbolic sine and cosine
functions. These functions are exponential growing without bound. Thus, the general
forms offered in many vibration text books permit the evaluation of only the first 12
modes (Gonçalves, et al., 2007). For numerical methods such as the FEM, the required
nodal points are increasing as the frequency goes high. It is more obvious in a complex
system than a single-body system. Consequently, the FEM requires very extensive
computation cost and may be impractical to obtain acceptable results. The evaluation of
eigensolutions of single-body structures at high frequencies is investigated by many
researchers (Mead, 1994; Tang, 2003; Moreles, et al., 2005; Gonçalves, et al., 2007). One
of the advantages of the proposed method is that the numerical instability is simply
35
avoided to evaluate eigenvalues at high frequency regimes. It is applicable to both single-
body and multi-body distributed parameter systems. The problem of evaluating
eigensolutions of complex distributed parameter systems is represented in this section.
The high-frequency analysis method is illustrated on eigensolutions of a stepped beam.
Example 2.3: Eigensolutions of a stepped beam
A stepped beam composed of uniform Euler-Bernoulli beam elements is fixed at the
left end, is free at the right end, is constrained by a spring (k), and carries a lumped mass
(m); see Fig. 6.1. For numerical simulation, choose that 10 and 1. km The system
parameters are assigned as follows:
1 1 1
2 2 2
3 3 3
Segment1: 1.2, 10, 1.3
Segment 2 : 2.0, 16, 0.7
Segment 3: 1.0, 10, 1.0
EI l
EI l
EI l
(2.55)
Table 2.3 gives the 100
th
-104
th
and 200
th
-204
th
natural frequencies of the stepped
beam, which are obtained by the proposed method and finite element method. As the
number of elements increases, the FEM solutions converge to the solutions by the
proposed method.
2.5.2 Transient response
Conventional analytical methods such as eigenfunction expansion require exact system
eigenfunctions to obtain solutions of distributed parameter systems. An exact estimation
of system eigenfunctions is difficult and can be achieved only in limited cases involving
self-adjoint systems with uniform parameter distributions and simple boundary
36
conditions. When the boundary conditions are eigenvalue-dependent, or when the system
is non-self-adjoint, closed-form solutions are beyond reaching for most cases. Although
the finite element method is capable of modeling and analyzing the system in such cases,
it is computationally intensive, especially for systems with high-frequency dynamics. The
proposed method does not depend on system eigenfunctions, is able to accurately predict
solutions of a system.
By Eq. (2.19), the transient response of the system is given as
0 0 0
( , ) ( , , ) ( , ) ( , ) ( )
t L t
b
x t x t d d x t d
G p H (2.56)
Here, the system Green’s functions G and H are expressed by an infinite transfer function
residue series as Eq. (2.20). For an excitation whose spectrum is specified in a high-
frequency region or in combination of lower and higher frequency regions, the system
response can be estimated by a few modes which are dominant to the system. By the
selection of modes in the frequency ranges of interest, the Green’s functions can be
modified as
0
2
1
1
( , , ) 2 Re ( ,0, ) (0, , )
2 Re ( ,0, ) (0, , )
k
k
N
st
k k b k
k
N
st
k k b k
kN
x t x s s e
x s s e
G R M
RM
(2.57a)
0 2
1
1
( , ) 2 Re ( ,0, ) 2 Re ( ,0, )
kk
N N
s t s t
k k k k
k k N
x t x s e x s e
H R R (2.57b)
where Re{} is the real part of {} ;
k
s is the kth eigenvalue of the system.
37
The modification of the Green’s functions in the proposed method not only provides
efficiency for asymptotic transient analysis, but also saves computational cost in high-
frequency analysis.
Example 2.4: Transient analysis of a non-uniform beam
The proposed method is illustrated on high frequency analysis of a non-uniform
damped beam with end mass, damper and spring; see Figure 2.3. The lumped mass is
subject to a transverse force ()
b
ft . The transverse displacement ( , ) w x t of the beam is
governed by the differential equation
2 2 2
2 2 2
( ) ( , ) ( ) ( , ) ( ) ( , ) 0, (0, ) x w x t d x w x t EI x w x t x L
t t x x
(2.58)
with the boundary conditions
0
2 2 2
2 2 2
( , ) 0, ( , ) 0
( , ) 0, ( , ) ( ) ( , ) ( )
x
b
x L x L
xL
w x t w x t
x
w x t m c k w x t EI x w x t f t
x t t x x
Assume that zero initial disturbances. Here, () x , () EI x and () dx are the linear
density, bending stiffness and viscous damping of the beam which are functions of x,
respectively, m the inertia of the lumped mass, c the coefficient of the damper and k the
coefficient of the spring. In the numerical simulation, the system parameters () x , () EI x
and () dx have exponential distributions and is proportionally to one another and the
values of the system parameter are as follows.
38
2 2 2
0 0 0
0 0 0
( ) , ( ) , ( )
16, 1.6, 1, 0.2, 1, 4, 4, 8
x x x
x e d x d e EI x EI e
d EI L m c k
The boundary disturbance on the lumped mass is given by
12
( ) sin sin
b L H
f t f t f t
where
12
0.01, 1, 45.2
L
ff rad/s and 1700.9
H
rad/s. Note that the excitation
frequencies
L
and
H
are close to the 5th and 27th damped frequency of the beam in
Table 2.4, respectively. Fig. 2.4 shows the frequency response of the beam which is
represented in closed-form as seen in section 2.2. It is generally difficult to obtain exact
solutions in high frequency ranges by conventional analytical method because numerical
instability issues. However, the proposed method does not suffer such instability as
shown in Figure 2.4. In the calculation of the beam transient response, the infinite series
in Green’s functions can be truncated by Eq. (2.57). The truncated Green’s functions
include the first 8 modes as fundamental and low frequency excitation modes, and
26~28th modes as high frequency excitation modes. Figure 2.5 shows the transient
response by the proposed method and it is compared the results by the FEM. It is seen
that the proposed method is efficient in mid- and high-frequency ranges and provides
accurate response of the system subject to the high-frequency excitation.
39
k
m
l
1
l
2
l
3
w(x,t)
x
Figure 2.2 A stepped beam
Table 2.3 The natural frequencies
k
of the stepped beam
k The proposed method
Finite Element Method
60 elements 150 elements 300 elements
100
101
102
103
104
200
201
202
203
204
32041.310
32730.176
33126.681
33983.872
34745.689
129183.568
129844.526
131471.916
133112.997
133819.389
42792.108
43670.715
44936.410
46190.080
47173.771
32427.783
33095.137
33541.557
34440.894
35185.404
149391.617
151686.986
153708.337
155069.364
157372.127
32067.862
32755.419
33155.181
34015.513
34776.365
130640.953
131364.594
133139.844
134727.923
135448.153
40
x
w(x,t)
f
b
(t)
m
c k
ρ(x), EI(x), d(x)
Figure 2.3 A non-uniform beam with tip mass, damper and spring
Table 2.4 The first ten eigenvalues of the beam (i1 )
k s
k
1
2
3
4
5
6
7
8
9
10
26
27
28
− 0.3090 + 1.2180i
− 0.1313 + 4.4926i
− 0.0833 + 13.1848i
− 0.0679 + 26.7789i
− 0.0612 + 45.2987i
− 0.0576 + 68.7481i
− 0.0555 + 97.1296i
− 0.0542 + 130.4441i
− 0.0533 + 168.6923i
− 0.0527 + 211.8746i
− 0.0504 + 1573.8944i
− 0.0504 + 1700.9659i
− 0.0503 + 1832.9721i
41
Figure 2.4 The frequency response of the beam at x = 0.5
42
(a) (b)
(c) (d)
Figure 2.5 The transient displacement of the beam at x = 0.5;
(b)-(d): expanded views of the circle in (a).
43
Chapter 3. Vibration Analysis and Control of
Bars, Shafts and Strings
3.1 Introduction
Dynamic response prediction is important in modeling and analysis of flexible
vibrating systems, such as buildings, bridges, turbines, machines and aerospace structures
(Rayleigh, 1945; Timoshenko, 1981; Morse, 1948; Caughey and O’Kelly, 1965;
Meirovitch, 1967). Due to operational requirements and performance specifications,
flexible vibrating systems in many applications have spatially non-uniform distributions
of geometric and material properties, and are combined with lumped masses and energy
dissipative devices such as viscous dampers. These design considerations mathematically
lead to non-uniformly distributed systems that are non-proportionally damped (Caughey
and O’Kelly, 1965; Meirovitch, 1967). In vibration analysis of this type of a distributed
system, closed-form solutions in general are difficult to obtain; approximate methods are
commonly used. This chapter is concerned with the transient and frequency response
which is exact and of closed-form for a class of non-uniformly distributed systems,
namely, bars in longitudinal vibration, shafts in torsional vibration, and taut strings in
transverse vibration.
44
3.2 System Descriptions
Consider a non-uniformly distributed system with a length L subject to an external
load ( , ) f x t ; see Fig. 3.1. ( ), ( ) and ( ) x d x x , which are functions of spatial coordinate x,
are the inertia, stiffness and viscous damping parameters of the system, respectively. The
displacement ( , ) w x t of the system in consideration is governed by the differential
equation
2
2
( , ) ( , ) ( , )
( ) ( ) ( ) ( , ), (0, )
w x t w x t w x t
x d x x f x t x L
t t x x
(3.1a)
along with the boundary conditions
2
10 2
2
10 2
at 0 : ( )
at : ( )
LL
RR
ww
x m a a w t
tx
ww
x L m b b w t
tx
(3.1b)
where
L
m and
R
m are end masses;
0
a ,
1
a ,
0
b and
1
b are constants that are properly
assigned to characterize different types of boundaries; and ()
L
t and ()
R
t are related to
prescribed boundary excitations (external load or possible foundation motion). Eq. (3.1a)
is a model for elastic bars in longitudinal vibration, circular shafts in torsional vibration,
and taut strings in transverse vibration. The system is also subject to the initial conditions
00
( ,0)
( ,0) ( ), ( ), 0,
wx
w x u x v x x L
t
(3.1c)
where
0
() ux and
0
() vx are given profiles of initial displacement and velocity of the
system.
45
In vibration analysis, if a vibrating continuum is proportionally damped, its transient
response can be expressed by an infinite series of the eigenfunctions that are associated
with the corresponding undamped self-adjoint system (Caughey and O’Kelly, 1965). The
non-uniform system described by Eq. (3.1) is proportionally damped if it does not have
end masses and dampers and if the damping parameter is proportional to the inertia
parameter:
0 0 0
0, 0
( ) ( ), ( ) ( ), ( ) ( )
LR
mm
x A x d x d A x x B x
(3.2)
where
0 0 0
,, d are constant, and () Ax and () Bx are known positive functions. The
non-uniform system in general is non-proportionally damped when the system has the
end masses and/or when the damping parameter () dx is not proportional to the inertia
parameter () x . It is well known that the eigenfunctions of a non-proportionally damped
vibrating system are complex and non-orthogonal. The boundary controller which
contains dissipation devices such as a damper causes the system to be non-proportionally
damped.
Figure 3.1 Schematic of a non-uniformly distributed system
46
3.3 Distributed Transfer Function Formulation
In this work, the boundary-initial value problem described in the previous section is
solved by the development of a new distributed transfer function formulation, which is an
extension distributed transfer function analysis of uniform systems (Yang and Tan, 1992;
Yang, 2005). To this end, take Laplace transform of Eq. (3.1) with respect to time, and
cast the resulting equations into a spatial state form
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ), (0, ) x s x s x s x s x L
x
Fp (3.3a)
ˆ ˆ (0, ) ( , ) ( )
b b b
s L s s MN (3.3b)
where the over-hat stands for Laplace transform with respect to time t, s is the Laplace
transform parameter, and
( , )
() 0
( , )
ˆ ˆ ( , ) , ( , ) , ( )
( , )
() 1 ()
L
eI
b
R
w x s
s
p x s
x s x s s
w x s
s x
x
p
2
01
( , )
( ) ( ) 1 ( )
( ) ( )
xs
x s d x s d x
x x dx
F (3.4)
2
01
2
01
00
,
00
L
bb
R
m s a a
m s b b
MN
with
0 0 0
0 0 0
0 0 0
ˆ
( , ) ( , ) ( ) ( ) ( ) ( ) ( )
ˆ ( ) ( ) (0) (0) (0)
ˆ ( ) ( ) ( ) ( ) ( )
eI
L L L L
R R R L
p x s f x s x su x v x d x u x
s s m su v d u
s s m su L v L d u L
(3.5)
47
Thus, the original boundary-initial value problem is converted to the s-domain state
equation (3.3a) subject to the boundary condition (3.3b).
As a special case, for a proportionally damped system satisfying the conditions (3.2),
the state and boundary matrices in Eq. (3.4) become
2
0
01
01
01
( , ) ( ) 1 ( )
( ) ( )
00
,
00
bb
xs A x dB x
B x B x dx
aa
bb
F
MN
(3.6)
where
22
0 0 0
( ) / s d s , and the boundary matrices ,
bb
MN are not functions of s.
Note that the matrices in Eq. (3.6) are the same in format as those for the corresponding
undamped system if
2
in the equation is has the form
22
00
/ s . This resemblance
in state representation will be utilized later on to determine system eigenvalues; see
Section 3.6.
3.4 Dynamic Analysis
3.4.1 The s-domain Solution
The s-domain solution of Eq. (3.3) can be expressed in terms of the distributed transfer
functions
0
ˆ ˆ
ˆ ˆ ( , ) ( , , ) ( , ) ( , ) ( ), , (0, )
L
b
x s x s x s d x s s x L
G p H (3.7)
where
ˆ
G and
ˆ
H are the distributed transfer functions of the system,
48
1
ˆ
( , ) (0, , ),
ˆ
( , , )
ˆ
( , ) ( , , ),
ˆ
( , ) ( ,0, ) ( )
b
b
x s s x
xs
x s L s x
x s x s s
HM
G
HN
HZ
(3.8)
with
( ) ( ,0, )
bb
s L s Z M N (3.9)
Equations (3.7) and (3.8) are similar to the transfer function formulation although
uniform and stepped systems are considered therein. The ( , , ) xs in Eqs. (3.8) and (3.9)
is the state transition matrix of the system.
3.4.2 Eigenvalue Problem
In the state form given in the section 3.3, the eigenvalue problem is described by
()
( , ) ( ), (0, )
x
x s x x L
x
F
(3.10a)
subject to the boundary condition with control forces
(0) ( ) 0
bb
L MN
(3.10b)
where s is an eigenvalue, and () x the eigenfunction associated with s. Equation (3.10) is
obtained by vanishing disturbances in Eq. (3.3). The solution of Eq. (3.10a) has the form
( ) ( , , ) x x s a
(3.11)
where a is a constant vector to be determined. Equation (3.11) is inserted into Eq.
(3.10b) to yield the transcendental characteristic equation of the system
det ( ) det ( , , ) 0
bb
s L s Z M N
(3.12)
49
where Eq. (3.9) has been used. The roots of Eq. (3.12)
are the eigenvalues of the system.
This implies that the system eigenvalues are also the poles of the distributed transfer
functions given in Eq. (3.8).
3.4.3 Frequency Response
Assume that the non-uniform system, which can be either proportionally or non-
proportionally damped, is subject to harmonic excitations.
( , ) ( )e , ( ) e , ( ) e
j t j t j t
L L R R
f x t F x t t
(3.13)
where 1 j , and is the excitation frequency. The steady-state response (the
frequency response) of the system, by Eq. (3.7), is given by
( , ) ( , e , (0, )
jt
ss
x t x j x L
(3.14)
where
0
0
()
ˆ ˆ
( , ( , , ) d ( , ) , (0, )
1 ()
L
L
R
F
x j x j x j x L
GH (3.15)
and
ˆ
( , , ) xj G and
ˆ
( , ) xj H are the frequency response functions of the system. With
Eq. (3.8), the frequency response functions can be evaluated in the exact and closed form.
3.4.4 Transient Response
The transient response of the damped non-uniform system is obtained through inverse
Laplace transform of Eq. (3.7), which produces the following Green’s function formula
50
0 0 0
0 0 0
0
00
0
() 0
( , )
( , ) ( , , ) ( , )
() 1 ()
0
( ) ( )
( , , ) ( ) ( , , ) ( ) ( , , ) ( )
1 ( ) ( )
(0) (0)
( , ) ( , )
()
t L t
L
R
L
LL
R
f
x t x t d d x t d
d
x t u x t v x t u d
t
m u m v
x t x t
m u L t
ηGH
G G G
HH
0
00
(0)
( ) ( )
L
RR
du
m v L d u L
(3.16)
where the matrix Green’s functions G and H are the inverse Laplace transforms of the
transfer functions
ˆ
G and
ˆ
H , respectively; and ( , ) xt is the inverse Laplace transform of
ˆ ( , ) xs , which by Eq. (3.4) contains both the displacement and strain parameters of the
system. The first two terms on the right-hand side of Eq. (3.16) represent the system
response caused by external and boundary excitations and the remaining terms describe
the contributions of initial disturbances.
Following (Yang, 2008, 2010), the Green’s functions of the system are determined as
1
1
( , , ) 2 Re ( ,0, ) (0, , )e
( , ) 2 Re ( ,0, ) e
k
k
st
k k b k
k
st
kk
k
x t x s s
x t x s
G R M
HR
(3.17)
where
1
Res ( ) , 1,2,...
k
k
ss
sk
RZ (3.18)
which is called as the transfer function residue associated with the kth eigenvalue
k
s .
According to Eqs. (3.9) and (2.15), an explicit expression of () s Z can be obtained, as
shown in Chapter 2.
51
3.5 Specific Distributed Systems
The key in the proposed method is to determine a state transition matrix for the non-
uniformly distributed systems. A state transition matrix of the system can be determined
by Eq. (2.7). Particular non-uniform systems have exact fundamental matrices which are
presented in this section. The state transition matrix of a non-uniformly distributed
system, in general, cannot be expressed in exact and closed-form and has to be evaluated
numerically when no existence of exact solutions as shown in section 2.2.
A fundamental matrix satisfying that ( , ) / ( , ) ( , ) x s x x s x s U F U , by the theory of
differential equations, is
12
12
( , ) ( , )
( , )
( , ) ( , )
u x s u x s
xs
u x s u x s
xx
U (3.19)
where
1
( , ) u x s and
2
( , ) u x s are two independent solutions of the s-domain differential
equation
2
( , )
( ) ( ) ( , ) ( ) 0, (0, )
u x s
x s d x s u x s x x L
xx
(3.20)
With
1
( , ) u x s and
2
( , ) u x s , the system transition matrix can be constructed by Eqs.
(2.7) and (3.19). Independent solutions of Eq. (3.20) for several types of damped systems
have been derived, and they are presented in sequel. For the derivation of these solutions,
see the Appendix A.
52
A. Homogeneous systems
For a homogeneous system, its inertia, stiffness, and damping parameters of the system
are proportional to one another, and defined by
0 0 0
( ) ( ), ( ) ( ), ( ) ( ) x A x d x d A x x A x (3.21)
where
0
,
0
d , and
0
are constants, and () Ax is a positive function of x. With (3.21), Eq.
(3.20) becomes
2 2
2 00
2
0
() ( , ) 1 ( ) ( , )
( , ) 0,
()
s d s u x s dA x u x s
u x s
x A x dx x
(3.22)
Independent solutions of (3.22) are obtained for the following seven types of A(x).
Type 1. Uniform: () A x A (costant)
12
,
xx
u e u e
(3.23)
Type 2. Linear: ( ) 1 A x ax
1 0 2 0
( 1/ ) , ( 1/ ) u I x a u K x a (3.24)
Type 3. Exponential:
2
( ) e
ax
Ax
For
22
0 a ,
2 2 2 2
12
e cosh( ), e sinh( )
ax ax
u a x u a x (3.25a)
and for
22
0 a ,
12
e , e
ax ax
u u x (3.25b)
Type 4. Conical:
2
( ) (1 ) A x ax
11
12
1 cosh( ), 1 sinh( ) u ax x u ax x
(3.26)
53
Type 5. Cosine square:
2
( ) cos ( ) A x ax b
For
22
0 a ,
2 2 2 2
12
cosh( ) sinh( )
,
cos( ) cos( )
a x a x
uu
ax b ax b
(3.27a)
and for
22
0 a ,
11
12
cos ( ), cos ( ) u ax b u x ax b
(3.27b)
Type 6. Catenoidal:
2
( ) cosh A x ax b
For
22
0 a ,
2 2 2 2
12
cosh( ) sinh( )
,
cosh cosh
a x a x
uu
ax b ax b
(3.28a)
and for
22
0 a ,
11
12
cosh , cosh u ax b u x ax b
(3.28b)
Type 7. Power function: ( ) (1 ) , 0
b
A x ax b
Write (1 ) / 2 nb . For n to be an integer
12
(1 ) ( 1/ ) , (1 ) ( 1/ )
nn
nn
u ax I x a u ax K x a (3.29a)
and for n to be a non-integer
12
(1 ) ( 1/ ) , (1 ) ( 1/ )
nn
nn
u ax I x a u ax I x a
(3.29b)
where ,
nn
IK are modified Bessel functions of order n.
B. Inhomogeneous systems
In this study, we examine a class of inhomogeneous systems, whose inertia and
stiffness parameters of the system are not proportional to each other, but whose damping
and inertia parameters are proportional to each other. This means that
0 0 0
( ) ( ), ( ) ( ), ( ) ( ) x A x d x d A x x A x
(3.30)
54
Substituting with Eq. (3.30) into Eq. (3.20) gives
2
2 00
0
( , )
( ) ( ) ( , ) 0,
s d s u x s
A x A x u x s
xx
(3.31)
Two types of parameter distributions are considered as follows.
Type 8. Power distribution: ( ) (1 ) , ( ) (1 )
bc
A x x A x x
, bc
If 2 bc , with the transformation 1 x , Eq. (3.31) becomes an Euler equation.
12
2 2 2
12
(1 ) , (1 ) , for ( 1) 4 / 0 u x u x c
(3.32a)
12
2 2 2
12
(1 ) , (1 ) ln(1 ), for ( 1) 4 / 0 u x u x x c
(3.32b)
where
1
and
2
are the roots of the characteristic equation
2 2 2
( 1) / 0 c .
If 2 bc , with the coordinate transformation
1
2
2
2
( , ) (1 ) ( , )
2
(1 )
( 2)
b
cb
u x s x z s
x
cb
Eq. (3.31) becomes a modified Bessel equation with order
1
2
b
n
cb
.
12
22
1
12
22
2
2
(1 ) (1 )
( 2)
2
(1 ) (1 )
( 2)
b c b
n
b c b
n
u x I x
cb
u x K x
cb
(3.33a)
for n to be an integer,
55
12
22
1
12
22
2
2
(1 ) (1 )
( 2)
2
(1 ) (1 )
( 2)
b c b
n
b c b
n
u x I x
cb
u x I x
cb
(3.33b)
and for n to be an integer.
Type 9. Exponential distribution : ( ) e
ax
Ax
, ( ) ,
bx
A x e a b
Forab , refer to Type 3 in homogeneous systems. Forab , with the coordinate
transformation
2
( , ) ( , ), with / ( )
2
n
ab
x
u x s z s n a a b
e
ab
Eq. (3.31) becomes a modified Bessel equation with order n. The independent solutions
of (3.31) are given by
2 2 2 2
12
2 2 2 2
,
nn
a b a b a b a b
x x x x
nn
u e I e u e K e
a b a b a b a b
(3.34a)
for n to be an integer,
2 2 2 2
12
2 2 2 2
,
nn
a b a b a b a b
x x x x
nn
u e I e u e I e
a b a b a b a b
(3.34b)
and for n to be a non-integer. Here ,
nn
IK are modified Bessel functions of order n.
The solutions presented in this section cover a variety of distributed systems in
engineering problems.
56
3.6 Boundary Control
3.6.1 Boundary Controllers
Consider the system with boundary controllers; see Figure 3.2. The boundary control
forces in Laplace-domain can be expressed by the form
10
ˆ
( ) ( ) (0, ) f s c s w s and
2
ˆ
( ) ( ) ( , )
L
f s c s w L s (3.35)
where
0
() cs and ()
L
cs are transfer functions describing the controller.
By the distributed transfer function formulation, the state equations are expressed by
Eq. (3.3). With the boundary controllers, the boundary matrices in Eq. (3.3b) are
modified as
2
0 0 1
2
01
00 ()
,
() 00
L
bb
RL
m s a c s a
m s b c s b
MN (3.36)
3.6.2 Stability of Root-Locus Analysis
It is well known that the stability is defined by determining the eigenvalues of the
system. Two root finding algorithms for solution of Eq. (3.12) are developed based on
system damping status, that is, proportional or non-proportional damping. Without loss of
generality, assume that the distributed system is under-damped. The system eigenvalues
are then in the complex conjugate form , 1,2,3,...
k k k
s j k
, 1 j . The root
finding algorithms are described as follows.
57
A. Root finding for systems with proportional damping
For a proportionally damped system satisfying conditions (3.21), its state
representation is the same in format as that for the corresponding undamped system, as
discussed in the previous section. By Eq. (3.12), the characteristic equation for both the
systems can
( ) det ( ) 0 s Z (3.37)
where
22
0 0 0
( ) / s d s for the damped system, and
22
00
/ s for the undamped
system. Denote the roots of Eq. (3.37) for the undamped system by
k
jb , 1,2,3,... k ,
1 j , where
k
b are positive numbers. The root finding algorithm for a proportionally
damped system thus takes two steps:
(i) Determine the roots
k
b of the undamped system, that is, ( ) 0
k
jb .
(ii) By letting
22
0 0 0
( ) ( ) /
k
jb s d s , the eigenvalues of the proportionally damped
system are obtained as follows
2
2
0 0 0
0 0 0
11
, 1,2,3,..., 1
24
k
k
d b d
s j k j
(3.38)
With this algorithm, only real roots
k
b need to be determined. Accurate and efficient real
root finding methods have been well developed, and can be adopted; see Reference
(Meirovitch, 1980) for instance.
58
B. Root finding for systems with non-proportional damping
For a non-proportionally damped system, solution of the characteristic equation (3.12)
requires root finding in the two-dimensional complex plane, which is different from that
for a proportionally damped system. In this study, a root locus technique is proposed,
which takes the following three steps.
Step 1. Modify the characteristic equation (3.12) by multiplying all the damping
coefficients in Eq. (3.33) by a non-dimensional parameter
00
(0 1) : ( ) ( ) c x c x ,
( ) ( )
LL
c x c x . This yields a new characteristic equation
( , ) 0 s (3.39)
which is a function of s and . When 0 , Eq. (3.39) becomes the characteristic
equation of the corresponding undamped system; when 1 , Eq. (3.39) is the same as
Eq. (3.12).
Step 2. Set 0 and find the roots (natural frequencies)
k
of the corresponding
undamped system ( ,0) 0, 1,2,...
k
jk . The solution process is similar to that for Eq.
(3.37).
Step 3. Determine the complex eigenvalues of the damped system via an iterative
procedure. Let the iteration step be 1/ N , where N is the number of iterations. To
determine the kth eigenvalue (
k k k
sj ), undertake the following two tasks.
(i) Set an initial guess as
(0)
kk
zj and find a root
(1)
k
z in the neighborhood
(0)
k
z in the
complex plane by solving ( , ) 0 z .
59
(ii) For 2,3,..., mN , set an initial guess
( 1) m
k
z
and determine a root
() m
k
z in the
neighborhood
( 1) m
k
z
in the complex plane from ( , ) 0 zm .
In solution of ( , ) 0 zm for 1 mN , standard solvers for simultaneous
nonlinear algebraic equations can be adopted. Hence, for 0, , 2 ,..., , N we
obtain a sequence of complex numbers
(0) (1) ( ) ( )
, , ,...,
mN
k k k k
z z z z , which forms a root locus.
Because 1 N , the eigenvalues of the damped system in the kth mode are given by
()
,
N
k k k k
s z s s
(3.40)
where
k
s is the complex conjugate of
k
s .
Once an eigenvalue s is known, the associate eigenfunction can be determined through
determination of a non-trivial vector a in Eq. (3.11). However, eigenfunctions are not
needed in the proposed solution method.
Figure 3.2 Schematic of a controlled non-uniform system
60
3.7 Numerical examples
The proposed transient analysis is illustrated on the following three examples: (1) A
non-uniform shaft subject to an external toque. (2) A non-uniform shaft subject to
boundary controllers.
Example 3.1: A non-uniform shaft in torsional vibration
In Fig. 3.3, a non-uniform shaft with an end disk is in torsional vibration. The rotation
(twist) of the shaft under a torque () t at x
is governed by
2
2
( , ) ( , )
( ) ( ) ( ) ( ), (0, )
x t x t
J x GJ x t x x x L
t x x
(3.41)
with the boundary conditions
0
2
2
(0, )
(0, ) (0) 0
( , ) ( , )
( , ) ( ) 0
DL
t
k t GJ
x
L t L t
I k L t GJ L
tx
(3.42)
Assume conical distributions for the system parameters
22
0 0 0 0
( ) (1 / ) , ( ) (1 / ) J x J x L GJ x G J x L (3.43)
The state matrix F
and the boundary matrices in Eq. (3.4) are
0
2
2
01
00 0
( , ) , ,
( ) ( )
() 00
( ) ( )
bb
DL
k
xs
J x s GJ x
I s k GJ L
GJ x GJ x
F M N (3.44)
with ( ) ( ) / GJ x dGJ x dx . By Eq. (3.19) and (3.26), a fundamental matrix of the bar can
be written as
61
22
cosh( ) 1 sinh( )
1 / 1 /
( , )
sinh( ) cosh( ) cosh( ) sinh( )
1 / (1 / ) 1 / (1 / )
xx
x L x L
xs
x x x x
x L L x L x L L x L
U (3.45)
where
00
/
s
.
For numerical simulation, choose the non-dimensional values of the system parameters
as
0 0 0 0 0
2, 240, 0.2, 1, 25, 10, 20
DL
J G J L I k k . The first 15 natural
frequencies of the shaft are listed in Table 3.1, where the proposed method is compared
with the FEM. Let the shaft be initially at rest. Consider an external torque
0
( ) T (1 )
t
te
applied to the shaft at point x
. Assume that
0
T 2, 0.1 and
/ 0.5 xL
. The forced vibration of the shaft is computed by Eq. (3.16), (3.17), and
using the first 40 modes, and plotted in Figs 3.4 and 3.5. Figures 3.4(a) and 3.4(b) show
the spatial distributions of the rotation and shear strain of the shaft at t = 0.1, 0.3, 0.5, 0.7
and 1.0. Figure 3.5 gives a time history (0 20) t of the shaft vibration at the end disk
location() xL which is compared with the result by the FEM. The proposed method is
efficient to analyze response of the system under high frequency excitations. When the
excitation
0.1
( ) sin3355.4
t
t e t
is applied at point x
, the transient response is shown in
Fig. 3.6 and compared with those by the FEM with 30, 150 and 300 elements. Here, the
excitation frequency, 3355.4 rad/s, is close to the 99th mode. The transient response is
computed by Eq. (3.16) with the modified Green's function.
62
40
1
100
98
( , , ) 2 Re ( ,0, ) (0, , )e
2 Re ( ,0, ) (0, , )e
k
k
st
k k b k
k
st
k k b k
k
x t x s s
x s s
G R M
RM
(3.46)
The truncation of modes is required to compute the Green's function. So, the proper
selection of modes in Eq. (3.46) reduces the computation rate as well as gives good
approximations of the Green's functions which are important to achieve transient
response.
63
Figure 3.3 A shaft with end disk in torsional vibration
Table 3.1 The natural frequencies
k
of the non-uniform shaft (in rad/s)
k
DTFM
Finite Element Method
30 elements 150 elements 300 elements
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1.05435
16.43223
51.37162
85.88642
120.34367
154.78185
189.21139
223.63627
258.05836
292.47864
326.89769
361.31586
395.73338
430.15040
464.56704
1.05435
16.43516
51.42755
86.13710
121.02574
156.22739
191.84847
227.98954
264.74915
302.22453
340.51092
379.70003
419.87766
461.12007
503.48906
1.05435
16.43234
51.37386
85.89644
120.37091
154.83952
189.31647
223.80953
258.32436
292.86572
327.43796
362.04523
396.69156
431.38089
466.11714
1.05435
16.43226
51.37218
85.88893
120.35048
154.79627
189.23765
223.67958
258.12484
292.57538
327.03271
361.49812
395.97280
430.45783
464.95428
64
(a)
(b)
Figure 3.4 Spatial profiles of the forced vibration of the shaft:
(a) Rotation distributions; (b) shear strain distributions
65
Figure 3.5 The forced vibration of the shaft at x = L;
solid line −: DTFM; dotted line −: FEM with 30 elements
66
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1
-0.5
0
0.5
1
1.5
2
2.5
3
3.5
x 10
-5
x
Rotation, (x,t=0.1)
30-Element FEM Model
Propoesed Method
(a)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-3
-2
-1
0
1
2
3
4
x 10
-5
x
Rotation, (x,t=0.1)
150-Element FEM Model
Propoesed Method
(b)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1
-0.5
0
0.5
1
1.5
2
2.5
3
3.5
x 10
-5
x
Rotation, (x,t=0.1)
Propoesed Method
300-Element FEM Model
(c)
Figure 3.6 The transient responses ( ,0.1) x under the external excitation
0.1
( ) sin3355.4 ( )
t
t e t x x
with / 0.5 xL
.
67
Examples 3.2: A non-uniform bar subject to boundary controllers
The displacement ( , ) u x t of a damped bar is governed by the partial differential
equation
2
2
( , ) ( , ) ( , )
( ) ( ) ( ) ( , ), (0, )
u x t u x t u x t
x d x EA x f x t x L
t t x x
(3.46)
where ( ), ( ), ( ) x d x EA x are the linear density, the viscous damping coefficient and the
longitudinal rigidity of the bar, respectively. The distributions of the system parameters
are of exponential form:
2 / 2 / 2 /
0 0 0
( ) , ( ) , ( )
x L x L x L
x e d x d e EA x EA e
(3.47)
A fundamental matrix of the bar, by Eq. (3.19) and (3.25a), is obtained as
/
1
cosh( ) sinh( )
( , )
cosh( ) sinh( ) sinh( ) cosh( )
xL
mx mx
m
x s e
mx m mx mx mx
L mL
U (3.48)
where
2 2 2
0 0 0
( ) / / m s d s EA L .
Assume that the bar is fixed at the left end (x = 0) and a control force
2
() ft is applied
at the right end (x = L).
2
(0, ) 0, ( ) ( , ) ( ) u t EA L u L t f t
x
(3.49)
In numerical simulation, the following non-dimensional values for the system
parameters are assigned:
0 0 0
1, 10, 1, 2, 1. EA d L In Eq. (3.36), a controller
is designed with a spring and a linear damper, i.e. ( ) (1 )
L
c s s with 0.04 . The
eigenvalues of the controlled bar are listed in the third column of Table 3.2. For
68
comparison, the eigenvalues of the corresponding undamped system
0
( 0 and 0) d
are listed in the first column of Table 3.2. The eigenvalues of the uncontrolled bar
0
( 0 and 0) d are computed by Eq. (3.38) and are listed in the second column of the
table. The root locus method proposed in section 3.6 is utilized to determine the
eigenvalues of this non-proportionally damped and controlled system. Figure 3.6 shows
the first three eigenvalue loci of this controlled bar with 0 ≤ μ ≤ 2.
The frequency response of the bar at the right end (x = L)
12
0
ˆ
( ) 1 0 ( , )
1
h L j
H (3.50)
is computed, and plotted against the frequency parameter in Figure 3.7.
An external force
0
( ) (1 e )
t
f t F
is applied to the bar at point
f
x . Let the bar be at
rest initially. The transient response of the bar is
0
1
0
2
( , ) Re ( ,0, ) (0, , ) ( )
1 ()
k k b f k k
k f
F
x t x s x s q t
x
RM (3.51a)
where
()
0
ee 1
( ) e (1 e )d
()
k
k
t st t
st k
k
k k k
s
qt
s s s
(3.51b)
For
0
2, 0.1and / 0.5
f
F x L , the forced response of the bar is computed with
first 60 terms in Eq. (3.16) and (3.17), and plotted in Fig. 3.8. Figures 3.8(a) and 3.8(b)
show the spatial distributions of the displacement and strain ( / ux ) of the bar at times t
= 0.3, 0.5, and 1.0.
69
Table 3.2 The eigenvalues s
k
of the non-uniform bar (i1 )
k Undamped Uncontrolled Controlled
1
2
3
4
5
6
7
8
9
10
15
20
50
9.61208i
17.28509i
26.37193i
35.89281i
45.58665i
55.36469i
65.18963i
75.04323i
84.91560i
94.80092i
144.32879i
193.93078i
491.84345i
-0.5 + 9.59906i
-0.5 + 17.27785i
-0.5 + 26.36719i
-0.5 + 35.88932i
-0.5 + 45.58391i
-0.5 + 55.36244i
-0.5 + 65.18771i
-0.5 + 75.04156i
-0.5 + 84.91413i
-0.5 + 94.79960i
-0.5 + 144.32792i
-0.5 + 193.93013i
-0.5 + 491.84319i
-1.4260 + 9.97397i
-2.3189 + 17.85894i
-2.7163 + 26.90588i
-2.9021 + 36.34817i
-2.9990 + 45.97287i
-3.0546 + 55.69624i
-3.0889 + 65.47831i
-3.1116 + 75.29836i
-3.1272 + 85.14369i
-3.1383 + 95.00687i
-3.1646 + 144.46665i
-3.1738 + 194.03406i
-3.1834 + 491.88443i
Figure 3.7 The eigenvalue loci of the bar with a controller (0 ≤ μ ≤ 2)
70
(a)
(b)
Figure 3.8 The frequency response of the bar at the right end: amplitude (a) and
phase (b); dash line –: uncontrolled, solid line –: with a controller (μ=0.04)
71
(a)
(b)
Figure 3.9 Spatial profiles of the forced vibration of the bar:
(a) displacement distributions; (b) strain distributions; dash
line −: uncontrolled, solid line −: controlled (μ=0.04)
72
Chapter 4. Analysis and Control of Acoustic Fields in a Duct
4.1 Introduction
The analysis and control of acoustic ducts with temperature gradients is a problem of
considerable interest in many engineering applications. This chapter is concerned with
analytical solutions of transient response of one-dimensional duct systems with axial
mean temperature gradients that are under arbitrary sound generations, boundary
disturbances and initial conditions. In this investigation, several types of non-uniform
temperature distributions are considered. By adopting the concept of time delay relation
(Yang and Mote, 1992), active noise control of the duct systems with non-uniform
temperature distributions is studied from a theoretical point of view. Through the
introduction of a specific time delay block(s) in the feedback control system, a closed-
loop control design which is capable of stabilizing the system is provided. Through
examples, the feasibility and efficiency of delayed feedback control method are verified.
4.2 System Description
For a perfect, inviscid and non-heat-conducting gas that is in a constant area duct, its
one-dimensional governing equations are as follows (Sujith, 1995, 1996; Kumar and
Sujith, 1997, 1998; Cummings, 1977)
Momentum:
x
u u p
uf
t x x
(4.1)
73
Continuity:
u
uq
t x x
(4.2)
Constitutive:
2
pp
u c u
t x t x
(4.3)
State: p RT (4.4)
where ρ and u are the density and velocity of the gas, p and T the acoustic pressure and
temperature in the duct, c the speed of sound, f
x
an external force per unit volume, q mass
influx per unit volume, and R the specific gas constant.
Write the parameters in the previous equations as
( , ) ( ) ( , ), ( , ) ( ) ( , ), ( , ) ( ) ( , ) x t x x t u x t u x u x t p x t p x p x t (4.5)
where the over-bar stands for steady solutions, and the prime small-amplitude time-
dependent quantities. Assume that the effect of the mean flow is negligible, which means
that the proposed method is valid only for mean Mach numbers less than 0.1 (Sujith et al.,
1995). With Eq. (4.1) to (4.5), a linearized wave equation with a mean temperature
profile T is obtained as follows (Sujith et al., 1995; Cummings, 1977)
22
2 2 2
1
,0
p dT p p
S x L
x T dx x c t
(4.6)
where c RT with /
PV
CC being the ratio of specific heats,
1
( , )
x
x
f qd
S S x t f
t x dx
(4.7)
and L is the length of the duct. The boundary conditions of the duct can be written as
74
10
10
( , )
at 0 : ( , ) ( )
( , )
at : ( , ) ( )
L
R
p x t
x a a p x t t
x
p x t
x L b b p x t t
x
(4.8)
where
0 1 0 1
, , and a a b b are constants that are properly assigned to characterize different
types of boundaries, ()
L
t and ()
R
t are related to prescribed boundary excitations.
Additionally, the system is subject to the initial conditions
00
( ,0)
( ,0) ( ), ( ), 0
px
p x p x v x x L
t
(4.9)
with
0
() px and
0
() vx being determined by the profiles of initial acoustic pressure and its
temporal derivative.
4.3 Distributed Transfer Function Formulation
In this investigation, a transfer function formulation is first derived, based on which the
s-domain response of the duct is obtained. Transient solutions are then determined
through inverse Laplace transform of the s-domain response. To obtain a transfer
function formulation for the boundary-initial problem, take Laplace transform of Eq.
(4.6) and (4.8) with respect to time, and cast the results in the following spatial state form
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ), (0, ) x s x s x s x s x L
x
Fq (4.10a)
ˆ ˆ (0, ) ( , ) ( )
b b b
s L s s MN (4.10b)
where the over-hat stands for Laplace transform, s is the Laplace transform parameter,
ˆ ( , ) xs is the state vector, and
75
01
2
01
2
00 2
ˆ ( , )
ˆ () 0
ˆ ˆ ( , ) , ( , ) ( , ) , ( )
ˆ ( , )
ˆ () 1
01
00
( , ) , ,
1
00
1
ˆ
( , ) ( , ) ( ) ( )
L
eI b
R
bb
eI
p x s
s
x s x s S x s s
p x s
s
x
aa
xs
s dT
bb
c T dx
S x s S x s sp x v x
c
q
F M N
(4.11)
In derivation of Eq. (4.10), the initial conditions (4.9) have been used. Thus, the
original boundary-initial problem is converted to an equivalent s-domain state equation
subject to a boundary condition.
4.4 Dynamic Response
The s-domain response can be constructed through use of state transition matrix. With
the state transition matrix, the solution of Eq. (4.10) takes the form
0
ˆ ˆ
ˆ ˆ ( , ) ( , , ) ( , ) ( , ) ( ), , (0, )
L
b
x s x s s d x s s x L
G q H (4.12)
where
ˆ
G and
ˆ
H are the distributed transfer functions of the duct, and they are given
1
ˆ
( , ) (0, , ),
ˆ
( , , )
ˆ
( , ) ( , , ),
ˆ
( , ) ( ,0, ) ( ,0, )
b
b
bb
x s s x
xs
x s L s x
x s x s L s
HM
G
HN
H M N
(4.13)
Eigenvalue Problem: The eigenvalue problem is the state form is expressed by the
homogeneous equation
()
( , ) ( ), (0, )
x
x s x x L
x
F
(4.14a)
76
(0) ( ) 0
bb
L MN (4.14b)
where s is an eigenvalue, and () x the eigenfunction associated with s. It is easy to see
that the previous equations are equivalent to the standard eigenvalue problem for the
acoustic system. The solution of Eq. (4.14a) has the form (Yang, 2005)
( ) ( , , ) x x s a
(4.15)
where a is a constant vector to be determined. Substituting of Eq. (19) into Eq. (18b)
yields the transcendental characteristic equation of the system
det ( , , ) 0
bb
Ls MN (4.16)
The roots of Eq. (4.17)
are the eigenvalues of the system. For a duct without damping, its
eigenvalues are given by
kk
sj for k = 1,2,…, where 1 j and
k
is the kth
eigenfrequency.
The eigensolutions can be used to compute the acoustic pressure and velocity of
harmonic time dependency, for the duct under no external excitations (S = 0 in Eq. (4.6)).
To this end, write the acoustic pressure and velocity as
( , ) ( ) , ( , ) ( )
j t j t
p x t P x e u x t U x e
(4.17)
Assume that the acoustic pressure at a point
*
x of the duct is nonzero, i.e.,
*
( ) 0 Px ,
*
0 xL . By Eqs. (4.16) and (4.17) and linearized momentum equation 0
up
tx
,
it can be shown that
1 * 2
*
1 * 1 *
( , , ) ( ) ( , , )
( ) ( ) , ( )
( , , ) ( , , )
x j P x x j
P x P x U x
x j j x j
aa
aa
ee
ee
(4.18)
77
where
1
10 e ,
2
01 e , and is an eigenfrequency of the duct.
Transient Solutions: For a duct with a mean temperature distribution, its transient
response is obtained through inverse Laplace transform of Eq. (4.12), which produces the
following Green’s function formula
0 0 0
00 2
0
() 0
( , ) ( , , ) ( , ) ( , )
() 1
0
1
( , , ) ( ) ( , , ) ( )
1
t L t
L
R
L
x t x t S d d x t d
x t p x t v d
ct
GH
GG
(4.19)
where the matrix Green’s functions G and H are the inverse Laplace transforms of the
transfer functions
ˆ
G and
ˆ
H , respectively; ( , ) xt is the inverse Laplace transform of
ˆ ( , ) xs , and
2
() c RT x .
The first two terms on the right-hand side of Eq. (4.20) represent the system response
caused by external and boundary generations, and the last term describes the
contributions of initial disturbances. The Green’s functions are determined as
1
1
( , , ) 2 Re ( ,0, ) (0, , )
( , ) 2 Re ( ,0, )
k
k
st
k k b k
k
st
kk
k
x t x s s e
x t x s e
G R M
HR
(4.20)
where
kk
sj , with 1 j and
k
being the kth natural frequency of the duct.
It can be shown that ( ,0, ) ( ,0, ) x j x j and that
kk
j RQ where 1 j , and
( ,0, ) xj and
k
Q are real matrices. So, Green’s function in Eq. (4.21) can be written
as
78
1
1
( , , ) 2 ( ,0, ) (0, , )sin
( , ) 2 ( ,0, ) sin
k k b k k
k
k k k
k
x t x s j t
x t x j t
G Q M
HQ
(4.21)
Finally, substituting Eq. (4.21) into Eq. (4.19) yields the transient response of the duct
, , ,
1
( , ) 2 ( ,0, ) ( ) ( ) ( )
k f k b k o k
k
x t x j t t t
Q I I I (4.22)
where
,
00
,
0
, 0 0
2
0
0
( ) (0, , )sin ( ) ( , )
1
()
( ) sin ( )
()
0
1
( ) (0, , ) ( ) cos ( )sin
1
tL
f k b k k
t
L
b k k
R
L
o k b k k k k
t j t S d d
t t d
t j p t v t d
c
IM Φ
I
IM Φ
(4.23)
In this solution process, neither approximation has been made.
4.5 Specific Duct Systems
Determination of the state transition matrix is the key in the proposed solution method.
A state transition matrix of the system can be determined by the fundamental matrix in
Eq. (2.7). Generally, a fundamental matrix of acoustic ducts with arbitrary temperature
gradients cannot be always expressed in exact and closed-form. If it is not possible to
obtain exact solutions, a state transition matrix of the systems has to be evaluated
79
numerically by Eq. (2.26). In this section, exact fundamental matrices of ducts with
particular temperature gradients are presented.
By the theory of differential equations, a fundamental matrix is in the form
12
12
ˆ ˆ ( , ) ( , )
( , )
ˆ ˆ ( , ) ( , )
p x s p x s
xs
p x s p x s
xx
U (4.24)
where
1
ˆ ( , ) p x s and
2
ˆ ( , ) p x s
are two independent solutions of the s-domain differential
equation
22
2
ˆ ˆ ( , ) 1 ( , )
ˆ ( , ) 0
p x s dT p x s
p x s
x T dx x T
(4.25)
where /sR and T is a positive function of x. With
1
ˆ ( , ) p x s and
2
ˆ ( , ) p x s , the
transition matrix of the duct can be constructed by Eqs. (2.5) and (4.26).
In this work, we consider five types of distributions for the mean temperature () Tx .
Through use of suitable coordinate transformations (Hildebrand, 1964; Mei, 1997), the
independent solutions of Eq. (4.27) for these temperature profiles are derived, and they
are presented in sequel.
Type 1. Uniform:
0
() T x T (costant)
00
//
12
ˆ ˆ ,
x T x T
p e p e
(4.26)
Type 2. Linear:
0
() T x T mx
1 0 2 0
22
ˆ ˆ , p I T p K T
mm
(4.27)
80
Type 3. Exponential:
2
0
()
mx
T x T e
00
1 1 1 1
ˆ ˆ ,
TT
p I p K
TT
m T m T
(4.28)
Type 4. Conical:
2
0
( ) ( ) T x T mx
22
For 1 4 / 0 m ,
12
12
ˆ ˆ , ln p T p T T
(4.29a)
22
For 1 4 / 0 m ,
12
12
ˆ ˆ , p T p T
(4.29b)
where
12
and are the roots of the characteristic equation
2 2 2
/0 m .
Type 5. Power function:
0
( ) ( ) , 0
c
T x T mx c
Write (1 ) / (2 ) with 2 n c c c . For n to be an integer
1 2 1 2
2 2 2 2
12
22
ˆ ˆ ,
(2 ) (2 )
c c c c
c c c c
nn
p T I T p T K T
m c m c
(4.30a)
and for n to be a non-integer
1 2 1 2
2 2 2 2
12
22
ˆ ˆ ,
(2 ) (2 )
c c c c
c c c c
nn
p T I T p T I T
m c m c
(4.30b)
where ,
nn
IK are modified Bessel functions of order n. When 2 c , the solutions are
given in Type 4.
The solutions presented in this section cover a variety of temperature distributions of
duct systems in engineering problems.
81
4.6 Acoustic Control of Duct Systems
4.6.1 Closed-Loop System
Consider the closed-loop duct system whose response is controlled by one point
actuator (speaker) located at
a
x based on the measurements from r point sensors
(microphones) located at
l
s
x , 1,2, , lr . In each sensor loop, there is a time delay
block of gain
l
and time delay constant 0
l
d
T . The controller has a transfer function of
the form ( ) ( ) / ( ) C s N s D s where the polynomials N(s) and D(s) have no common
roots, and μ is a gain parameter. To guarantee the controllability and observability for all
modes, the actuator and sensors must avoid the modal points of the system eigenfunctions
and the controller must avoid pole-zero cancellation (Yang, 1991). In the closed-loop
control, the control source located at
a
x is given as
1
ˆ
ˆ ( , ) ( ) ( ) ( , )
d
l
l
r
Ts
c a l s
l
S x s x x C s e p x s
(4.31)
By Eq. (4.12), the acoustic pressure of the closed-loop duct system in s-domain is
0
ˆ ˆ
ˆ ˆ ( , ) ( , , ) ( , ) ( , )
L
o c eI
p x s g x s S s S s d
(4.32)
where ˆ ( , , )
o
g x s is the 1-row and 2-column element of the distributed transfer function
Ĝ given in closed form, and can be rewritten as
( , , )
ˆ ( , , )
()
o
o
o
K x s
g x s
Ds
(4.33)
82
The denominator ()
o
Ds is the characteristic equation of an open-loop duct system, i.e.
det ( , , )
bb
Ls MN .
Substitution of Eq. (4.31) into Eq. (4.32), solving for ˆ ( , )
l
s
p x s by letting
l
s
xx and
replacing ˆ ( , )
l
s
p x s by the solution gives the system response (Yang, 1991)
0
ˆ ˆ ( , ) ( , , ) ( , )
L
cl eI
p x s g x s S s d
(4.34a)
where
( , , )
ˆ ( , , )
()
cl
cl
cl
K x s
g x s
Ds
(4.34b)
1
1
()
( , , ) ( ) ( , , ) { ( , , ) ( , , )
()
( , , ) ( , , ) }
d
l
l
d
l
l
r
Ts
cl o o a o s l
l
o
r
Ts
o o s a l
l
Ns
K x s D s K x s K x x s K x s e
Ds
K x s K x x s e
(4.34c)
1
( ) ( ) ( ) ( ) ( , , )
d
l
l
r
Ts
cl o o s a l
l
D s D s D s N s K x x s e
The characteristic equation of the control system is ( ) 0
cl
Ds .
4.6.2 Time-Delay Relations and Stability
In Ref. (Yang, 1992), it is shown that the system response at one point can be exactly
predicted from the measurement at other point(s) of the system and the destabilizing
effect of non-collocation of sensors and actuators is eliminated through introduction of
specific time delay block(s) in the control system. The time delay relations of ducts with
83
axial temperature gradients are represented for two different boundary conditions as
follows.
CASE A: Two Open Ends
(a) For a duct with uniform temperature (Type 1), the time delay relation for Eq. (4.6) is
given by (Yang, 1992)
( ) ( )
kd
jT
k a k s
v x v x e
with , , (0, ) and
s a s a d
L
x x L x x L T
c
(4.35)
In this case, Eq. (4.6) of the duct is same as that of uniform second order systems such as
strings, bars and shafts (Meirovitch, 2001).
(b) Define the temperatures
1
TT at x = 0,
2
TT at x = L,
a
TT at
a
xx ,
s
TT at
s
xx and , (0, )
sa
x x L . For a duct with a linear temperature profile (Type 2), the time
delay relation of Eq. (4.6) is given by
( ) ~ ( )
kd
jT
k a k s
v x v x e
with
12
21
and
a s d
TT
T T T T T
a
(4.36)
where
2
m
aR . See Appendix B for derivations.
(c) For a duct with a conical temperature profile (Type 4), the time delay relation for Eq.
(4.6) is given by
( ) ~ ( )
kd
jT
k a k s
v x v x e
with
12
21
ln( / )
and
a s d
TT
T T T T T
m
(4.37)
See Appendix C for details.
84
CASE B: Closed-Open End
(a) For a duct with a uniform temperature (Type 1), the time delay relation for Eq. (4.6) is
given by (Yang, 1992)
3 12
1 2 3
( ) ( ) ( ) ( )
kd k d k d
jT j T j T
k a k s k s k s
v x v x e v x e v x e
(4.38a)
where
1 2 3 1 2 3 1 2 3
, , , , , (0, )
a s s s s a s s a s s s
x x x x x x x x x x x x L and
2 3 1 3 3
1 2 3
,,
s s s s a s
ddd
x x x x x x
TTT
c c c
(4.38b)
(b) Define the temperatures
1
s
TT at
1
s
xx ,
2
s
TT at
2
s
xx and
3
s
TT at
3
s
xx . A
duct with a linear temperature profile (Type 2) has the time delay relation of Eq. (4.6) as
3 12
1 2 3
( ) ~ ( ) ( ) ( )
kd k d k d
jT j T j T
k a k s k s k s
v x v x e v x e v x e
(4.39a)
where
1 2 3 1 2 3 1 2 3
, , , , , (0, )
a s s s s a s s a s s s
T T T T x x x x x x x x L and the time
delays are
2 3 1 3 3
1 2 3
,,
s s s s a s
ddd
T T T T T T
TTT
a a a
(4.39b)
See Appendix D for derivations.
(c) For a duct with a conical temperature profile (Type 4), the time delay relation for Eq.
(4.6) is given by
3 12
1 2 3
( ) ~ ( ) ( ) ( )
kd k d k d
jT j T j T
k a k s k s k s
v x v x e v x e v x e
(4.40a)
where
2 1 3 1 2 3 1 2 3
, , , , , (0, )
a s s s s a s s a s s s
T T T T x x x x x x x x L and the time delays are
2 3 1 3 3
1 2 3
ln( / ) ln( / ) ln( / )
,,
s s s s a s
ddd
T T T T T T
TTT
m m m
(4.40b)
85
See Appendix E for derivations.
From Eq. (4.34c) and ( ) 0
cl
Ds it is seen that the closed-loop eigenvalues depend on
the gain parameter μ which is real and positive. Based on the analysis by Yang and Mote
(1991, 1992), the control system (4.34) is asymptotically stable if
0
Re 0
k
ds
d
for all k, and for (0, )
cr
(4.41)
where ()
k
s is the kth root locus of the closed-loop system originating from the open-
loop poles
k
j at 0 . According to ref. (Yang and Mote, 1991) and Eq. (4.34b, c),
1
0
( ) ( , , )
( ) ( )
kd
l
l
k
r
jT
k o s a k l
kl
ok
sj
N j K x x j e
ds
d d
D s D j
ds
(4.42)
Let the intersection with the imaginary axis be
k
jb . The critical value
cr
is given by
1
()
inf
( ) ( , , )
kd
l
l
ok
cr r
k
jb T
k o s a k l
l
D jb
C jb K x x jb e
(4.43)
where
cr
is the first intersection of the kth root locus with the imaginary axis with
increasing μ .
4.6.3 Perfect Noise Cancelling Controller
The objective of noise control is to attenuate sound at a specific (reference)
location()
r
xx . Assume zero initial conditions and a point source. By Eq. (4.34), the
acoustic pressure of the closed-loop system in the s-domain is given by
86
ˆ ˆ ( , ) ( , , ) ( , )
r cl r
p x s g x s S s (4.44a)
where
( , , )
ˆ ( , , )
()
cl r
cl r
cl
K x s
g x s
Ds
(4.44b)
ˆ ˆ ˆ ˆ ˆ ( , , ) ( , , ) ( ){ ( , , ) ( , , ) ( , , ) ( , , )}
ˆ ( ) 1 ( ) ( , , )
cl r o r o r a o s o r o s a
cl o s a
K x s g x s C s g x x s g x s g x s g x x s
D s C s g x x s
The output ˆ ( , )
r
p x s is zero if the noise is perfectly cancelled. Letting ˆ ( , )
r
p x s be zero in
Eq. (4.44), the transfer function of a controller can be determined by ( , , ) 0
cl r
K x s .
ˆ ( , , )
()
ˆ ˆ ˆ ˆ ( , , ) ( , , ) ( , , ) ( , , )
or
o r o s a o r a o s
g x s
Cs
g x s g x x s g x x s g x s
(4.45)
As an example, a duct with uniform temperature is closed at the left end ( /0 px at
x = 0) and open at the right end ( 0 p at x = L). Consider the time delay control design
in Figure 4.1(c). By the DTFM, the transfer function of an open-loop system is given by
sinh ( )cosh
,
cosh
ˆ ( , , )
cosh sinh ( )
,
cosh
o
Lx
x
L
g x s
xL
x
L
(4.46)
where / sc .
The external point source is located at x = 0. By Eq. (4.45), the transfer function of a
controller is given by
1
1
( ) ( ) 2
22
()
(1 )
a s a s
sT
x x x x sT
s
C s e
e e c e
(4.47)
where
s m r
x x x .
87
The stability of a closed-loop system can be examined by Eq. (4.41). Let the open-loop
poles be
k
j . The transfer function of a controller can be written as
2
()
sin ( )
k
k
k
as
Cj
c x x
c
(4.48)
where μ is real and positive.
From Eq. (4.48) it is seen that the transfer function of a controller is always real.
1
2
0
cosh sinh ( )
2
sin ( ) sinh
21
cos cos
sin ( )
k
kk
sa
jT kk
kk
a s k
kk
sa
k
as
jj
c x L x
ds
cc
e
j
d
c x x j L L
cc
c
j x x
L c c
xx
c
(4.49)
Thus, we know that the closed-loop system is marginally stable.
0
Re 0
k
ds
d
for all k, and for (0, ) (4.50)
From this fact, the closed-loop system is marginally stable if the transfer function of a
controller has only real part. It is difficult for a perfect noise cancelling controller to be
implemented practically because it is infinite dimensional. Through use of finite
approximation of the controller, a finite dimensional controller can be determined for
practical applications of active noise control in a duct (Pota and Kelkar, 2000). Similarly,
a finite dimensional controller for ducts with non-uniform temperature gradients can be
obtained by the proposed method. The closed-loop system can be marginally stable when
the transfer function of the controller has only real part.
88
C(s)
x
x
s x
a
d
Ts
e
(a)
C(s)
Source 2
Source 1
x
x
a
x
s3
1 d
Ts
e
2 d
Ts
e
3 d
Ts
e
x
s2 x
s1
+
-
+
Time delay
(b)
C(s)
x
x
s x
a
(c)
x
r
Figure 4.1 Schematic of noncollocated control of ducts
89
4.7 Examples
Consider a dry-air flow duct of finite length, which is closed at the left end (x = 0) and
open at the right end (x = L). In the numerical simulation, three types of temperature
distributions in a duct are considered:
(a) Uniform distribution (Type 1):
00
, with 300K. T T T
(b) Linear distribution (Type 2):
00
, with 700K and 50. T T mx T m
(c) Exponential distribution (Type 3):
2
00
, with 700K and 0.053.
mx
T T e T m
The duct length is assumed to be L= 8m. The boundary conditions of the duct are
0 : / 0
:0
x p x
x L p
(4.51)
For each temperature distribution, the state matrix in Eq. (4.11) is
2
2
01
For uniform distribution,
0
s
c
F (4.52a)
2
2
0
01
For linear distribution,
sm
c T mx
F (4.52b)
2
2
01
For exponential distribution,
2
s
m
c
F (4.52c)
The boundary matrices, according to Eq. (4.34), are of the form
0 1 0 0
,
0 0 1 0
bb
MN (4.53)
90
Examples 4.1: Eigensolutions of the ducts of three types
The eigenfrequencies of the duct with the above-mentioned three temperature
distributions are computed by using Eq. (4.17) and the first six eigenfrequencies are listed
in Table 4.1. It is seen that the temperature distribution can significantly alter the system
eigenvalues. According to Eq. (4.19), the relative acoustic pressure amplitude of the duct
can be computed as follows
1
1
( , , ) ()
(0) (0, , )
xj Px
Pj
a
a
e
e
(4.54)
where is an eigenfrequency of the duct, and
*
0 x (the closed end of the duct) has
been used. Consider the sixth eigenfrequency in each case of the mean temperature
distributions; namely,
6
= 750.14 rad/s for the uniform distribution,
6
= 947.38 rad/s
for the linear distribution, and
6
= 919.62 rad/s for the exponential distribution. Figure
4.2 plots the relative acoustic pressure amplitude for the selected eigenfrequencies. The
acoustic pressure amplitudes with uniform and exponential distributions match the results
given in Ref. (Sujith et al., 1995).
Example 4.2: Transient acoustic pressure due to boundary excitation
Consider the duct with the linear temperature distribution. A sinusoidal wave is applied
at the left end of the duct
0
0
( , )
sin( )
Lc
x
p x t
t
x
(4.55)
91
This boundary condition can be used to describe sound generation by a speaker in an
enclosed cabinet. By Eqs. (4.24) and (4.25), the transient response of the duct is
0
1
0
1
( , ) 2 ( ,0, ) sin ( )sin( )
0
t
L k k c
k
x t x j t d
Q (4.56)
Assign
0
0.02
L
Pa/m and 947
c
rad/s. The transient acoustic pressure under the
boundary excitation is computed by Eq. (4.24) with the first 60 terms (modes). The
spatial profiles of the acoustic pressure at times t = 0.004, 0.008, 0.012, 0.018, 0.03 and
0.2s are plotted in Figure 4.3, which shows the propagation of sound waves along the
duct. Because the excitation frequency
c
is close to the sixth eigenfrequency of the duct
(see Table 4.1), the pattern of the waves is similar to the sixth mode shape
(eigenfunction) of the duct as shown in Figure 4.2.
Example 4.3: Transient acoustic pressure subject to an excitation at an interior point
For the duct with the exponential excitation, sound waves are generated by a point
source at
f
xx ,
/
0
( , ) (1 ) ( )
t
f
S x t S e x x
. By Eq. (4.7), ( , ) S x t is related external
forces and mass influx. With Eq. (4.24), the transient response of the duct is
/
0
1
0
0
( , ) 2 ( ,0, ) (0, , ) sin ( ) (1 )
1
t
k f k k
k
x t S x j x j t e d
QM
b
(4.57)
For simulation, choose
0
0.02 S kgm
-3
sec
-2
, 10sec and 4
f
xm , and the first 60
terms (modes) in Eq. (4.23). Figure 4.4 shows the spatial profiles of the acoustic pressure
92
at times t = 0, 0.05, 0.1, 0.15 and 0.2 sec. In Fig. 4.5, the acoustic pressure at the closed
end (x = 0) is plotted against time.
Example 4.4: Active noise control of ducts with two types of boundaries
(1) A duct with two open ends
A duct is open at the left end ( 0 p at x = 0) and open at the right end ( 0 p at x =
L). Consider the time delay control design in Figure 4.1(a). One sensor and one actuator
are located at
s
x and
a
x , respectively. For a duct with a uniform temperature, the transfer
function of an open-loop system is given by Eqs. (4.13) and (4.33)
( , , )
ˆ ( , , )
()
o
o
o
K x s
g x s
Ds
(4.58a)
1
sinh ( )sinh ,
( , , )
1
sinh sinh ( ),
( ) sinh
o
o
L x x
K x s
x L x
D s L
(4.58b)
where / sc .
Let the open-loop poles be
k
j . By the time delay relation Eq. (4.35) and Eq. (4.42),
the stability of a closed-loop system is investigated with the geometrical property of root
loci (Yang and Mote, 1991)
20
0
0
()
sin
()
kk
ka
kk
ds N j c
j p x
d p L D j
(4.59)
93
where
0
/
kk
pc . Note that c is constant for a duct with uniform temperature the sin
2
term in Eq. (4.59) is non-zero and positive.
By Eq. (4.41), it is seen that the control system under noncollocated sensor and
actuator with a specific time delay, is stable for (0, )
cr
if for (0, )
()
Im 0
()
Nj
Dj
(4.60)
Let the intersection with the imaginary axis be
k
jb . By Eq. (4.43), the critical gain
value
cr
is found be
22
()
inf min( , )
( )sin
kk
cr s a
k
k
k
b I b
a x x
b
cM b a
c
(4.61a)
where
()
2 2 1
( ) ( ) ( ) ( )
()
( ) ( ) ( ), ( ) tan
()
j
C j R jI M e
I
M R I
R
(4.61b)
For a duct with a linear temperature distribution
0
() T x T mx , the transfer function
of an open-loop system is given by Eqs. Eq. (4.35) and (4.42). Let the open-loop poles
be
k
j . With the time delay relation Eq. (4.36) and Eq. (4.41), the stability of a closed-
loop system is investigated with the geometrical property of root loci as shown in
Appendix F.
3/4
2 2
1/4
0
21
()
sin
()
a k a k
k
kk
s
TT ds aT N j j
d D j a
T T T
(4.62)
94
where
0
/
kk
pc . Note that the sin
2
term in Eq. (47) is non-zero and positive.
By Eq. (4.42), it is seen that the control system under noncollocated sensor and
actuator with a specific time delay, is stable for (0, )
cr
if for (0, )
()
Im 0
()
Nj
Dj
(4.63)
In numerical simulation, the locations of the sensor and actuator are assigned.
For uniform temperature distribution: 0.12 ,
s a s
x L x L x
For linear temperature distribution:
2
2 1 0
0.12 , {( ) }/
s a s
x L x T T T T m
In a certain region of the parameter μ (0 0.01) and under a controller () C s s , the
root loci of the first four eigenfrequencies of the control systems are drawn in Fig. 4.6(a).
(2) A duct with closed-open end
A duct is closed at the left end ( /0 px at x = 0) and open at the right end ( 0 p
at x = L): see Figure 4.1(c). Three sensors (microphones) and one actuator (speaker) are
located at
1 2 3
,,
s s s
x x x and
a
x , respectively. The transfer function
( , , )
ˆ ( , , )
()
o
o
o
K x s
g x s
Ds
of
an open-loop system is given by
1
sinh ( )cosh ,
( , , )
1
cosh sinh ( ),
( ) cosh
o
o
L x x
K x s
x L x
D s L
(4.64)
where / sc .
95
Noncollation of sensors and actuators often face the destabilizing effect of system.
With the specific time delay block(s) in the control system, the destabilizing effect can be
abbreviated. By the time delay relation Eq. (4.38) and Eq. (4.41), the stability of a closed-
loop system is investigated with the geometrical property of root loci (Yang and Mote,
1991)
1
1
3 2
23
00
00
0
0 0 0 0
2
20
()
{sin ( )cos
( )sin
cos sin ( ) cos sin ( ) }
()
cos
()
kd
kd kd
jT
kk
k s k a
k k k
jT jT
k s k a k s k a
k
ka
kk
ds N j c
p L x p x e
d jp L D j p L
p x p L x e p x p L x e
Nj c
j p x
L D j
(4.65)
where
0
/
kk
pc . Note that the cos
2
term in Eq. (4.65) is non-zero and positive.
By Eq. (4.42), it is seen that the control system under noncollocated sensor and
actuator with a specific time delay, is stable for (0, )
cr
if for (0, )
()
Im 0
()
Nj
Dj
(4.66)
The parameter μ is not arbitrarily large for the noncollocated sensor and actuator. The
acoustic pressure in a duct will be destabilized if
cr
. By Eq. (4.43), the critical gain
parameter
cr
is determined (Yang and Mote, 1991). With the same procedures of a duct
with a uniform temperature, the stability condition of a duct with a linear temperature
distribution,
0
() T x T mx , can be determined. In a duct with non-uniform temperature
profiles, a critical gain parameter
cr
is much influenced by the temperatures at the
location of sensor and actuator.
96
For a duct with a uniform temperature, the transfer function of a controller is given by
() C s s . The locations of sensors and an actuator are given as follows:
1 2 3
1 2 3
0.92 , 0.24 , 0.12
s s s
a s s s
x L x L x L
x x x x
(4.67)
For a duct with a linear temperature distribution, i.e.
0
() T x T mx , the locations of
sensors are same as Eq. (4.67) and that of an actuator is determined by
1 2 3
2
0 s s s
a
T T T T
x
m
(4.68)
where
1
s
T ,
2
s
T and
3
s
T are the temperatures at the points
1 2 3
, and
s s s
x x x x , respectively.
According to Eqs. (4.41) and (4.43), the duct systems are stable for (0, )
cr
. The
root loci of the first four eigenfrequencies of the control systems in a certain region of the
parameter μ, i.e. 0 0.01 , are drawn in Fig. 4.6(b). With the stability verified,
consider the frequency response of the duct to the harmonic excitation at the closed end.
The response is computed at the point of the first sensor
1
()
s
xx . Figure 4.7 shows the
frequency response of a duct with a linear temperature distribution under the controller
and no controller. The resonance peak amplitudes are reduced by the proper controller.
97
Table 4.1 The eigenfrequencies
k
(rad/s) of the duct with three temperature distributions
Temperature distribution 1
2
3
4
5
6
Uniform distribution 68.19 204.58 340.97 477.36 613.75 750.14
Linear distribution 78.98 256.26 429.54 602.30 774.88 947.38
Exponential distribution 76.66 248.75 416.96 584.65 752.18 919.62
Figure 4.2 Relative acoustic pressure amplitude at the sixth eigenfrequency;
dot-line −: uniform temperature distribution, dash-line −: exponential
temperature distribution, solid-line −: linear distribution
98
(a) t = 0.004s (b) t = 0.008s
(c) t = 0.012s (d) t = 0.018s
(e) t = 0.03s (f) t = 0.2s
Figure 4.3 Transient waves in the duct due to a boundary excitation at the left end (x=L)
99
Figure 4.4 Spatial distributions of the acoustic pressure due to internal sound generation
Figure 4.5 Time history for internal sound generation at the closed end (x = 0)
100
Figure 4.6 Root loci of the duct under the noncollocated controller; (a) two open ends,
(b) closed-open end; o: uniform temperature profile; *: linear temperature profile
Figure 4.7 The frequency response of a duct with linear temperature profile;
dot-line −: without control; solid-line −: with the time-delayed control
101
Chapter 5. Dynamic Analysis of Transversely Vibrating Beams
5.1 Introduction
In this chapter, we present an analytical method in a systematic approach to obtain an
exact and closed form transient solution of a beam. Two models for the transversely
vibrating uniform beam are considered. The two beam models are the Euler-Bernoulli
and Timoshenko. The proposed method is valid for both self-adjoint and non-self-adjoint
systems, and is convenient for computer coding. The high efficiency of the method is
justified in the various numerical examples.
5.2 System Descriptions
5.2.1 Euler-Bernoulli beam model
The transverse displacement of a uniform Euler-Bernoulli beam is governed by
24
24
( , ) ( , )
( , ), (0, )
w x t w x t
EI f x t x L
tx
(5.1)
where w(x,t) is the transverse displacement,
the linear density and EI the bending
stiffness. The boundary conditions of the system are of the general form
23
0 1 ,1 23
2
2 3 ,2 2
(0, ) (0, )
at 0 : (0, ) ( )
(0, ) (0, )
()
LL
L
w t w t
x m a w t a t
tx
w t w t
a a t
xx
(5.2a)
102
23
0 1 ,1 23
2
2 3 ,2 2
( , ) ( , )
at : ( , ) ( )
( , ) ( , )
()
RR
R
w L t w L t
x L m b w L t b t
tx
w L t w L t
b b t
xx
(5.2b)
where
i
a and ( 0,1,2,3)
i
bi are constants that are properly assigned to characterize
different types of boundaries;
L
m and
R
m are end masses, and
,1 ,2 ,1
( ), ( ), ( )
L L R
t t t and
,2
()
R
t
are related to prescribed boundary excitations. The system is subject to the initial
conditions
00
( ,0)
( ,0) ( ), ( ), (0, )
wx
w x u x v x x L
t
(5.3)
where
0
() ux and
0
() vx are given profiles of initial displacement and velocity of the
system.
5.2.2 Timoshenko beam model
Consider a beam with attached lateral springs ,
LR
kk and torsional springs ,
LR
;
see Fig. 5.1. Based on Timoshenko beam theory, the transverse deflection and the angle
of rotation due to bending are denoted by ( , ) w x t and ( , ) xt , respectively. The kinetic,
potential energy and the potential work done by the external force ( , ) f x t are given by
2 2 2 2 22
0
2
2 2 2
0
22
1 1 1 1 1
(0, ) (0, ) ( , ) ( , ) ( , ) ( , )
2 2 2 2 2
1 1 1
( , ) ( , ) ( , ) (0, ) (0, )
2 2 2
11
( , ) ( , )
22
L
L L R R
L
s L L
RR
T dx m w t I t m w L t I L t w x t I x t
V EI x t k w x t x t dx k w t t
k w L t L t
103
0
( , )
L
W f x t wdx
(5.4)
where
L
m and
R
m are end masses;
L
I and
R
I are rotational moment of inertia; I
the
rotary inertia and
s
k the shear stiffness; the dot (∙) denotes the derivative with respect to
time and the prime (ʹ) the derivative with respect to spatial coordinate x. The shear
stiffness is denoted by
s
k AG where G is the shear modulus, A the cross-section
area, and a “shape factor” that depends on the geometry of the cross-section. The value
2 / 3
is often associated with rectangular cross-sections, while the value 3/ 4
is
similarly associated with circular cross-sections (Bottega, 2006).
Introducing Eq. (5.4) into the extended Hamilton’s principle
2
1
( ) 0
t
t
T V W dt
(5.5)
yields
22
11
22
11
22
11
00
0
0
( ) ( , ) ( )
( ) ( )
tt LL
ss
tt
tt
s R R s L L
x L x
tt
tt
R R L L
x L x
tt
k w w f x t wdxdt EI k w I dxdt
xx
k w k w m w wdt k w k w m w wdt
EI I dt EI I
0 dt
(5.6)
From Eq. (5.6), it then follows that the governing equations of motions for the vibration
of the uniform Timoshenko beam are
22
22
22
22
( , ) ( , ) ( , )
( , ), (0, )
( , ) ( , ) ( , )
( , ) 0
s
s
w x t x t w x t
k f x t x L
t x x
x t x t w x t
I EI k x t
t x x
(5.7)
104
Equation (5.7) is a model for Timoshenko beam which adds the effects of shear
distortion and rotary inertia to the Euler-Bernoulli beam. From the last four terns of Eq.
(5.6), the boundary conditions of the beam are of the general form
2
0 1 ,1 2
2
2 3 ,2 2
(0, ) (0, )
at 0 : (0, ) (0, ) ( )
(0, ) (0, )
(0, ) ( )
LL
LL
w t w t
x m a w t a t t
tx
tt
I a t a t
tx
(5.8a)
2
0 1 ,1 2
2
2 3 ,2 2
( , ) ( , )
at : ( , ) ( , ) ( )
( , ) ( , )
( , ) ( )
RR
RR
w L t w L t
x L m b w L t b L t t
tx
L t L t
I b L t b t
tx
(5.8b)
where
,1 ,2 ,1
( ), ( ), ( )
L L R
t t t and
,2
()
R
t are related to prescribed boundary excitations
(external load/ moment or possible foundation motion);
i
a and ( 0,1,2,3)
i
bi are
constants that are properly assigned to characterize different types of boundaries.
The system is subject to the initial conditions
00
00
( ,0)
( ,0) ( ), ( ), (0, )
(0, )
(0, ) ( ), ( )
wx
w x u x v x x L
t
t
t x x
t
(5.9)
where
0 0 0
( ), ( ), ( ), u x v x x and
0
() x are prescribed functions that describe the initial
displacement, velocity, rotation, and rotation rate.
105
Figure 5.1 Schematic of a beam model
5.3 Distributed Transfer Function Formulation
5.3.1 Euler-Bernoulli Beam
Taking Laplace transform to Eqs. (5.1)- (5.3) with respect to time gives
4
2
4
( , )
( , ) ( , ), (0, )
eI
w x s
s w x s EI p x s x L
x
(5.10)
3
2
1 0 ,1
3
2
2 3 ,2
2
(0, )
(0, ) ( )
(0, ) (0, )
()
LL
L
ws
a m s a w s s
x
w s w s
a a s
x x
(5.11a)
3
2
1 0 ,1
3
2
2 3 ,2
2
( , )
( , ) ( )
( , ) ( , )
()
RR
R
w L s
b m s b w L s s
x
w L s w L s
b b s
x x
(5.11b)
where the over-bar stands for Laplace transform, s is the Laplace transform parameter,
and
00
( , ) ( , ) ( ) ( )
eI
p x s f x s su x v x
106
,1 0 0 ,2 ,2
,1 0 0 ,2 ,2
( ) ( ) (0) (0) , ( ) ( )
( ) ( ) ( ) ( ) , ( ) ( )
L L L L L
R R R R R
s s m su v s s
s s m su L v L s s
(5.12)
By defining the spatial state vector
23
23
( , ) ( , ) ( , )
ˆ ( , ) ( , ) , 0,
T
w x s w x s w x s
x s w x s x L
x xx
(5.13)
Eqs. (5.10) - (5.12) are rewritten in a first-order state form
ˆ ˆ ˆ ( , ) ( ) ( , ) ( , ), (0, )
ˆ ˆ (0, ) ( , ) ( )
b b b
x s s x s x s x L
x
s L s s
Fp
MN
(5.14)
where
2
0 1 0 0
0
0 0 1 0
0
( , )
ˆ ( ) , ( , )
0 0 0 1
0
1
000
eI
p x s
s x s
EI
s
EI
Fp (5.15)
2
,1
01
,2
23
2
,1 01
,2 23
() 0 0 0 0 00
() 0 0 0 0 00
, , ( )
() 00 0 0 0 0
() 00 0 0 0 0
L
L
L
b b b
R R
R
s
m s a a
s
aa
s
s m s b b
s bb
MN
5.3.2 Timoshenko Beam
By taking Laplace transform to Eqs. (5.7) - (5.9) with respect to time,
2
2
,1 2
2
2
,2 2
( , ) ( , )
( , ) ( , ), (0, )
( , ) ( , )
( , ) ( , ) ( , )
s eI
s eI
x s w x s
s w x s k p x s x L
xx
x s w x s
I s x s EI k x s p x s
xx
(5.16)
107
2
0 1 ,1
2
2 3 ,2
(0, )
(0, ) (0, ) ( )
(0, )
(0, ) ( )
LL
LL
ws
m s a w s a s s
x
s
I s a s a s
x
(5.17a)
2
0 1 ,1
2
2 3 ,2
( , )
( , ) ( , ) ( )
( , )
( , ) ( )
RR
RR
w L s
m s b w L s b L s s
x
Ls
I s b L s b s
x
(5.17b)
and
,1 0 0
,2 0 0
( , ) ( , ) ( ) ( )
( , ) ( ) ( )
eI
eI
p x s f x s su x v x
p x s I s x x
,1 ,1 0 0
,2 ,2 0 0
( ) ( ) (0) (0)
( ) ( ) (0) (0)
L L L
L L L
s s m su v
s s I s
(5.18)
,1 ,1 0 0
,2 ,2 0 0
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
R R R
R R R
s s m su L v L
s s I s L L
By defining the spatial state vector
( , ) ( , )
ˆ ( , ) ( , ) ( , ) , 0,
T
w x s x s
x s w x s x s x L
xx
(5.19)
Eqs. (5.16) - (5.18) is rewritten in a first-order state form
ˆ ˆ ˆ ( , ) ( ) ( , ) ( , ), (0, )
ˆ ˆ (0, ) ( , ) ( )
b b b
x s s x s x s x L
x
s L s s
Fp
MN
(5.20)
where
108
2
,1 ,2
2
0 0 1 0
0 0 0 1 00
00
ˆ ( ) , ( , ) 0 0 1
10
01
00
eI eI
s s
s
s
pp
s
s x s
kI k
I s k
k
EI EI
Fp (5.21)
2
0 1 1
2
23
2
0 1 1
2
23
0 0 0 0 0
0 0 0 0 00
,
0 0 0 0 0
00 0 0 0 0
L
L
bb
R
R
m s a a a
I s a a
m s b b b
I s b b
MN
with
,1 ,2 ,1 ,2
( ) ( ) ( ) ( ) ( )
T
b L L R R
s s s s s .
In this section, the state equation and boundary conditions of Euler-Bernoulli beam and
Timoshenko beam are formulated in matrix forms.
5.4 Dynamic Analysis
5.4.1 System Response
The s-domain solution of Eqs. (5.14) and (5.20) is given as
0
ˆ ˆ
ˆ ˆ ( , ) ( , , ) ( , ) ( , ) ( ), , (0, )
L
b
x s x s x s d x s s x L
G p H (5.22)
where
ˆ
G and
ˆ
H are the distributed transfer functions of the system,
()
( )(
1
( ) ( )
ˆ
( , ) ,
ˆ
( , , )
ˆ
( , ) ,
ˆ
( , )
s
b
sL
b
s x s L
bb
x s e x
xs
x s e x
x s e e
F
F
FF
HM
G
HN
H M N
(5.23)
109
Assume that the system is subject to harmonic excitations. The steady-state response
(the frequency response) of the system, by Eq. (5.22), is given by
( , ) ( , e , (0, )
jt
ss
x t x j x L
(5.24)
where
0
0
()
ˆ ˆ
( , ( , , ) d ( , ) , (0, )
1 ()
L
L
R
F
x j x j x j x L
GH (5.25)
and
ˆ
( , , ) xj G and
ˆ
( , ) xj H are the frequency response functions of the system. With
Eq. (5.23), the frequency response functions can be evaluated in the exact and closed
form.
The transient solution of the system is obtained through inverse Laplace transform of
Eq. (5.22), which produces the following Green’s function formula
0 0 0
( , ) ( , , ) ( , )d d ( , ) ( ) d
t L t
b
x t x t x t
G p H (5.26a)
where the matrix Green’s functions G and H are the inverse Laplace transforms of the
transfer functions
ˆ
G and
ˆ
H , respectively; ( , ) xt p and ()
b
t are the inverse Laplace
transforms of ˆ( , ) xs p and ()
b
s ;and ( , ) xt is the inverse Laplace transform of ˆ ( , ) xs .
By a residue formula, the Green’s functions of the system are determined as
( ) ( )
1
()
1
( , , ) 2 Re e
( , ) 2 Re e
k k k
kk
s x s s t
kb
k
s x s t
k
k
x t e e
x t e
FF
F
G R M
HR
(5.26b)
where
k
R is the residue of distributed transfer function.
110
5.4.2 Eigenvalue Problem
By vanishing disturbance terms in Eqs. (5.14) and (5.20), the eigenvalue problem is
described by
()
( , ) ( ), (0, )
x
x s x x L
x
F
(5.27a)
(0) ( ) 0
bb
L MN (5.27b)
where s is an eigenvalue, and () x the eigenfunction associated with s.
The solution of Eq. (5.27a) has the form
()
()
sx
xe
F
a where a is a constant vector to
be determined. Substitution it into Eq. (5.27b) yields
()
det ( ) det 0
sL
bb
se
F
Z M N (5.28)
Eq. (5.28) is the characteristic equation of and its roots are the eigenvalues of the system.
In the basic kinematic assumptions for Euler-Bernoulli beam theory it is assumed that
shearing and rotary inertia can be neglected. However, the importance of shearing and
rotary inertia increases for higher order resonant modes. Euler beam theory has lower
effective mass and is stiffer than that of Timoshenko beam. Thus, the eigenvalues of
Timoshenko beam are lower than those of Euler-Bernoulli beam in Fig. 5.2.
The radius of gyration is related to the rotary inertia and given as /.
G
r I A If it is
assumed that the rotary inertia is small ( 0
G
r ), and then in the limit of G the
shearing becomes negligible, Timoshenko beam theory are identical to Euler beam theory.
It means that if
G
r ( / )
G
rL is small, Euler beam theory can be used. Conversely,
Timoshenko beam theory should be used.
111
Shearing and rotary inertia become increasingly important for higher mode and thus
Euler beam theory fails to describe the eigenvalues of the higher modes as shown in Fig.
5.2. It is shown that Timoshenko beam theory can be used for / 0.1 hL and Euler beam
theory can be used for / 0.05 hL (Han et al., 1999) where h is the height of a beam.
Example 5.1: A cantilever beam
A beam in transverse vibration is fixed at the left end, and free at the right end. The
vibration of an Euler-Bernoulli beam and a Timoshenko beam is governed by Eqs. (5.1)
and (5.7), respectively. By Eq. (5.28), the natural frequencies are the non-negative roots
of the real-valued transcendental equation
()
( ) det( )
jL
bb
e
F
MN (5.29)
with sj and 1 j .
For numerical simulation, the physical parameters of the beam are as follows: 1,
4
0.02, 10, 10 and 1
s
I EI k L
. The first 5 natural frequencies and several higher-
mode frequencies of Timoshenko beam are listed in Table 5.1, where the proposed
method in compared with the finite element method (FEM). The first 4 modes of the
vibration of the systems obtained by
()
()
sx
xe
F
a are plotted in Fig. 5.3. It is shown that
the eigenfunctions of Euler-Bernoulli and Timoshenko beam models are quite different as
the mode increases. Let the beams be initially at rest. Consider an external load
0
( ) (1 )
t
f t F e
applied to the beam at the point
f
x . Thus, the external force term in
Eq. (5.26a) is
112
0
3
1
( , )
t
s
Fe
xt
k
p e (5.30)
where
T
3
0 0 1 0 e . Assume that
0
0.2, 1.0 and / 0.5
f
F x L . The forced
vibration of beam is computed with first 10 terms in Eq. (5.26) and plotted in Fig. 5.4 and
5.5. Figure 5.4 shows the spatial distributions of the displacements and the rotation of the
beam at t = 0.1, 0.3, 0.5, 0.7, and 1.0. Figure 5.5 gives a time history (0 10) t of the
vibration at the location (x = L), where compared with the solution by Euler-Bernoulli
beam model.
113
Figure 5.2 The first natural frequency curves of the two models of a beam with
the modulus ratio / 0.1 E kG and
4
0
/ EI L ;
solid-line −: Euler-Bernoulli beam, dash-line −: Timoshenko beam
Table 5.1 The natural frequencies
k
of the Timoshenko beam (in rad/s)
k
DTFM
(Timoshenko)
Finite element method
30 elements 150 elements 300 elements
1
2
3
4
5
10
15
20
10.61197
53.76193
122.43730
196.37605
272.09203
636.23061
878.81341
1148.25963
10.61327
53.87063
123.01847
198.02467
275.52947
662.95265
932.80324
1286.67016
10.61202
53.76627
122.46049
196.44174
272.22881
637.31341
880.56008
1153.99624
10.61198
53.76301
122.44309
196.39247
272.12622
636.50134
879.24539
1149.69128
114
(a) (b)
(c) (d)
Figure 5.3 The first four mode shapes of beams; solid-line −: Timoshenko beam,
dash-line −: Euler-Bernoulli beam
115
(a)
(b)
Figure 5.4 Spatial profiles of the forced vibration of the Timoshenko beam:
(a) Displacement distributions, (b) Rotation distributions
116
(a)
(b)
Figure 5.5 The forced vibration of the beams at x = L; (a) Displacement (b) Rotation;
solid-line −: Timoshenko beam, dash-line −: Euler-Bernoulli beam
117
5.5 Spinning Timoshenko Beam
o
x
y
z
f(x,t)
L
Figure 5.6 A spinning beam and its inertial frame oxyz
5.5.1 System Description
Consider a beam in Figure 5.6. The equation of motion based on an inertial frame oxyz
for a spinning, uniform Timoshenko beam of a finite length (L) is given by (Katz ,1988;
Zu and Han, 1992, 1994; Tan and Kuang, 1995)
22
22
22
22
( , ) ( , ) ( , )
( , ), (0, )
( , ) ( , ) ( , ) ( , )
( , ) 0
s
xs
w x t x t w x t
k f x t x L
t x x
x t x t x t w x t
I j J EI k x t
t t x x
(5.31)
where is the spin rate,
z
J is the polar mass moment of inertia and other parameters are
defined in the previous section. The transverse deflections along oy, oz directions and
their corresponding bending angles are represented by, respectively
,
y z y z
w w jw j (5.32)
The boundary and initial conditions of the system are the same form as Eqs. (5.8) and
(5.9), respectively.
118
5.5.2 Dynamic Analysis
In this section, the dynamic response of a spinning Timoshenko beam is obtained by
the same procedures which is used in a simple Timoshenko beam analysis. First, a spatial
state representation is devised by the distributed transfer function formulation. Second,
the s-domain solution is obtained through use of the distributed transfer function. With
the given a spatial state representation and the s-domain solution, the applications are
taken as follows. Eigenvalue problem is defined by eliminating external terms in a spatial
state representation. Frequency response is obtained by the s-domain solution under the
assumption of harmonic excitations. Transient solution is determined by inverse Laplace
transform to the s-domain solution with a residue formula. In this work, the procedures
are also applied to other systems in the same way.
By the distributed transfer function formulation, Eq. (5.31) and the boundary
conditions can be expressed in the same form as Eq. (5.20)
ˆ ˆ ˆ ( , ) ( ) ( , ) ( , ), (0, )
ˆ ˆ (0, ) ( , ) ( )
b b b
x s s x s x s x L
x
s L s s
Fp
MN
(5.33)
where
2
2
0 0 1 0
0 0 0 1
() 0 0 1
00
s
xs
s
s
s
k
I s j J s k
k
EI EI
F (5.34)
The characteristic equation of the state matrix () s F is
42
0 bc
where
119
2
2
2
,
xs
x
ss
I I s i J s k
i J s s
b s c
k EI EI k EI
(5.35)
and the roots of the equation are
22
1,2 3,4
44
,
22
b b c b b c
(5.36)
The case where
2
4 b c b and
2
4 b c b would corresponding to lower and higher
frequency vibration, respectively (Zu and Han, 1992).
Eigenvalue Problem. By vanishing the excitation terms in Eq. (5.33), the eigenvalue
problem of the system can be defined as
( ) ( , ) ( ), (0, )
(0) ( ) 0
bb
x x s x x L
x
L
F
MN
(5.37)
where s is an eigenvalue, and () x the eigenfunction associated with s. By Eq. (5.37)
and
()
( ) (0)
sx
xe
F
, the characteristic equation of the system is
()
det 0
sL
bb
e
F
MN (5.38)
The root , 1,2,...
k
sk , is the kth eigenvalue of the system. The eigenfunctions associated
with
k
s is given by
()
()
k
sx
kk
xe
F
a (5.39)
where
k
a is a nonzero solution of the homogeneous equation ( ) 0
k
s Za .
Frequency and Transient Response. As explained in section 5.4, the s-domain solution of
the system is devised as Eq. (5.22) by the distributed transfer function method. The
120
frequency and transient response can be obtained by Eqs. (5.25) and (5.26), respectively.
The proposed solution method is illustrated on the following example.
Example 5.2: A Spinning Timoshenko Beam
Consider a spinning Timoshenko beam loaded a pointwise impulse force at the point x
f
.
Let the beam be initially at rest. The beam is fixed at the left end, and free at the right end.
A drill is an example for fixed-free boundary conditions. The spatial state representation
of the beam is expressed by
ˆ ˆ ˆ ( , ) ( ) ( , ) ( , ), (0, )
ˆ ˆ (0, ) ( , ) ( )
b b b
x s s x s x s x L
x
s L s s
Fp
MN
(5.40)
where ( ) 0
b
s and
2
0
2
0 0 1 0
0 0 0 1 0
0
ˆ ( ) , ( , ) ( ) ( ) 0 0 1
1
0
00
f
s s
xs
s
I s
s x s t x x
k k
I s j J s k
k
EI EI
Fp (5.41)
1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
,
0 0 0 0 0 0
0 0 0 0 0 0 0
bb
ss
kk
EI
MN (5.42)
For numerical simulation, the non-dimensional parameters of the beam are assigned as
follows:
4
1, 0.02, 10, 10 and 1
s
I EI k L
. Table 5.2 lists the first 8 eigenvalues
of the system which computed by Eq. (5.38), and shows that when a beam is put into a
121
spinning motion (0 200) , its at-rest natural frequency splits into two components:
forward and backward possessions. Although the critical speed was not discussed in this
work, either forward or backward prossession resonance cannot be obtained for a system
with imperfections (i.e., unbalanced mass or initial deflections) or other excitations (Zu
and Han, 1992).
Let the system be subject to an impulse load at the point (x
f
/L=0.5). Using the first 10
modes of the Green’s function, figure 5.8 shows the forced response of a spinning
Timoshenko beam with a spin rate Ω = 50 at time t = 0, 1, 3, 7, 10 (ms). It is seen that the
response of a beam with a backward spin is larger than that of a beam with a forward spin.
In Fig. 5.9, time history of the forced response of a beam for three spin rates (Ω = 50, 100
and 200) at the point (x = L) is illustrated.
122
Table 5.2 The natural frequencies
kk
sj of the spinning Timoshenko beam, 1 j
k ω
k
(Ω=50) ω
k
(Ω=100) ω
k
(Ω=200)
1
2
3
4
5
6
7
8
+15.4072
+77.2266
+156.3130
+236.8079
+316.6539
+394.5874
+470.8338
+545.6969
−7.0748
−38.6491
−99.0029
−166.6836
−238.2643
−310.4627
−382.8991
−455.0778
+20.4850
+104.4248
+196.8064
+282.4252
+364.7925
+444.4238
+521.7185
+597.2216
−4.9175
−29.4898
−80.8878
−141.1684
−207.5906
−276.0148
−345.7799
−415.9498
+30.1007
+160.7161
+290.4563
+391.5616
+477.7182
+558.6432
+636.5041
+712.0732
−2.8729
−20.0983
−57.2962
−104.5371
−160.5166
−220.6503
−284.1069
−349.2773
(+): Forward, (−): Backward
Figure 5.7 Loci of the first 3 eigenvalues of the spinning beam with
a varying spin rate Ω; solid-line –: Forward, dash-line −: Backward
123
(a)
(b)
Figure 5.8 Spatial profiles of the forced vibration of the spinning Timoshenko beam
(Ω = 50), (a) displacement (b) rotation; solid-line –: Forward, dash-line –: Backward
124
Figure 5.9 The time history (0 1) t of displacement (at x = L) of the Timoshenko
beam for a (forward) spin rate Ω; dot-line −: Ω = 50, dash-line −: Ω = 100,
solid-line −: Ω = 200
125
5.6 Timoshenko Beam Column
x
f(x,t)
N
x
N
x
L
Figure 5.10 Schematic of a beam column
5.6.1 System Description
See Figure 5.10. By using the principle of virtual work and Hamilton’s principle for
non-conservative systems, the governing equation of motion for a Timoshenko beam
column applied an axial load N
x
is written as (Kounadis, 1980)
2 2 2
2 2 2
22
22
( , ) ( , ) ( , ) ( , )
( , ), (0, )
( , ) ( , ) ( , )
( , ) 0
sx
s
w x t x t w x t w x t
k N f x t x L
t x x x
x t x t w x t
I EI k x t
t x x
(5.43)
Under the assumption of an axial force in a cross-section is normal to the cross-section,
the motion of a Timoshenko beam column is governed by (Kounadis, 1980; Aristizabal-
Ochoa, 2004)
22
22
22
22
( , ) ( , ) ( , ) ( , )
( , ), (0, )
( , ) ( , ) ( , )
( , ) 0
sx
sx
w x t x t w x t x t
k N f x t x L
t x x x
x t x t w x t
I EI k N x t
t x x
(5.44)
126
In addition, Eq. (5.44) cannot be derived by means of a variational approach.
Arboleda- Aristizabal-Ochoa (2008) derived slope-deflection equations for stability of
Timoshenko beam column structures with semi-rigid connections based on Haringx
model. However, Blaauwendraad (2008) presented that the Haringx theory is inadequate
for shear-weak beam columns.
5.6.2 Dynamic Analysis
The eigensolution, frequency and transient solution of a Timoshenko beam column is
obtained through the same procedures mentioned in the previous section. The spatial state
representation, by the distributed transfer function formulation, is written as
ˆ ˆ ˆ ( , ) ( ) ( , ) ( , ), (0, )
ˆ ˆ (0, ) ( , ) ( )
b b b
x s s x s x s x L
x
s L s s
Fp
MN
(5.45)
where
2
2
0 0 1 0
0 0 0 1
() 00
00
s
s x s x
s
s
k s
s
k N k N
I s k
k
EI EI
F (5.46)
Eigenvalue Problem. The system eigenvalues are obtained by
()
det ( ) det 0
sL
se
F
Z M N
bb
(5.47)
Eq. (5.47) is called by the characteristic equation and has infinite root , 1,2,...
k
sk . The
mode shapes corresponding to
k
s is
127
()
()
k
sx
kk
xe
F
a (5.48)
where
k
a is a nonzero solution of the homogeneous equation ( ) 0
k
s Za .
The s-domain solution is given as
0
ˆ ˆ
ˆ ˆ ( , ) ( , , ) ( , ) ( , ) ( ), , (0, )
L
b
x s x s x s d x s s x L
G p H (5.49)
where
ˆ
G and
ˆ
H are the distributed transfer functions of the system and defined in the
previous section.
By substituting sj in Eq. (5.20),
ˆ
( , , ) xj G and
ˆ
( , ) xj H are the frequency
response functions of the system which can be evaluated in the exact and closed form.
The transient solution of the system is obtained through inverse Laplace transform of the
s-domain solution, Eq. (5.49).
0 0 0
( , ) ( , , ) ( , )d d ( , ) ( ) d
t L t
b
x t x t x t
G p H (5.50)
where the system Green’s functions G and H are the inverse Laplace transforms of the
transfer functions Ĝ and Ĥ, respectively; ( , ) xt p and ()
b
t are the inverse Laplace
transforms of ˆ( , ) xs p and ()
b
s ; and ( , ) xt is the inverse Laplace transform of ˆ ( , ) xs .
Example 5.3: An axially loaded Timoshenko Beam
Here, consider a uniform Timoshenko beam axially loaded by force N
x
. Let the beam
be initially at rest and simply supported at both ends. The boundary conditions and
governing equation are expressed by Eq. (5.45) where
128
1 0 0 0 0 0 0 0
0 0 0 0 0 0 0
,
0 0 0 0 1 0 0 0
0 0 0 0 0 0 0
bb
EI
EI
MN (5.51)
2
2
0 0 1 0
0 0 0 1
() 00
00
s
s x s x
s
s
k s
s
k N k N
I s k
k
EI EI
F (5.52)
For numerical illustration, the geometric and material parameters of the beam are the
length L = 1, the bending stiffness EI = 180, the linear density ρ = 0.84, the rotary inertia
I
ρ
= 7×10
-6
and the shear rigidity k
s
= 6.75×10
6
. Figure 5.11 plots loci of the first 3 natural
frequencies of the beam column, and shows the first critical force (N
cr1
) where the first
natural frequency locus intersects with x-axis. Let the system be subject to an external
load
0
( ) (1 )
t
f t F e
at the point x
f
. Assume that F
0
= 1, σ = 1.0 and x
f
/L = 0.5. Using
the first 10 modes of the Green’s function, the forced vibration of displacement of a beam
column for several axial forces N
x
is computed for the time (0 2) t . Figure 5.12 shows
the time history of a beam column at the location (x/L = 0.8), and Figure 5.13 plots the
spatial profiles of transient response of a beam with a specific axial force (N
x
= 20).
129
Figure 5.11 Loci of the first 3 natural frequencies of the Timoshenko
beam column for a varying axial force
Figure 5.12 The time history (0 2) t of displacement (at x/L = 0.8) of
the Timoshenko beam column for an axial force N
x
; solid-line −: N
x
= 80,
dash-line −: N
x
= 40, dot-line −: N
x
= 20.
130
(a)
(b)
Figure 5.13 Spatial profiles of the forced vibration of the Timoshenko beam column
when N
x
= 20: (a) displacement (b) rotation
131
5.7 Axially Moving and Loaded Timoshenko Beam
5.7.1 System Description
By following the ref. (Lee et al. 2004), the total kinetic and potential energies and the
virtual work due to the external force ( , ) f x t are given by
2 2 2
0
2 2 2
0
1
( ( ) ) ( )
2
1
()
2
L
L
sx
T c w cw I c dx
V EI k w N w dx
(5.53a)
0
( , )
L
W f x t wdx
(5.53b)
The extended Hamilton’s principle is
2
1
( ) 0
t
t
T V W dt
(5.54)
By substituting Eqs. (5.10) and (5.12) to the extended Hamilton’s principle, the equation
of motion for the uniform Timoshenko beam can be derived as
2
2
2 ( , ), (0, )
20
sx
s
c w cw w k w N w f x t x L
I c c EI k w
(5.55)
where c is the axial velocity, the dot (∙) denotes the derivative with respect to time and the
prime (ʹ) the derivative with respect to spatial coordinate x. The shear force Q(x) and the
bending moment M(x) are given as
( ) ( )
()
sx
Q x k w c w cw N w
M x EI I c c
(5.56)
If c is equal to zero ( 0 c ), Eq. (5.55) becomes Eq. (5.43).
132
5.7.2 Dynamic Analysis
The governing equation, boundary and initial conditions are written in the spatial state
form
ˆ ˆ ˆ ( , ) ( ) ( , ) ( , ), (0, )
ˆ ˆ (0, ) ( , ) ( )
b b b
x s s x s x s x L
x
s L s s
Fp
MN
(5.57)
where
22
12
,
sx
k c N EI I c
and
2
1 1 1
2
2 2 2
0 0 1 0
0 0 0 1
2
0
()
2
0
s
s
s
k s cs
s
I s k I cs
k
F (5.58)
With the distributed transfer functions
ˆ
G and
ˆ
H , the s-domain solution is given as
0
ˆ ˆ
ˆ ˆ ( , ) ( , , ) ( , ) ( , ) ( ), , (0, )
L
b
x s x s x s d x s s x L
G p H (5.59)
Let all external disturbances be zero in Eq. (5.57). The characteristic equation of the
system can be derived as det ( ) det( ( ,0, )) 0
bb
s L s Z M N Φ , from which the system
eigenvalues s
k
, k = 1,2,…, are determined. The corresponding eigenfunctions are obtained
as
()
()
k
sx
kk
xe
F
a (5.60)
where
k
a is a non-trivial solution of the equation ( ) 0
k
s Za .
133
The transient solution of the system is obtained through inverse Laplace transform of
the s-domain solution, i.e.
1
ˆ ( , ) ( , ) x t x s
. By Eq. (5.10) and replacing sj , the
frequency response functions
ˆ
( , , ) xj G and
ˆ
( , ) xj H are evaluated in the exact and
closed form.
Example 5.4: An axially moving Timoshenko Beam
Let a beam be initially at rest. The beam is subject to the external force,
1.0
( , ) (1 ) ( 0.5)
t
f x t e x
. The beam is simply supported at both ends. The spatial
state representation of the beam is expressed by Eqs. (5.57) and (5.58), and boundary
matrices are assigned as
2
2
0 0 0 0 1 0 0 0
0 0 0 0 00
,
1 0 0 0 0 0 0 0
00 0 0 0 0
bb
I cs I c EI
I cs I c EI
MN (5.61)
The boundary matrices are derived from the boundary conditions that the displacement
and moment are zero at both ends. The analysis of the beam contains the frequency-
dependent boundary conditions and the non-self adjoint effect. For numerical simulation,
the non-dimensional parameters of the beam are assigned as follows: ρ = 1, I
ρ
= 0.02, EI
= 10, k
s
= 10
4
, N
x
= 20 and L = 1. Table 5.3 lists the eigenvalues of a beam for various
moving speeds (c = 2, 4 and 6) which are computed by det ( ) 0 s Z . The moving speeds
are below the critical speed, and the eigenvalues have imaginary values. The real part of
eigenvalues which makes the system unstable are observed within certain ranges over the
134
critical speed. (Lee et al., 2004) Figure 5.14 shows the time history of displacement of a
beam at the location (x/L = 0.8) for 02 t .
Table 5.3 The eigenvalues
kk
sj of the moving Timoshenko beam 1 () j
k ω
k
(c=2) ω
k
(c=4) ω
k
(c=6)
1
2
3
4
5
6
7
8
24.4676
89.7286
164.8930
240.9584
316.0267
389.9938
463.0712
535.4738
21.4085
86.3023
159.9017
234.2165
307.5319
379.7622
451.1273
521.8282
15.8182
80.5231
151.6003
223.0073
293.4044
431.2431
499.1103
566.4860
Figure 5.14 The time history (0 2) t of displacement (at x/L = 0.8) of
the axially moving Timoshenko beam with a speed, c;
solid-line −: c = 6, dash-line −: c = 4, dot-line −: c = 2
135
Chapter 6. Dynamic Analysis of Stepped
Bars, Shafts and Strings
6.1 Introduction
This chapter is concerned with transient vibration analysis of stepped second-order
distributed dynamic systems. The stepped system in consideration is an assemblage of
one-dimensional distributed components (like elastic bars, torsional shafts and taut
strings), and lumped masses. Recently, Yang (2010) developed the analytical method for
transient vibration analysis of stepped systems composed of uniformly distributed
components. This study is on the extent of his investigation (Yang, 2010). The keys in the
transient solution method are to determine an exact state transition matrix and system
eigenvalues. In this study, the analytical method is developed to obtain state transition
matrix and eigenvalues which are required to determine transient solutions of stepped
systems. The proposed method for transient vibration analysis does not depend on system
eigenfunctions, is able to accurately predict jumps in stress and strain distributions. The
method is numerically efficient because its utility only involves simple operations of two-
by-two matrices.
6.2 System Description
The stepped and non-uniformly distributed system in consideration is an assembly of N
one-dimensional, serially non-uniform components are interconnected; see Figure 6.1,
136
where ,
i
x 1,2, , 1 iN , are the interior nodes at which adjacent components are inter-
connected,
0
and
N
xx are the boundary nodes representing the ends of the system, and
1 i i i
l x x
is the length of the ith component. Set
0
0 x , so that
12 ii
x l l l for
1,2, , in . The displacement ( , )
i
w x t of the ith component is governed by the wave
equation
2
1 2
( , ) ( , )
( ) ( ) ( , ), ( , )
ii
i i i i i
w x t w x t
x x f x t x x x
t x x
(6.1)
where ()
i
x is an inertia parameter, ()
i
x a stiffness parameter, and ( , )
i
f x t an external
load. At the nodes of the stepped system there may be spring constraints and mounted
lumped masses.
The interface between the ith and (i+1)th components at node
i
x is described by the
matching conditions
1
2
1
1 2
( , ) ( , )
( , ) ( , ) ( , )
( ) ( , ) ( ) ( ) ( )
i i i i
i i i i i i
i i i i i i i i i i i
w x t w x t
w x t w x t w x t
m x k w x t x q t k z t
t x x
(6.2)
for 1,2, , 1 iN , where
i
m is a lumped mass or a rigid disk,
i
k the coefficient of a
spring; ()
i
qt an external force applied at the lumped mass, and ()
i
zt a foundation motion
of the node. Also, the system is subject to the boundary conditions
2
0 1 0 2
2
10 2
at : ( )
at : ( )
LL
n R R
ww
x x m a a w t
tx
ww
x x m b b w t
tx
(6.3)
137
where
L
m and
R
m are end masses;
0
a ,
1
a ,
0
b and
1
b are constants that are properly
assigned to characterize different types of boundaries; and ()
L
t and ()
R
t are related to
prescribed boundary excitations.
The initial conditions of the system are assigned as
0, 0, 1
( ,0)
( ,0) ( ), ( ), ,
i
i i i i i
wx
w x u x v x x x x
t
(6.4)
for 1,2, , iN , where
0,
()
i
ux and
0,
()
i
vx are given profiles of initial displacement and
velocity of the ith component.
Figure 6.1 Schematic of an n-component stepped distributed dynamic system
6.3 Distributed Transfer Function Formulation
Laplace transform of Eqs. (6.1) to (6.4) with respect to time gives
2
,1
( , )
( ) ( , ) ( ) ( , ), ( , )
i
i i i eI i i i
w x s
x s w x s x p x s x x x
xx
(6.5)
2 1
1 0 1
2
10
(0, )
( ) (0, ) ( )
( , )
( ) ( , ) ( )
LL
n
R n R
ws
a m s a w s s
x
w L s
b m s b w L s s
x
(6.6)
138
where the over-bar stands for Laplace transform, s is the Laplace transform parameter,
and
, 0, 0,
0,1 0 0,1 0
0, 0,
( , ) ( , ) ( ) ( ) ( )
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
eI i i i i i
L L L
R R R n n n n
p x s f x s x su x v x
s s m su x v x
s s m su x v x
(6.7)
Take Laplace transform of Eq. (6.2) with respect to time,
1
( , ) ( , )
i i i i
w x s w x s
2 1
1 0,
( , ) ( , )
ˆ ( ) ( ) ( ) ( , ) ( )
i i i i
i i i i i i i i e i
w x s w x s
x x m s k w x s q s
xx
(6.8)
where
0, 0, 0,
ˆ ( ) ( ) ( ) ( ( ) ( ))
e i i i i i i i i i
q s q s k z s m su x v x
Equation (6.1) is rewritten in a first-order state form
1
ˆ ˆ ( , ) ( , ) ( , ) ( , ), ( , )
i i i i i i
x s x s x s x s x x x
x
Fp (6.9)
where
,
2
01
( , )
( , ) 0
ˆ ( , ) , ( , ) , ( , )
( , ) ( ) ( )
1 ()
( ) ( )
i
eI i
i i i
i ii
i
ii
w x s
p x s
x s x s x s
w x s x s x
x
x xx
η Fp (6.10)
where ( ) ( ) / x d x dx .
For convenience of analysis, define a global domain
0 1 1 2 1
( , ) ( , ) ( , )
NN
x x x x x x
U U U
and convert Eq. (6.9) to a global state equation
ˆ ˆ ( , ) ( , ) ( , ) ( , ), x s x s x s x s x
x
Fp (6.11)
where the global quantities
139
ˆ ˆ ( , ) ( , )
i
x s x s , ( , ) ( , )
i
x s x s FF , ( , ) ( , )
i
x s x s pp
The boundary and matching conditions are also cast into a global state form
0
ˆ ˆ ( , ) ( , ) ( )
ˆ ˆ ( , ) ( , ) ( ), 1, 2, , 1
b b N b
i i i i
x s x s s
x s x s s i N
MN
T ν
(6.12)
where
2
01
2
01
00 ()
, , ( )
()
00
L
L
b b b
R R
s
m s a a
s
m s b b s
MN γ (6.13)
0,
2
1
11
10
0
ˆ ()
, ( )
()
()
1
( ) ( )
ei
ii
i i i i
ii
i i i i
qs
s
m s k x
x
xx
T ν
In this study, the state Eq. (6.11) is solved through use of the state transition matrix
which is the unique solution of (Chen, 1999).
( , , ) ( , ) ( , , ), , x s x s x s x
x
Φ F Φ (6.14)
Following Ref. (Yang, 2010) and the properties of state transition matrix, it can be
shown that the state transition matrix for the stepped system is
1 0 1
1 1 1 1 2 1 1 1 1
( , , ) , ,
( , , )
( , , ) ( , , ) ( , , ), , , 2
i i i i i i i i
x s x x x
xs
x x s x x s x s x x x i N
TT
(6.15)
where ( , , )
ii
xs is the state transition matrix of the ith component.
Note that the state transition matrices of non-uniform components are represented in
chapter 3.5. With the state transition matrix given in Eq. (6.15), the s-domain response of
the stepped system is taken as
140
1
0
1
ˆ ˆ ˆ
ˆ ˆ ( , ) ( , , ) ( , ) ( , ) ( ) ( , , ) ( )
N
L
b i i
i
x s x s x s d x s s x x s s
G p H G ν (6.16)
where
0
ˆ
( , ) ( , , ),
ˆ
( , , )
ˆ
( , ) ( , , ),
b
bN
x s x s x
xs
x s x s x
HM
G
HN
(6.17a)
1
0
ˆ
( , ) ( , , ) ( ) x s x x s s
HZ (6.17b)
with
0
( ) ( , , )
b b N
s x x s Z M N (6.18)
The matrices
ˆ
G and
ˆ
H are called the distributed transfer functions of the stepped system,
and () s Z the boundary impedance matrix.
As shall be seen in the subsequent sections, the eigensolutions and transient response
of the stepped system are determined by the spatial state equations and s-domain solution.
6.4 Dynamic Analysis
6.4.1 Eigensolutions
Define the eigenvalue problem of the stepped system by vanishing the excitation terms
in Eqs. (6.11) and (6.12),
( ) ( , ) ( ), x x s x x
x
ψ F ψ (6.19)
subject to the conditions
( ) ( ), 1,2, , 1
i i i
x x i N ψ T ψ (6.20a)
141
0
( ) ( ) 0
b b N
xx M ψ N ψ (6.20b)
where s is an eigenvalue, and () x ψ
the eigenfunction associated with s. Write
0
( ) ( , , ) x x x s ψ Φ a (6.21)
with a being a vector to be determined, which satisfies Eqs. (6.19) and (6.20a). Plugging
Eq. (6.21) into the boundary condition (6.20b) gives
0
( ( , , )) 0
b b N
x x s MN Φ a (6.22)
The characteristic equation of the stepped system is then
0
det ( ) det( ( , , )) 0
b b N
s x x s Z M N Φ (6.23)
The roots of Eq. (6.23) can be written as , 1, 1,2, ,
kk
s j j k where
k
is the
kth natural frequency of the system. From Eqs. (6.23) and (6.15), the natural frequencies
are the non-negative roots of the real-valued transcendental equation.
( ) det( ( )) 0
bb
M N W (6.24)
where
1 1 1 1 2 1 1 1 0
( ) ( , , ) ( , , ) ( , , )
N N N N N N N
x x j x x j x x j
W Φ T Φ T Φ (6.25)
Although eigenfunctions are not needed in the proposed transient analysis, they are
useful in free vibration analysis. The eigenfunction (mode shape) associated with
k
j is
given by
0
( ) ( , , )
k k k
x x x j ψ Φ a (6.26)
where
k
a is a non-zero solution of the homogeneous equation ( ) 0
k
j Z .
142
6.4.2 Computation of eigenvalues
Eigenvalues of the system are determined as the non-negative roots of the real-valued
transcendental equation, Eq. (6.24) and infinite. Geometrically, a root of Eq. (6.24) is the
value of ω at which the graph of the equation intersects the ω-axis (See Fig. 6.2). Missing
some eigenvalues may cause significantly different result in transient vibration analysis.
The possibility of missing any is increased when eigenvalues form irregular distribution
of points in frequency plane. Stepped distributed systems and multilayered smart
structures may have close grouping of eigenvalues which causes increasing the risk of
missing some eigenvalues (Lee et al., 2002). As an example, see Figure 6.2 (a); the
characteristic equation of the three-segment shaft in Reference (Yang, 2010). In such
cases, very fine intervals of ω are required to avoid the possibility of missing any root of
transcendental function because the function is very badly behaved for some intervals.
Nevertheless, the standard root-finding methods always carry the risk of missing some of
the eigenvalues in the course of computation (Wu, 1997; Lee et al., 2002), and the
computation cost is increased drastically. This section is concerned with determination of
eigenvalues (natural frequencies) of the stepped system.
Without loss of generality, it is assumed that eigenvalues of the stepped system are
distinct (no repeated eigenvalues) since the rank of ()
k
j Z in Eq. (6.23) is always equal
to 1. The roots of function () form a more or less irregular distribution of the points in
the plane. Let us pass to the derivative ( ) / dd
and consider its roots; they are in
general different from those of () and form a new distribution of points. There is a
sequence of distributions of points formed by the roots of the successive derivatives.
143
Depending on the analytic character of the function, a recognizable trend for functions of
certain sorts is observed (Polya, 1943). In the stepped systems, the roots of ( ) / dd
have almost regular distribution although those of ()
have irregular distribution in Fig.
6.2. If the roots of function form a regular distribution, the roots can be determined by
standard root-finding techniques without any difficulty.
With the assumption that the eigenvalues of the system are distinct, the bounds of
natural frequencies of the stepped system can be determined with satisfying the condition
1
( ) ( ) 0
kk
where
k
is the non-negative roots of ( ) / dd and only one
eigenvalue located within the bound
1
( , )
kk
.
1
, 1,2,
kkk
k
(6.27)
When a root of the characteristic equation is equal to the root of its
derivative )
kk
, the system has repeated eigenvalues at the point. The proposed
formula can be extended to a system which has repeated eigenvalues. However it is not
considered in this study. In a summary, the determination of bounds of natural
frequencies takes two steps:
Step1: compute roots of ( ) / dd .
Step 2: determine the bounds
k
with satisfying
1
( ) ( ) 0
kk
With the intervals (bounds) described in this section guarantee to contain a root - if
k
and
1 k
are the endpoints of the interval then ()
k
must differ in sign from
1
()
k
,
the standard root-bracketing algorithms are used in this study. This ensures that the
function crosses zero at least once in the interval. In addition, the risks of intensive
144
computation rate and missing any can be alleviated to determine eigenvalues of the
system. With proper scaling, such as
1
1
( ) / (1 ),
n
n constant, the functions ()
and ( ) / dd are bounded in the entire region [0, ) (Yang, 2010). Hence, standard
root-searching techniques are applicable to Eq. (6.24) for accurate solution of the
eigenvalues in the entire region. In the proposed transient analysis, the derivative of
characteristic function ( ) / d s ds Z is also required to determine transfer function residues
as shall be seen in the subsequent section.
(a) (b)
Figure 6.2 (a) Characteristic function of the stepped shaft from ref. (Yang, 2010),
(b) its derivative of the characteristic function with respect to ω
6.4.3 Transient response
In this research, exact transient solutions are obtained by the distributed transfer
function formulation obtained in section 6.3 and by a residue formula for inverse Laplace
transform. Inverse Laplace transform of Eq. (6.16) yields Green’s function formula
145
1
1
0
1
1
00
1
0, 0,
0
1
( , ) ( , , ) ( , )
1
() 0
1
( , , )( ( ) ( )) ( , )
() 1
0
( , , ) ( ) ( , , ) ( )
1
i
i
i
n
tx
i ST
x
i
i
n
tt
L
i i i i ST
i R i
x
i
ii ST
x
i
x t x t f d d
t
x x q k z d x t d
t
x t u x t v d
t
η G
GH
GG
1
1
0, 0,
1
1
0,1 0 0,1 0
0, 0,
0
( , , ) ( ) ( , , ) ( )
1
( ) ( )
( , ) ( , )
( ) ( )
i
n
i
n
i
i i i i i i ST
i
i
LL
R n n R n n
m
x x t u x x x t v x
t
m u x m v x
x t x t
m u x m v x t
GG
HH
(6.28)
for , x where ( , ) xt is the inverse Laplace transform of ˆ ( , ) xs ; and the matrix
Green’s functions G and H are the inverse Laplace transforms of the transfer functions
ˆ
G and
ˆ
H , respectively. The first four terms on the right-hand side of Eq. (6.28)
represent the effects of external loads, boundary excitations and initial disturbances, and
determination of the transient response of a stepped system are the contributions of the
lumped masses due to the initial disturbances.
Following the Ref. (Yang, 2010), the Green’s functions of the system are given by
1
1
( , , ) 2 Re ( ,0, ) (0, , )e
( , ) 2 Re ( ,0, ) e
k
k
st
k k b k
k
st
kk
k
x t x s s
x t x s
G R M
HR
(6.29)
The transient response of the stepped distributed system can be determined by Eq.
(6.10) if the transfer function residues are known. The transfer function residue can be
expressed as
146
adj ( )
()
tr adj ( )
k
k
k
k
sj
j
ds
j
ds
Z
R
Z
Z
(6.30)
The above transient solution is of the exact form; no discretization or approximation has
been made. This transient solution method is numerically efficient because it only
involves operations of two-by-two matrices, regardless of the number of components.
Furthermore, the method treats different system configurations (number of components,
physical parameters, constraints, lumped masses and boundary conditions) in a symbolic
manner, and avoids tedious derivations and messy expressions that are encountered in
many analytical methods.
6.5 Example
The proposed transient analysis is demonstrated on a three-segment bar in longitudinal
vibration.
Example 6.1: A Three-Segment Bar in Longitudinal Vibration
A stepped bar is composed of three non-uniform and uniform components; see Fig. 6.3.
The longitudinal displacement ( , )
i
u x t of the ith components is governed by
2
1
2
( , ) ( , )
( ) ( ) ( ) ( ), ( , )
ii
i i f i i
u x t u x t
x EA x f t x x x x x
xx t
(6.31)
The system is subject to the boundary conditions
0
( , ) 0 u x t ,
3
3 3 3
( , ) 1
( ) ( , ) ( )
b
u x t
EA x u x t z t
kx
(6.32)
147
The matching conditions at internodes x
i
, (i=1, 2) are
1
( , ) ( , )
i i i i
u x t u x t
,
1
1
( , ) ( , )
( ) ( )
i i i i
i i i i
u x t u x t
EA x EA x
xx
(6.33)
Assume that the linear density and the longitudinal rigidity of components have the
exponential, uniform, and conical distributions.
2
2
1 0,1 2 0,2 3 0,3 3
2
2
1 0,1 2 0,2 3 0,3 3
( ) , ( ) , ( ) 1 /
( ) , ( ) , ( ) 1 /
ax
ax
x e x x x l
EA x EA e EA x EA EA x EA x l
(6.34)
where
0,i
and
0,i
E (i=1,2,3) are constant.
The state transition matrix can be computed by a fundamental matrix. By Eqs. (3.19),
(3.25), (3.23) and (3.26), a fundamental matrix of the 1
st
, 2
nd
, and 3
rd
component can be
written as
1
1
cosh( ) sinh( )
( , )
cosh( ) sinh( ) sinh( ) cosh( )
ax
xx
x s e
a
a x x x x
U (6.35a)
2
()
2
1
cosh( ) sinh( )
( , )
sinh( ) cosh( )
sx
xx
x s e
xx
F
U (6.35b)
33
3
22
3 3 3 33
cosh( ) 1 sinh( )
1 / 1 /
( , )
sinh( ) cosh( ) cosh( ) 1 sinh( )
1 / 1 / (1 / ) (1 / )
xx
x l x l
xs
x x x x
x l L x l l x l x l
U (6.35c)
where
0, 0,
/
ii
s
EA
, and
22
a .
148
By the matching conditions (6.2) at internodes, the interface matrix T
i
in Eq. (6.12) is
given by
1
10
, 1,2 ()
0
()
i ii
ii
i EA x
EA x
T (6.36)
For numerical simulation, the non-dimensional values of the system parameters are
assigned as follows
0,1 0,1 1
0,2 0,2 2
0,3 0,3 3
Segment 1: 2.0, 600, 1.0
Segment 2 : 1.5, 300, 0.7
Segment 3: 1.2, 250, 1.3
0.5, 0.2, 100
EA l
EA l
EA l
ak
(6.37)
The natural frequencies of the stepped bar in the first 20 modes and several higher
modes are computed by Eq. (6.24). The natural frequencies by the proposed method are
compared with the results by a finite element method (FEM), and are listed in Table 6.1.
The FEM with 150 and 300 elements gives good results on the first 20 natural
frequencies; however the errors in the FEM predictions grow when higher-mode natural
frequencies are computed. By Eq. (6.26), the first five mode shapes and the 14th mode
shape are plotted in Figure 6.4.
Let the bar be initially at rest. In Fig. 6.3, a longitudinally vibrating elastic bar, which
is constrained by the spring at the right end, is subject to a foundation motion
0
( ) sin
b
z t Z t . Choose
0
Z = 0.2, =215. Through use of the state transition matrix,
the transient response of the stepped bar is computed by using the first 40 modes in Eq.
(6.28) and (6.29). Figure 6.5 (a)-(d) show the traveling wave characteristics of the
149
boundary disturbance at earlier times. Once the traveling wave reach the left boundary,
the bar vibration becomes a mixture of incoming and outgoing waves, as seen in Fig.
6.5(e). After some time, the vibration of the system eventually settles in a pattern shown
in Fig. 6.5(f). This pattern is close to the 14
th
mode shape of the stepped bar, because the
boundary excitation frequency is near the 14
th
natural frequency of the system: see
Table 6.1.
150
u(x,t)
Z
b
(t)
k
x
x
0
x
2
x
1
x
3
l
1
l
2
l
3
Figure 6.3 A three-segment bar in longitudinal vibration
151
Table 6.1 The natural frequencies
k
of the stepped bar (in rad/s)
k
DTFM
Finite Element Method
30 elements 150 elements 300 elements
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
25
30
50
100
200
9.97433
24.82032
42.68476
55.66625
71.38696
88.67533
104.62378
118.76719
135.29514
153.12845
166.31815
183.06114
199.72937
215.91168
229.49608
247.59928
263.82140
277.73483
294.36693
311.65082
389.43511
470.71774
787.81519
1583.72050
3176.01275
9.97865
24.84742
42.80590
55.97478
72.05003
89.98216
106.60045
121.89577
139.94601
159.66931
174.57008
194.67447
214.54536
233.54349
252.13929
275.43968
296.14684
315.31900
340.01564
361.95807
465.57926
580.91559
9.97450
24.82141
42.68961
55.67856
71.41346
88.72752
104.70273
118.89156
135.48063
153.39033
166.64490
183.52410
200.32533
216.61894
230.39898
248.74270
265.16004
279.26398
296.31607
313.88364
393.73615
478.63255
825.31224
1852.75141
9.97437
24.82059
42.68597
55.66933
71.39359
88.68838
104.64352
118.79828
135.34151
153.19391
166.39980
183.17682
199.87831
216.08843
229.72157
247.88498
264.15596
278.11652
294.85361
312.20889
390.50686
472.69533
797.18050
1658.00015
3715.75581
152
(a) Mode 1 (b) Mode 2
(c) Mode 3 (d) Mode 4
(e) Mode 5 (f) Mode 14
Figure 6.4 Mode shapes of the stepped bar: (a) mode 1; (b) mode 2;
(c) mode 3; (d) mode 4; (e) mode 5; (f) mode 14.
153
(a) t = 0.05s
(b) t = 0.15s
(c) t = 0.2s
(d) t = 0.3s
(e) t = 1.5s
(f) t = 3.0s
Figure 6.5 Spatial profiles of the transient displacement of the bar subject to
a sinusoidal boundary excitation
0
( ) sin
b
z t Z t , with
0
Z = 0.005, = 215.
154
Chapter 7. Multi-Body Structures
7.1 Introduction
The multi-body structure in consideration is composed of multiple continua that are
interconnected at certain points and/or in certain regions. Each continuum may be under
constraints, combined with lumped parameter systems, and subject to external and
boundary disturbances. The generalized formulation of multi-body distributed parameter
systems is represented in Chapter 2. In the analysis, the multi-body distributed parameter
system is first divided into a number of subsystems; the distributed transfer functions of
the system are determined in exact and closed form by synthesis of a spatial state form. In
this chapter, the proposed method is applied to eigenvalue problem of various multi-body
distributed parameter systems. The method deals with the system in a unified and
systematic manner, while it is convenient for computer coding. System eigenvalues
contain information required in modeling and understanding the systems. Eigenvalues are
often plotted versus a system parameter by a family of “root” loci. Eigenvalue loci
approach closely, and rapidly veer away. This case, called “curve veering”, has been
usually observed in approximate solutions associated with discretized models. Even if
exact solutions of a system are available, the curve veering can be observed in a certain
structure (Perkins and Mote, 1986). The understanding of such eigenvalue loci is
important for a design of system.
155
7.2 Stepped Beams
7.2.1 System Description
The stepped distributed system in consideration is an assembly of N one-dimensional,
Euler-Bernoulli beams are interconnected serially, where ,
i
x i =1,2,…,N-1 are the interior
nodes at which adjacent components are inter-connected,
0
x and
N
x are the boundary
nodes representing the ends of the system, and
1 i i i
l x x
is the length of the ith
component. The displacement ( , )
i
w x t of the ith component is governed by (Yang, 2005)
24
1 24
( , ) ( , )
( , ) ( , ), ( , )
ii
i i i i i i i
w x t w x t
EI w x t f x t x x x
tx
(7.1)
where
i
is an inertia parameter,
i
EI a stiffness parameter,
i
a coefficient of the elastic
foundation and ( , )
i
f x t an external load.
Also, the system is subject to the boundary conditions
23
0 0 1 ,1 23
2
2 3 ,2 2
(0, ) (0, )
at : (0, ) ( )
(0, ) (0, )
()
LL
L
w t w t
x x m a w t a t
tx
w t w t
a a t
xx
(7.2a)
23
0 1 ,1 23
2
2 3 ,2 2
( , ) ( , )
at : ( , ) ( )
( , ) ( , )
()
n R R
R
w L t w L t
x x m b w L t b t
tx
w L t w L t
b b t
xx
(7.2b)
where
i
a and ( 0,1,2,3)
i
bi are constants that are properly assigned to characterize
different types of boundaries;
L
m and
R
m are end masses, and
,1 ,2 ,1
( ), ( ), ( )
L L R
t t t , and
156
,2
()
R
t
are related to prescribed boundary excitations. The initial conditions of the
system are assigned as
0, 0, 1
( ,0)
( ,0) ( ), ( ), ,
i
i i i i i
wx
w x u x v x x x x
t
(7.3)
for 1,2, , iN , where
0,
()
i
ux and
0,
()
i
vx are given profiles of initial displacement and
velocity of the ith component.
Assume that a constraint is imposed at the interior nodes at which adjacent components
are inter-connected. Figure 7.1 shows six types of constraints considered, where k
r
is the
coefficient of a torsional spring, k
t
is the coefficient of a translational spring, m and I are
the inertia and moment of inertia of a lumped mass and I
D
is the moment of inertia of a
hinged rigid disk. The interface with the constraint at node x
i
is described by the
matching conditions as follows.
(C1) Hinge:
1
22
11
1 22
( , ) 0, ( , ) 0
( , ) ( , ) ( , ) ( , )
,
i i i i
i i i i i i i i
ii
w x t w x t
w x t w x t w x t w x t
EI EI
x x x x
(7.4a)
(C2) Sliding constraint:
1
1
33
1
1 33
( , ) ( , )
( , ) ( , ), 0, 0
( , ) ( , )
i i i i
i i i i
i i i i
ii
w x t w x t
w x t w x t
xx
w x t w x t
EI EI
xx
(7.4b)
157
(C3) Spring constraint:
1
1
22
1
122
33
1
133
( , ) ( , )
( , ) ( , ),
( , ) ( , ) ( , )
( , ) ( , )
( , )
i i i i
i i i i
i i i i i i
i i r
i i i i
i i t i i
w x t w x t
w x t w x t
xx
w x t w x t w x t
EI EI k
x x x
w x t w x t
EI EI k w x t
xx
(7.4c)
(C4) Hinged disk:
1
1
2 2 3
1
1 2 2 2
( , ) ( , )
( , ) 0, ( , ) 0,
( , ) ( , ) ( , ) ( , )
i i i i
i i i i
i i i i i i i i
i i r D
w x t w x t
w x t w x t
xx
w x t w x t w x t w x t
EI EI k I
x x x t x
(7.4d)
(C5) Sliding mass:
1
1
3 3 2
1
1 3 3 2
( , ) ( , )
( , ) ( , ), 0, 0
( , ) ( , ) ( , )
( , )
i i i i
i i i i
i i i i i i
i i t i i
w x t w x t
w x t w x t
xx
w x t w x t w x t
EI EI k w x t m
x x t
(7.4e)
(C6) Mass and springs:
1
1
2 2 3
1
1 2 2 2
3 3 2
1
1 3 3 2
( , ) ( , )
( , ) ( , ),
( , ) ( , ) ( , ) ( , )
( , ) ( , ) ( , )
( , )
i i i i
i i i i
i i i i i i i i
i i r D
i i i i i i
i i t i i
w x t w x t
w x t w x t
xx
w x t w x t w x t w x t
EI EI k I
x x x t x
w x t w x t w x t
EI EI k w x t m
x x t
(7.4f)
for 1,2, , 1 iN . Note that Eq. (8.4c) with vanishing spring coefficients describes the
matching condition at interior node (x
i
) without a constraint.
158
7.2.2 Distributed Transfer Function Formulation
To this end, take Laplace transform of Eq. (7.1) with respect to time, and cast the
resulting equations into a spatial state form
1
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ), ( , )
i i i i i i
x s x s x s x s x x x
x
Fp (7.5a)
10
ˆ ˆ ( , ) ( , ) ( )
b b N N b
x s x s s MN (7.5b)
The matching conditions, Eq. (7.4), are Laplace transformed and cast into a state form
(Yang, 2010)
1
ˆ ˆ ( , ) ( , ) 0, 1,2, , 1
i i i i i i
x s x s i N
AB (7.5c)
where
i
A and
i
B are interface matrices describing the matching conditions at the ith
node, s is the Laplace transform parameter, and
23
23
( , ) ( , ) ( , )
ˆ ( , ) ( , )
T
i i i
ii
w x s w x s w x s
x s w x s
x x x
,1
,2 ,
,1
2
,2
0 1 0 0
0
0 0 1 0
0
0 0 0 1 ( , ) , ( , ) ,
0
()
000 1
L
L eI i
i i b
R i
ii
R
i
p
x s x s
EI
s
EI
Fp (7.6)
2
01
23
2
01
23
0 0 0 0 00
0 0 0 0 00
,
00 0 0 0 0
00 0 0 0 0
R
bb
L
m s a a
aa
m s b b
bb
MN
with
159
, 0. 0,
,1 ,1 0,1 0 0,1 0 ,2 ,2
ˆ
( , ) ( , ) ( ) ( )
ˆ ˆ ( ) ( ) ( ) ( ) , ( ) ( )
eI i i i i i
L L L L L
p x s f x s su x v x
s s m su x v x s s
(7.7)
,1 ,1 0, 0, ,2 ,2
ˆ ˆ ( ) ( ) ( ) ( ) , ( ) ( )
R R R N N N N R R
s s m su x v x s s
Thus, the original boundary-initial value problem is converted to the s-domain state
equation (7.5a) subject to the boundary (7.5b) and matching conditions (7.5c).
By defining a global domain
0 1 1 2 1
( , ) ( , ) ( , )
NN
x x x x x x
U U U for convenience of
analysis and the new state vector,
12
ˆ ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
T
T T T
N
x s x s x s x s (7.8a)
Eq. (7.5a) is rewritten in a global domain
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ), x s x s x s x s x
x
Fp (7.8b)
where
12
12
ˆ ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
( , ) ( , ) ( , ) ( , )
T
T T T
N
N
x s x s x s x s
x s diag x s x s x s
p p p p
F F F F
(7.8c)
and
12 N
diag a a a is a diagonal matrix whose diagonal entries
12 N
a a a .
With a global state vector, the boundary and matching conditions, Eq. (7.5b) and (7.5c),
can be expressed in a single state equation
ˆ ˆ ( , ) ( , ) ( )
G s G e G
x s x s s MN (7.9a)
where
160
1 0 2 1 1
1 1 2 2
ˆ ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
ˆ ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
T
T T T
s N N
T
T T T
e N N
T
Gb
x s x s x s x s
x s x s x s x s
00
(7.9b)
11
11
,
bb
GG
NN
M 0 0 0 0 N
0 A B 0 0
MN
00
0 0 A 0 0 B 0
with 0 is a zero matrix.
Note that the state equation of a stepped beam, Eq. (7.8b) and (7.9a), has the same
form as that of a single-body beam.
7.2.3 Applications
By vanishing the excitation terms in Eq. (7.8) and (7.9), define the eigenvalue problem
of the system
()
( , ) ( ),
x
x s x x
x
F
(7.10)
subject to the conditions
( ) ( ) 0
G s G e
xx MN (7.11)
where and
se
xx are representing the end nodes of the subsystems, s is an eigenvalue, and
() x the eigenfunction associated with s. Assume that
( ) ( , , )
s
x x x s a (7.12)
with a being a constant vector to be determined and ( , , ) xs the state transition matrix
which is the unique solution of (Chen, 1999)
161
( , , ) ( , ) ( , , ), , x s x s x s x
x
F (7.13)
The () x , which satisfies Eq. (7.12), is plugged in the boundary condition (7.11), yields
( ) ( , , ) 0
G G G e s
s x x s Z a M N a = (7.14)
The characteristic equation of the system then is det ( ) 0
G
s Z , from which the system
eigenvalues , 1,2,
k
sk , are determined. The corresponding eigenfunction is obtained
by ( ) ( ,0, )
k k k
x x s a where
k
a is a nonzero solution of the homogeneous equation
( ) 0
Gk
s Za . Without loss of generality, assume that eigenvalues of a stepped beam are
distinct for 0 s . Let the matrix ()
G
s Z be a 4N-by-4N square matrix. The rank of matrix
()
Gk
s Z is equal to 4N-1 because of det ( ) 0
Gk
s Z and adj ( ) 0
Gk
s Z . The algorithm
in Section 6.4.2 which reduces computation rate and the risk of missing eigenvalues is
applicable for determination of eigenvalues of a stepped beam.
As shown in the previous section, the state equations of a stepped beam have the same
forms as those of a single-body beam. The frequency and transient response of the system
can be obtained through the s-domain solution of the system. The proposed method is
illustrated on the following example.
Example 7.1: A three-segment beam with hinged supports
Consider a three-segment beam with hinged supports at the interior nodes; see Fig. 7.2.
The beam is fixed at the left end (point x
0
) and free at the right end (point x
3
).
162
24
1 24
( , ) ( , )
( , ), ( , )
ii
i i i i i
w x t w x t
EI f x t x x x
tx
(7.15)
with the boundary conditions
1 0 1 0
23
3 3 3 3 3 3 23
( , ) 0, ( , ) 0
( , ) 0, ( , ) 0
w x t w x t
x
EI w x t EI w x t
xx
(7.16)
The matching conditions at the point x
1
and x
2
are
11
22
11 22
( , ) 0, ( , ) 0, ( , ) ( , )
( , ) ( , ), 1,2
i i i i i i i i
i i i i i i
w x t w x t w x t w x t
xx
EI w x t EI w x t i
xx
(7.17)
The governing equation, boundary and matching conditions are expressed in state
forms by Eqs. (7.8) and (7.9). The state matrix and boundary matrices are given by
1 2 3
( , ) ( , ) ( , ) ( , ) x s diag x s x s x s F F F F (7.18a)
11
11
22
,
GG
M 0 0 0 0 N
M 0 A 0 N B 0 0
0 0 A 0 B 0
(7.18b)
with
11
3
2
3
0 1 0 0
0 0 0 0 1 0 0 0
0 0 1 0
0 0 0 0 0 1 0 0
0 0 0 1 ( , ) , ,
0 0 0 0000
000 0 0 0 0000
i
i
i
xs
EI
s
EI
EI
F M N (7.19a)
163
1
1 0 0 0 0 0 0 0
0 1 0 0 0 1 0 0
, , 1,2
0 0 0 0 0 0
0 0 0 0 1 0 0 0
ii
ii
i
EI EI
AB (7.19b)
For numerical simulation, the system parameters are assigned as
1 1 1
2 2 2
3 3 3
2.0, 10, 1.0
1.5, 8, 0.7
1.2, 6, 1.3
EI L
EI L
EI L
(7.20)
The first 8 natural frequencies of the beam are computed by det ( ) 0
G
s Z and listed in
Table 7.1, where the proposed method is compared with the FEM. By Eq. (7.12), the first
four mode shapes are plotted in Figure 7.3. Assume that the stepped beam is subject to
harmonic excitation at the point ( 0.2)
f
x . The frequency response of the stepped beam
at the right end () xL is plotted against the frequency parameter ω.
164
(C1) Hinge (C2) Sliding constraint
k
r
k
t
k
r
I
D
(C3) Spring constraint (C4) Hinged disk
k
t
m
m,I
k
t
k
r
(C5) Sliding mass (C6) Mass and springs
Figure 7.1 Specification of constraints for stepped beams (Yang, 2005)
165
w(x,t)
x
x
0
(1) (2)
x
1 x
3
x
2
(3)
l
1
l
2
l
3
Figure 7.2 A stepped beam with hinged constraints
(a) (b)
(c) (d)
Figure 7.3 Mode shapes of the stepped beam: (a) mode 1;
(b) mode 2; (c) mode 3; (d) mode 4
166
Table 7.1 Natural frequencies of the stepped beam (rad/s)
k DTFM
Finite Element Method
N=30 N=60 N=150
1
2
3
4
5
6
7
8
3.8343
24.8808
39.0778
65.0575
84.0452
123.9619
150.6502
216.0138
3.8343
24.8810
39.0786
65.0611
84.0531
123.9878
150.6959
216.1451
3.8343
24.8808
39.0778
65.0577
84.0456
123.9635
150.6531
216.0221
3.8343
24.8808
39.0778
65.0575
84.0452
123.9619
150.6502
216.0140
Figure 7.4 The frequency response of the stepped beam at x = L:
amplitude (a) and phase (b)
167
7.3 Pointwise Connected Systems
A multi-body structure is often composed of subsystems which are connected
pointwise through some mechanism. The systems can be expressed by the governing
equation, boundary and matching conditions. The matching conditions are described by
the force equilibrium and geometric compatibility at the internodes. By the distributed
transfer function formulation, the state equations in a global domain have the same form
as Eqs. (2.45) and (2.46). The pointwise connection of subsystems is described through
modification of the matrices,
G
M and
G
N . See the following example in detail.
Example 7.2: A coupled beam structure
See Fig. 7.5. Consider a coupled Euler-Bernoulli beam with pointwise connections at
the nodes x
2
and x
3
. The beam is fixed at the points (x
1
and x
4
). The displacement ( , )
i
w x t
of the subsystems is governed by (Yang, 1994)
24
24
( , ) ( , ) ( , ), , 1,2
i i i i i i
w x t EI w x t f x t x i
tx
(7.21)
where ( , )
i
f x t is an external disturbance;
i
EI and
i
are the bending stiffness and linear
density, respectively.
The system is subject to the boundary conditions
1 1 1
23
2 1 2 1 2 3 2 23
: ( , ) 0, ( , ) 0
: ( , ) 0, ( , ) ( , ) ( , )
x x w x t w x t
x
x x EI w x t EI w x t kw x t kw x t
xx
(7.22a)
168
23
3 2 3 2 3 2 3 23
4 4 4
: ( , ) 0, ( , ) ( , ) ( , )
: ( , ) 0, ( , ) 0
x x EI w x t EI w x t kw x t kw x t
xx
x x w x t w x t
x
(7.22b)
Defining the state vector ˆ ( , )
i
xs , the Laplace transformed boundary conditions are given
11 1 1 12 2 3 11 1 2 1
22 2 3 21 1 2 22 2 4 2
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( )
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( )
b
b
x s x s x s s
x s x s x s s
M M N
M N N
(7.23a)
where
1
()
b
s and
2
()
b
s is prescribed by the boundary disturbances, and
23
23
ˆ ( , ) , 1,2
T
i i i
ii
w w w
x s w i
x x x
11 12 11
1
1
0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0
,,
0 0 0 0 0 0 0 0 0 0 0
00 0 0 0 0 0 0 0
EI
k EI k
M M N (7.23b)
2
2
22 21 22
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
,,
0 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0
EI
k EI k
M N N
The governing equation and boundary conditions are rewritten in state forms, Eq.
(2.45) and (2.46).
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ),
ˆ ˆ ( , ) ( , ) ( )
G s G e G
x s x s x s x s x
x
x s x s s
Fp
MN
(7.24a)
where
169
12
1 1
2 1
( , ) ( , ) ( , )
ˆ () ( , )
ˆ ( , ) , ( )
ˆ () ( , )
T T
b
G
T T
b
x s diag x s x s
s xs
x s s
s xs
F F F
p
p
p
(7.24b)
11 12 11
22 21 22
,
GG
M M N 0
MN
0 M N N
By det( ( , , )) 0,
G G e s
x x s MN the natural frequencies of the system are computed and
listed in Table 7.2. It is seen that the results match well with those by FEM.
170
w(x,t)
x
x
1
(1)
(2)
ρ=1, EI=1
1 1
k=1
x
2
x
3
x
4 ρ=1, EI=1
Figure 7.5 A coupled beam structure
Table 7.2 The natural frequencies of the coupled beam (in rad/s)
k DTFM
Finite Element Method
N=8 N=16 N=20
1
2
3
4
5
6
7
8
9
10
3.5160
4.4948
22.0345
22.2177
61.6972
61.7622
120.9019
120.9350
199.8595
199.8796
3.5161
4.4951
22.0602
22.2440
62.1749
62.2412
122.6576
122.6899
228.1374
228.1569
3.5160
4.4948
22.0363
22.2195
61.7347
61.7999
121.1728
121.2061
201.0159
201.0362
3.5160
4.4948
22.0352
22.2184
61.7129
61.7780
121.0171
121.0503
200.3633
200.3835
171
7.4 Distributed Connected Systems
(II)
(I)
Distributed
connection
Figure 7.6 A distributed connected system
Consider subsystems (I) and (II) connected throughout their domain; see Figure 7.6.
For example, a series of beams with this type of connections such as elastic layers can be
a model for the analysis of carbon nano-tubes. Assembly of the coupled subsystems leads
to coupled integral equations, which cannot be solved in exact and closed form. To
overcome this issue, a state equation for the coupled subsystems is presented in Ref.
(Yang, 1994).
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ) ( , ), (0,1)
i i i i ci
x s x s x s x s x s x
x
F p Q (7.25a)
ˆ ˆ (0, ) (1, ) ( )
i i i i
s s s MN (7.25b)
where ( , )
ci
xs Q and ( , )
ci
xs Q are distributed forces applied at the subsystem (I) and (II), it
requires to maintain the connection in Fig. 7.6.
Assume that the connection is uniformly distributed; the distributed forces are given as
1 11 12 1
2 21 22 2
ˆ ( , ) ( ) ( ) ( , )
ˆ ( , ) ( ) ( ) ( , )
c
c
xs s s x s
xs s s x s
Q CC
Q CC
(7.26)
172
where ()
ij
s C are the matrix transfer functions describing the connection. Substitution of
(7.26) into (7.25) yields
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , )
ˆ ˆ ( , ) ( , ) ( )
G s G e G
x s x s x s x s
x
x s x s s
Fp
MN
(7.27a)
with
1 11 12
21 2 22
1 1 1
2 2 2
11
22
( , ) ( ) ( )
( , )
( ) ( , ) ( )
ˆ ˆ ( , ) ( , ) ( )
ˆ ˆ ( , ) , ( , ) ( )
ˆ ˆ ( , ) ( , ) ( )
,
G
GG
x s s s
xs
s x s s
x s x s s
x s x s s
x s x s s
F C C
F
C F C
p
p
p
M 0 N 0
MN
0 M 0 N
+
+
(7.27b)
Because the state equation of the coupled subsystems has the same form as that of an
uncoupled subsystem, the eigensolutions and response can be determined by the same
solution procedures of an uncoupled subsystem.
Example 7.3: Elastically connected beams
In Figure 7.7, two Euler-Bernoulli beams are connected through an elastic layer. When
the lengths of two beams are the same, i.e.
12
ll , the system is regarded as the assembly
of two substructures. Let ( , ), 1,2
i
w x t i represent the transverse displacement of the ith
beam which has the inertia parameter
i
, stiffness parameter
i
EI and the stiffness k per
length between the beams. Both beams are fixed at the left end and free at the right end.
The transverse displacement of the 1
st
and 2
nd
beam is governed by (Kelly and Srinivas,
2009)
173
24
1 1 1 2 1 1 1 24
24
2 2 2 1 2 2 2 24
( , ) ( ( , ) ( , )) ( , ) ( , )
( , ) ( ( , ) ( , )) ( , ) ( , )
w x t k w x t w x t EI w x t f x t
tx
w x t k w x t w x t EI w x t f x t
tx
(7.28)
Assume zero initial conditions. Taking Laplace transform to Eq. (7.28), the equation
and boundary conditions are expressed in a state form, Eq. (7.27) with
1 11 12
21 2 22
( ) ( ) ( )
()
( ) ( ) ( )
s s s
s
s s s
F C C
F
C F C
(7.29a)
with
2
0 1 0 0
0
0 0 1 0
0
0 0 0 1 ( ) , ( , ) for 1,2
0
000 1
i
ii
i
i
i
f
s x s i
EI
s
EI
Fp
(7.29b)
11 22
12
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
( ) , ( ) 0 0 0 0 0 0 0 0
0 0 0 0 0 0
ss
kk
EI EI
CC (7.29c)
12 11 21 22
( ) ( ), ( ) ( ) s s s s C C C C
For the numerical simulation, the system parameters are assigned as
1 1 1
2 2 2
0.5, 20, 1
0.5, 25, 1, 200
EI l
EI l k
(7.30)
The eigenvalue of the system are computed by det( ( , , )) 0
G G e s
x x s MN . The natural
eigenvalues of the system are listed in Table 7.3. For comparison, the results by FEM are
listed in the table.
174
Denote the change of the length
12
l l l where
12
ll . The length of an elastic
foundation is the same as that of the second beam. If the lengths of two beams are not the
same, i.e.
12
ll , the system is regarded as being composed of three substructures where it
is connected through a partial elastic foundation. The eigenvalue loci of the system are
plotted versus the change of the length l in Fig. 7.8. The loci are zoomed into the small
scales in Fig. 7.9. As the stiffness of the elastic foundation approaches to zero (k→0), the
system can be considered as two uncoupled subsystems. Eigenvalue loci of uncoupled
subsystems are plotted for comparison (dashed-curves in the figure). They form as
bounds for natural frequencies of the coupled system. Curve veering is observed in exact
solutions associated with structures connected through an elastic foundation according to
l . Eigenvalue loci approach closely and rapidly veer away. At the veering points, curve
veering indicate coupling between the eigenfunctions with regard to system parameter.
175
k
EI
1
, ρ
1
EI
2
, ρ
2
l
2
l
1
Figure 7.7 An elastically connected beam
Table 7.3 The natural frequencies of the elastically connected beam (in rad/s)
k DTFM
Finite Element Method
N=8 N=16 N=20
1
2
3
4
5
6
7
8
9
10
23.4853
36.8925
140.6699
157.1899
390.7143
436.7282
764.9117
855.1402
1264.1807
1413.3619
23.4860
36.8930
140.8310
157.3696
393.7316
440.1023
776.0121
867.5517
1443.0061
1613.2991
23.4853
36.8925
140.6810
157.2023
390.9513
436.9933
766.6240
857.0548
1271.4934
1421.5380
23.4853
36.8925
140.6745
157.1950
390.8135
436.8392
765.6401
855.9547
1267.3666
1416.9240
176
Figure 7.8 Loci of the natural frequencies versus the change of the beam length
12
l l l ; solid-curve −: beams with an elastic foundation,
dash-curve −: the uncoupled subsystems
(a) (b)
Figure 7.9 loci of the natural frequencies (the extended view in boxes in Fig. 7.8)
177
7.5 Combined Systems
y
1
x
1
x
2
x
n-1
Distributed system
Connection (C
w
, C
v
)
•
•
•
•
•
•
y
2
y
n-1
z
1
z
2
z
p
• • •
• • •
Lumped system
Figure 7.10 Connection of a distributed subsystem
and a lumped parameter system
Figure 7.10 shows a lumped system is connected to the distributed system at multiple
points. Let the vector of the displacements of the distributed system at point
k
x be () s w ,
and the vector of the displacements of the lumped system at point
k
y be () s v , k = 1,2, …,
n-1. Denote the vector of the constraint forces at point
k
x be () s Q . The constraint forces
at points
k
x are related to the displacements at the connecting points, and described as
(Yang, 1994)
( ) ( ) ( ) ( ) ( )
vw
s s s s s Q C v C w (7.31)
where ()
v
s C and ()
w
s C are known matrices describing the connection.
With the given transformation matrix
v
T , the relationship between () s v and () s u can
be described by
178
( ) ( )
v
ss v T u = (7.32)
Let the vector of the displacements of the lumped system at point
k
z be ()
k
us . In the
Laplace transform domain, the dynamic equation of the lumped system is in the matrix
form
( ) ( ) ( ) ( )
u
s s s s A u T Q f (7.33)
where () s A and T are given matrices, and ()
u
s f is the vector of the external forces
applied to the lumped system.
Substitution Eq. (7.32) to Eq. (7.31) and eliminating () s u with Eq. (7.33) give
( ) ( ) ( ) ( )
T
s s s s Q C w f (7.34)
where I is identity matrix and
1
1
1
11
( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
v v w
T v v v v u
s s s s
s s s s s s
C C T A T C
C T A T C T A
I
f I f
(7.35)
For the connection of a lumped system and a distributed system at a single point, the
matrix forms of Eqs. (7.31) - (7.35) are expressed as scalars.
In general, the matching conditions are described by the geometric compatibility and
force equilibrium at the internodes. With the constraint forces obtained by Eq. (7.34), the
force balances at the point
k
x can be determined. The matching conditions are needed to
consider only with the force balances. In the state forms Eqs. (2.45) and (2.46), such
connection of subsystems is described through modification of the matrices and
GG
MN .
The solutions of the system are determined by the same procedures as seen in chapter 7.2.
179
The proposed method is applied to determine eigensolutions of combined beams in the
following examples.
Example 7.4: A combined beam with an oscillator
Consider a simply supported Euler-Bernoulli beam connected to a mass-spring lumped
system atxa ; see Fig. 7.11. The combined beam can be regarded as two subsystems (I)
and (II) serially connected at the internode (x
1
). The combined beam is the same as the
stepped beams with the constraint. Thus, the state equation can be expressed by Eqs. (7.8)
and (7.9). Let the transverse displacement of two subsystems (I) and (II) be
1
( , ) w x t and
2
( , ) w x t , respectively. In the Laplace domain, the matching conditions at the internode is
given as
1 1 2 1 1 1 2 1
22
2 1 1 1 22
33
2 1 1 1 33
( , ) ( , ), ( , ) ( , )
( , ) ( , )
( , ) ( , ) ( )
w x s w x s w x s w x s
xx
EI w x s EI w x s
xx
EI w x s EI w x s Q s
xx
(7.36)
where the constraint force () Qs , by Eqs. (7.34) and (7.35), is
2
11 2
( ) ( , )
mks
Q s w x s
ms k
(7.37)
with
2
( ) , ( ) and ( )
vw
A s ms C s k C s k .
Let the system parameters be 1, 1, 2 EI L . The parameters of the oscillator are
selected as 32, 60, 1 m k a . By Eq. (7.14), the natural frequencies of the system are
computed and listed in Table 7.4 and compared with the results by FEM. Figure 7.12
180
shows the loci of the first three natural frequency as a function of the a spring of the
stiffness k. It is reported that curve veering usually occurs in spatially asymmetric
systems with intermediate elastic constraints. However, curve veering also exists on the
spatially symmetric beam with an oscillator as shown in Fig. 7.12. Furthermore, it is
shown that eigenfunctions associated with the eigenvalues on each locus are interchanged
in a rapid but continuous way at the veering. The switch of eigenfunctions is an important
characteristic of curve veering.
181
EI, ρ
L
m
k
x
0
x
2
x
1
x=a
(I) (II)
u(t)
Figure 7.11 A beam connected with an oscillator
182
Table 7.4 The natural frequencies of the beam with an oscillator (in rad/s)
k DTFM
Finite Element Method
N=8 N=16 N=20
1
2
3
4
5
6
7
8
9
10
0.4078
7.6937
9.8696
23.6385
39.4784
62.1818
88.8264
121.1530
157.9137
200.0103
0.4078
7.6949
9.8722
23.6730
39.6342
62.7683
90.4495
24.9575
175.2712
218.3350
0.4078
7.6938
9.8698
23.6408
39.4887
62.2215
88.9407
121.4390
158.5369
201.2540
0.4078
7.6937
9.8697
23.6395
39.4826
62.1982
88.8739
121.2723
158.1753
200.5357
Figure 7.12 Loci of the eigenvalues versus the coefficient of a spring k
183
Example 7.5: A beam with elastically mounted rigid body
The proposed method is applied to eigenvalue problem of an Euler-Bernoulli beam
mounted with a rigid body; see Figure 7.13. The combined beam can be considered the
assembly of serially connected three beams (I), (II) and (III). The combined beam is fixed
at the both ends. The transverse displacement ( , )
i
w x t of the subsystems (I)-(III) is
governed by (Yang, 2005)
24
24
( , ) ( , ) ( , ), 1,2,3
i i i
w x t EI w x t f x t i
tx
(7.38)
The motion of the rigid body is described by Eqs. (7.31), (7.32) and (7.33)
1 1 1 1 1 1
2 2 2 2 2 2
( ) 0 ( ) 0 ( , )
( ) 0 ( ) 0 ( , )
F s k z s k w x s
F s k z s k w x s
(7.39a)
1
2
() 1 ( )
() 1 ( )
g
zs b z s
zs a b s
(7.39b)
2
1
2
2
() ( ) 1 1 0
() () 0
g
Fs zs Ms
Fs s b a b Is
(7.39c)
The matching conditions at internodes
12
( , ) xx are
11
22
122
33
133
( , ) ( , ), ( , ) ( , )
( , ) ( , )
( , ) ( , ) ( )
i i i i i i i i
i i i i
i i i i i
w x t w x t w x t w x t
xx
EI w x t EI w x t
xx
EI w x t EI w x t F t
xx
(7.40)
for 1,2. i
184
By Eq. (7.34) and (7.39), the constraint forces are expressed by the displacements
terms of the beam only. Then, substitution it into Eq. (7.40) in the Laplace transform
domain gives simplified matching conditions with the displacements terms of the beam.
The state equation of the system has the same form as that of stepped beams, Eqs. (7.8)
and (7.9). For a numerical simulation, the system parameters are assigned as
1 1 2 2
1, 1, 2
0.4, 0.1, 0.05
0.9, 10, 0.2, 10
EI L
M I b
l k l k
(7.41)
The eigenvalues of the system are computed by Eq. (7.14) and listed in Table 7.5. The
results are compared with those by the FEM.
185
x
1
x
2
z
1
z
2
G
x
0
x
3
M, I b
k
1
k
2
(I)
(II)
(III)
Figure 7.13 A beam with elastically mounted rigid body
Table 7.5 The natural frequencies of the beam mounted a rigid body (in rad/s)
k DTFM
Finite Element Method
N=20 N=40 N=100
1
2
3
4
5
6
7
8
9
10
1.3935
4.2316
9.3092
15.5188
30.4819
50.0489
74.6964
104.3174
138.8019
178.3228
1.3935
4.2316
9.3092
15.5191
30.4838
50.0575
74.7248
104.3942
138.9811
178.6986
1.3935
4.2316
9.3092
15.5188
30.4820
50.0495
74.6982
104.3223
138.8134
178.3471
1.3935
4.2316
9.3092
15.5188
30.4819
50.0489
74.6964
104.3176
138.8022
178.3234
186
7.6 Truss/Frame Structures
u
w
X
Y
L
α
x
y
Figure 7.14 Schematic of a frame member
Consider an elastic plane frame member consisting of n prismatic elements that can
sustain transverse and longitudinal deformations. Let X and Y be the local coordinates.
The longitudinal displacement ( , ) u x t and transverse displacement ( , ) w x t of such a
member, as shown in Fig. 7.14, are governed by the differential equations (Chang, 2010)
22
22
24
24
( , ) ( , )
( , ), (0, )
( , ) ( , )
( , )
a
t
u x t u x t
EA f x t x L
tx
w x t w x t
EI f x t
tx
(7.42)
where
is the inertia parameter; EA and EI are the longitudinal stiffness and bending
stiffness of a member, respectively; L is the member length; ( , )
a
f x t is an axial external
load; and ( , )
t
f x t
is a transverse external load. The boundary conditions have general
forms
187
2
1 0 1 2
23
2 3 2 23
2
4 5 3 2
at 0 : ( )
()
()
LL
LL
L
uu
x m a a u t
tx
ww
m a w a t
tx
ww
a a t
xx
(7.43a)
2
1 0 1 2
23
2 3 2 23
2
4 5 3 2
at : ( )
()
()
RR
RR
R
uu
x L m b b u t
tx
ww
m b w b t
tx
ww
b b t
xx
(7.43b)
where
L
m and
R
m are end masses; ()
Lj
t and ()
Rj
t ( 1,2,3) j are related to prescribed
boundary excitations; and
i
a and
i
b ( 0,1,2,3) i are constants that are properly assigned
to characterize different types of boundaries. The member is also subject to the initial
conditions
00
( ,0)
( ,0) ( ), ( ), 0,
aa
ux
u x u x v x x L
t
(7.44a)
00
( ,0)
( ,0) ( ), ( )
tt
wx
w x u x v x
t
(7.44b)
where
0 0 0
( ), ( ), ( )
a a t
u x v x u x and
0
()
t
vx are given profiles of initial states of the member.
By the distributed transfer function formulation, the governing equation, boundary and
initial conditions can be expressed as a state form in Laplace transform domain
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ), (0, ) x s x s x s x s x L
x
Fp (7.45a)
ˆ ˆ (0, ) ( , ) ( )
b b b
s L s s MN
(7.45b)
188
where
2
2
0 1 0 0 0 0
/ 0 0 0 0 0
0 0 0 1 0 0
( , )
0 0 0 0 1 0
0 0 0 0 0 1
0 0 / 0 0 0
s EA
xs
s EI
F (7.46a)
23
23
ˆ ( , )
T
i i i i
ii
u w w w
x s u w
x x x x
(7.46b)
The matrices
b
M
and
b
N
contain those of the boundary operators in Eq. (7.43); ˆ( , ) xs p
and ()
b
s are vectors to contain the external and the boundary disturbances applied to the
frame with the initial states, respectively.
For the nodes where frame components interconnected, a coordinate transformation
matrix for each element is used to analyze nodal displacements and forces in a global
coordinate system because the matching conditions require geometric conditions and
force balances at the nodes where adjacent elements interconnected. In this section, a
frame component is described by partial differential equations, initial conditions and
boundary conditions. The proposed method in Section 2.4 is applicable to multi-body
structures with any number of components and arbitrary boundary conditions, and is
capable to obtaining exact solutions of multi-body flexible structures. The method is
illustrated on the frame structure composed of two elements.
189
Example 7.6: Truss/Frame structure
As shown in Fig. 7.15, consider the frame composed of two elements which have the
inertia parameter
i
, the longitudinal stiffness EA
i
, bending stiffness EI
i
and the same
length L. The longitudinal displacement ( , )
ii
u X t and transverse displacement ( , )
ii
w X t of
the ith element are governed by
22
,1 22
24
, 24
( , ) ( , ) ( , ), ( , ), 1,2
( , ) ( , ) ( , )
i i i i i i a i i i i i
i
i i i i i i t i i
i
u X t EA u X t f X t X x x i
tX
w X t EI w X t f X t
tX
(7.47)
Assume that the frame is rigidly connected at the nodes
02
and xx . The system is subject
to the boundary conditions
1 0 1 0 1 0
1
2 2 2 2 2 2
2
( , ) 0, ( , ) 0, ( , ) 0
( , ) 0, ( , ) 0, ( , ) 0
u x t w x t w x t
X
u x t w x t w x t
X
(7.48)
The central joint, i.e. the node
1
() x , is rigidly connected. The matching conditions for the
displacements and internal forces at the node
1
x , are given as
1 2 2
1 2 2
cos sin
sin cos
u u w
w u w
(7.49a)
3
1 2 2
1 2 2 3
1 2 2
33
1 2 2
1 2 2 33
1 2 2
cos sin 0
sin cos 0
u u w
EA EA EI
X X X
w u w
EI EA EI
X X X
(7.49b)
22
1 2 1 2
12 22
1 2 1 2
,
w w w w
EI EI
X X X X
(7.49c)
190
By the distributed transfer function formulation, the state equation of each element is
1
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ), ( , ), 1,2
i i i i i i i i i i i
i
X s X s X s X s X x x i
X
Fp (7.50a)
where
23
23
2
2
ˆ ( , )
0 1 0 0 0 0
/ 0 0 0 0 0
0 0 0 1 0 0
( , )
0 0 0 0 1 0
0 0 0 0 0 1
0 0 / 0 0 0
T
i i i i
i i i i
i i i i
ii
ii
ii
u w w w
X s u w
X X X X
s EA
Xs
s EI
F
(7.50b)
The boundary and the matching conditions are
1 1 0 1 2 2 1
1 1 1 1 2 1 2
ˆ ˆ ˆ ( , ) ( , ) ( )
ˆ ˆ ˆ ( , ) ( , ) ( )
x s x s s
x s x s s
MN
BA
(7.51)
Define the global domain
0 1 1 2
( , ) ( , ) x x x x U . By the transfer function formulation in
Section 2.4, Eqs. (7.50) and (7.51) are rewritten as
ˆ ˆ ˆ ( , ) ( , ) ( , ) ( , ), x s x s x s x s x
x
Fp (7.52a)
ˆ ˆ ( , ) ( , ) ( )
G s G e G
x s x s s MN (7.52b)
where
s
x and
e
x represent the end nodes of a component, and
1 1 1
2 2 2
12
ˆ ˆ ( , ) ( , )
ˆ ˆ ( , ) , ( , )
ˆ ˆ ( , ) ( , )
( , ) ( , ) ( , )
TT
TT
X s x s
x s x s
X s x s
x s diag x s x s
p
p
p
F F F
(7.53a)
191
11
11
,
GG
M 0 0 N
MN
0 A B 0
(7.53b)
with
36
11
36
1 0 0 0 0 0
0 0 1 0 0 0 1 0 0 0 0 0
,
0 0 0 1 0 0 0 0 1 0 0 0
0 0 0 1 0 0
0
MN
0
(7.54a)
22
1
22
2
cos 0 sin 0 0 0
sin 0 cos 0 0 0
0 cos 0 0 0 sin
0 sin 0 0 0 cos
0 0 0 1 0 0
0 0 0 0 0
EA EI
EA EI
EI
A (7.54b)
1
1
1
1
1 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 1 0 0
0 0 0 0 0
EA
EI
EI
B (7.54c)
For a numerical simulation, the system parameters are assigned as follows:
7
1 1 1
6
2 2 2
78, 2 10 , 2000, 1
27, 7 10 , 700,
4
EA EI L
EA EI
(7.55)
The first 8 natural frequencies of a frame structure composed of two members are
computed by det( ( , , )) 0
G G e s
x x s MN and listed in Table 7.6, where the proposed
method is compared with the FEM. The corresponding eigenfunction (mode shape) is
192
obtained by ( ) ( ,0, )
k k k
x x s a where
k
a
is a nonzero solution of the homogeneous
equation ( , , ) 0
G G e s
x x s M N a , and are plotted in Figure 7.16.
X
Y
α
x
2
x
0
x
1
X
2
X
1
Figure 7.15 A frame structure
Table 7.6 The natural frequencies of the frame structure (rad/s)
k DTFM
Finite element method
N=8 N=16 N=20
1
2
3
4
5
6
7
8
77.8560
110.8257
247.1231
287.2990
452.9538
542.4383
706.8876
789.4352
77.9048
110.9655
248.6106
289.5352
460.8170
554.2979
720.7471
795.8150
77.8591
110.8346
247.2210
287.4496
453.5092
543.4245
708.9127
790.7941
77.8572
110.8293
247.1635
287.3614
453.1895
542.8512
707.7606
790.2631
193
(a) (b)
(c) (d)
Figure 7.16 Mode shapes of the frame structure: (a) mode 1;
(b) mode 2; (c) mode 3; (d) mode 4
194
Chapter 8. Conclusions
The DTFM-based method has been developed for analysis and control of general one-
dimensional distributed parameter systems. Through use of system transition matrix, the
method is generalized and implemented for an eigenvalue problem and frequency
response of general distributed systems, especially for non-uniform systems. With the
distributed transfer function formulation and a residue theorem for inverse Laplace
transform, the method gives an exact and closed-form transient solution of the system
subject to arbitrary external, initial and boundary disturbances in a systematic manner.
Furthermore, the method does not depend on system eigenfunctions for analysis in
frequency and time domain. The key in the proposed method is to determine the state
transition matrix and transfer function residues required to estimate an exact time-domain
response. By the augmented formulation and Jacobi’s formula, the explicit formula is
represented for evaluation of the state transition matrix and the transfer function residues
needed in the transient analysis as has been shown in section 2.3.
A non-uniformly distributed system is described by the partial differential equation
whose coefficients depend explicitly on the spatial variable x, initial conditions over the
domain of the system and boundary conditions at the end points of the domain. Exact
solutions are limited to several cases of non-uniform systems. If exact solutions are
possible, the exact state transition matrix which plays important roles in the proposed
analysis can be determined. With the state transition matrix, the proposed solution
195
method delivers an exact and closed-form solution of a non-uniformly distributed system
in time and frequency domains. For arbitrary non-uniform systems whose exact solutions
are not available, the new technique to evaluate the approximate state transition matrix is
introduced. It is shown that the technique delivers accurate solutions for eigenvalue
problems of a non-uniform system.
The proposed method is different from existing series solution methods in two major
aspects. First, unlike eigenfunction expansion method, the method in transient analysis
does not depend on system eigenfunctions. Also in computation, it does not require
dealing with spatial integrals for normalization of system eigenfunctions. Second, many
series solution methods need different derivations for different system configurations.
The method adopts a symbolic formulation which systematically treats different
configurations, thus does not require tedious derivations. The DTFM-based method is not
only efficient in numerical computation but also convenient in computer coding.
In mid- and high-frequency analyses, it is apparently seen that if the spectrums of
external disturbances are narrow-banded, a few terms in the Green’s functions (selection
of modes) are enough to deliver accurate transient solutions. This feature makes the
method potentially useful for analysis of dynamic behaviors of distributed systems in
mid- and high-frequency ranges, where numerical methods may not be applicable.
The DTFM-based method is valid for both self-adjoint and non-self-adjoint systems.
Unlike existing series solution methods, the transient solution method does not need the
knowledge of system eigenfunctions. This thrust of the method is useful in time-domain
analysis. Avoiding cumbersome computations related with eigenfunctions, the method
196
gives exact and closed-form transient solutions of distributed systems with various effects,
such as a spinning beam and an axially loaded and moving beam.
This work presents the new distributed transfer function formulation for modeling and
analysis of complex distributed parameter systems that are composed of multiple
continua combined with lumped parameter systems and constraints. This formulation
makes the method applicable to complex distributed parameter systems without losing
any advantages of the DTFM-based method. Thus, the proposed method is extended for
analysis of general complex distributed systems. The method is especially capable of
handling systems with multiple branches and closed loops, distributed disturbances, and
multipoint connection of distributed and lumped systems, which otherwise may be
difficult to deal with by other methods. Furthermore, the method gives exact solutions of
truss/frame structures in longitudinal and transverse vibration.
In analysis of complex distributed parameter systems, the 'curve veering' phenomenon
is found to occur when exact solutions are available. The curve veering is represented for
two distributed systems; elastically connected beams and a combined beam with
oscillators. The curve veering of elastically connected beams is observed and much
influenced by spatially asymmetry. The curve veering is obsevered in a combined beam
with an oscillator and influenced by an intermediate elastic constraint. It is seen that
spatially asymmetry and elastic constraints are the cause of curve veerings of distributed
systems as have been shown in sections 7.4 and 7.5.
The eigenvalues of a distributed system are calculated by finding the roots of the
characteristic equation which is given by the transcendental form. Differently from the
197
spectral element model (Lee et al., 2002; Lee, 2009) and the finite element formulation,
the condensed characteristic equation is given by the DTFM-based method. By
introduction of specific bounds of eigenvalues, the root-finding technique is developed
for the determination of eigenvalues of stepped systems, such as stepped bars, shafts,
strings and beams combined with lumped masses and constraints. The distribution of
eigenvalues of a complex distributed parameter system is in general irregular. In such
cases, the proposed technique alleviates the risk of intensive computation rate as well as
the possibility of missing eigenvalues.
In this investigation, the boundary vibration control of non-uniform bars, shafts and
strings has been presented. By the distributed transfer function formulation, an exact and
closed-form solution of the controlled system can be determined. The active noise control
of ducts with axial temperature gradients has been presented. A closed-form expression
of the feedback control is given in this research. It is seen that the destabilizing effect of
non-collocation of sensors and actuators can be eliminated with the introduction of
specific time-delay relations. Although this study has focused on ducts of two different
types of boundary conditions, the distributed transfer function concept and the time delay
approach can be utilized in duct systems of other types of boundary conditions.
198
References
Abrate, S., 1995, Vibration of non-uniform rods and beams, Journal of Sound and
Vibration, 185, pp. 703-716.
Aristizabal-Ochoa, J. D., 2004, Timoshenko beam-column with generalized end
conditions and nonclassical modes of vibration of shear beams, Journal of Engineering
Mechanics, 130(10), pp. 1151-1159.
Aristizabal-Ochoa, J. D., 2008, Slope-deflection equations for stability and second-order
analysis of Timoshenko beam-column structures with semi-rigid connections,
Engineering Structures, 30, pp. 2517-2527.
Bapat, C. N. and Bhutani, N., 1994, General approach for free and forced vibrations of
stepped systems governed by the one-dimensional wave equation with non-classical
boundary conditions, Journal of Sound and Vibration, 172, pp. 1-22.
Bapat, C. N., 1995, Vibration of rods with uniformly tapered sections, Journal of Sound
and Vibration, 185, pp. 185-189.
Bergman, L. A. and McFarland, D. M., 1988, On the vibration of a point supported linear
distributed system, ASME Journal of Vibration, Acoustics, Stress, and Reliability in
Design, 110, pp. 485-492.
Bergman, L. A. and Nicholson, J. W., 1985, Forced vibration of a damped combined
linear system, Journal of Vibration, Acoustics, Stress, and Reliability in Design, 107, pp.
275-281.
Blaauwendraad, J., 2008, Timoshenko beam-column buckling. Does Dario stand the test?,
Engineering Structures, 30, pp. 3389-3393.
Bottega, W. J., 2006, Engineering Vibrations, CRC Taylor & Francis, Florida
Brown, J. W. and Churchill, R. V., Mei, C. C., 2003, Complex Variables and
Applicatoins, McGraw Hill, New York.
Butkovisky, A. G., 1983, Structural theory of distributed systems, Halsted press, John
Wiley and Sons, New York.
199
Caughey, T. K., O’Kelley, M. E. J., 1965, Classical normal modes in damped linear
systems, ASME Journal of Applied Mechanic, 32, pp. 583-588.
Chang, C. H., 2010, Mechanics of Elastic Structures with Inclined Members: Analysis of
Vibration, Buckling and Bending of X-Braced Frames and Conical Shells, Springer, New
York.
Chen, C.T., 1999, Linear Systems Theory and Design, 3rd ed., Oxford University Press,
New York.
Chen, Y., 1963, On the vibration of beams or rods carrying a concentrated mass, ASME
Journal of Applied Mechanics, 30, pp. 310-311.
Cummings, A., 1977, Ducts with axial temperature gradients: an approximate solution for
sound transmission and generation, Journal of Sound and Vibration, 51, pp. 55-67.
Cummings, A., 1978, Sound generation and transmission in flow ducts with axial
temperature gradients, Journal of Sound and Vibration, 57, pp. 261-279.
Dimarognas, S. A. and Paipeties, S. A., 1983, Analytical Methods in Rotor Dynamics,
Applied Science Publishers, New York.
Dokumaci, E., 1998, An approximate analytical solution for plane sound wave
transmission in inhomogeneous ducts, Journal of Sound and Vibration, 217, pp. 853-867.
Dokumaci, E., 1998, An exact transfer matrix formulation of plane sound wave
transmission in inhomogeneous ducts, Journal of Sound and Vibration, 217, pp. 869-882.
Eisenberger, M., 1991, Exact longitudinal vibration frequencies of a variable cross-
section rod, Applied Acoustics, 34, pp. 123-130.
Eriksson, L. J., Allie, M. C. and Greiner, R. A, 1987, The selection and application of an
IIR adaptive filter for use in active sound attenuation, IEEE Transactions on. Acoustic,
Speech and Signal Processing, 35, pp. 433-437
Fedotov, I., Joubert, S., Marais, J. and Shatalov, M., 2006, Another approach to vibration
analysis of stepped structures, Electronic Transactions on Numerical Analysis, 24, pp.
66-73.
Fedotov, I., Joubert, S., Marais, J. and Shatalov, M., 2008, Analysis for an N-stepped
Rayleigh bar with sections of complex geometry, Applied Mathematical Modeling, 32,
pp. 1-11.
200
Gonçalves, P.J.P., Brennan, M. J. and Elliottt, S. J., 2007, Numerical evaluation of high-
order modes of vibration in uniform Euler-Bernoulli beams, Journal of Sound and
Vibration, 301, pp. 1035-1039.
Gorman, D. J., 1975, Free Vibration Analysis of Beams and Shafts, Wiley, New York.
Han, S. M., Benaroya, H. and Wei, T., 1999, Dynamics of transversely vibrating beams
using four engineering theories, Journal of Sound and Vibration, 225, pp. 935-988.
Hildebrand, F. B., 1964, Advanced Calculus for Applications, Prentice-Hall, Englewood
Cliffs, NJ.
Hsueh, W. J., 1999, Free and forced vibrations of stepped rods and coupled systems,
Journal of Sound and Vibration, 226, pp. 891-904.
Hu, J. S., 1995, Active sound attenuation in finite length ducts using close-form transfer
function models, ASME Journal of Dynamic Systems, Measurement, and Control, 117,
pp. 143-154.
Hull, A. J., Radcliffe, C. J. and Southward, S. C., 1990, Global active noise control of a
one-dimensional acoustic duct using a feedback controller, Journal of Vibration and
Acoustics, 115, pp. 488-494.
Kapur, A., Cummings, A. and Mungur, P., 1972, Sound propagation in a combustion can
with axial temperature and density gradients, Journal of Sound and Vibration, 25, pp.
129-138.
Karthik, B., Kumar, B. M. and Sujith, R. I., 2000, Exact solutions to one-dimensional
acoustic fields with temperature gradient and mean flow, Journal of the Acoustical
Society of America, 108, pp. 38-43.
Katz, R., Lee, C. W., Ulsoy, A. G., and Scott, R. A., 1988, The dynamic response of a
rotating shaft subject to a moving load, Journal of Sound and Vibration, 122, pp. 131-148.
Kelly, S. G. and Srinivas, S., 2009, Free vibrations of elastically connected stretched
beams, Journal of Sound and Vibration, 326, pp. 883–893.
Kounadis, A. N., 1980, On the derivation of equations of motion for a vibration
Timoshenko column, Journal of Sound and Vibration, 73(2), pp. 177-184.
Kumar, B. M. and Sujith, R. I., 1997, Exact solution for one-dimensional acoustic fields
in ducts with a quadratic mean temperature profile, Journal of the Acoustical Society of
America, 101, pp. 3798-3799.
201
Kumar, B. M. and Sujith, R. I., 1998, Exact solution for one-dimensional acoustic fields
in ducts with polynomial mean temperature profiles, Journal of Vibration and Acoustics,
120, pp. 965-969.
Kumar, B. M. and Sujth, R. I., 1997, Exact solutions for the longitudinal vibration of
non-uniformed rods, Journal of Sound and Vibration, 207, pp. 721-729.
LaFontaine, R. F. and Shepherd, I. C., 1983, An experimental study of a broadband active
noise attenuator for cancellation of random noise in ducts, 91(3), pp. 351-362
Langley, R. S., 1991, Analysis of beam and plate vibrations by using the wave equation,
Journal of Sound and Vibration, 150, pp. 47-65.
Langley, R. S., 1995, Mode localization up to high frequencies in coupled one-
dimensional subsystems, Journal of Sound and Vibration, 185, pp. 79-91.
Langley, R.S. and Bardell, N.S., 1998, A review of current analysis capabilities
applicable to the high frequency vibration prediction of aerospace structures, The
Aeronautical Journal, 102, pp. 287-297.
Lee, U., 2009, Spectral Element Method in Structural Dynamics, John Wiley and Sons,
New York.
Lee, U., Kim, J. and Oh, H., 2004, Spectral analysis for the transverse vibration of an
axially moving Timoshenko beam, Journal of Sound and Vibration, 271, pp. 685-703.
Lee, U., Kim, J. Shin, J. and Leung, A.Y.T., 2002, Development of a Wittrick-Williams
algorithm for the spectral element model of elastic-piezoelectric two-layer active beams,
International Journal of Mechanical Sciences, 44, pp. 305-318.
Li, H. and Stone, B. J., 1999, Time domain modeling of a reciprocating engine, Journal
of Mechanical Systems and Signal Processing, 13, pp. 169-178.
Li, Q. S., 2000, Exact solution for free longitudinal vibrations of non-uniform rods,
Journal of Sound and Vibration, 234, pp. 1-19.
Li, Q. S., 2003, Torsional vibration of multi-step non-uniform rods with various
concentrated elements, Journal of Sound and Vibration, 260, pp. 637–651.
Li, Q. S., Fang, J. Q. and Jeary, A. P., 1998, Calculation of vertical dynamic
characteristics of tall buildings with viscous damping, Int. Journal of Solids Structures,
35, pp. 3165-3176.
202
Li, Q., Cao, H. and Li, G., 1996, Static and dynamic analysis of straight bars with
variable cross section, Computers & Structures, 59, pp. 1185-1191.
Lin, Y. K., 1969, Dynamics of beam type periodic structures, ASME Journal of
Engineering for Industry, Nov., pp. 1133-1141.
Lueg, P., 1936, Process of silencing sound oscillations, Patent No. 2,043,416.
Mead, D. J. and Yaman, Y., 1990, The harmonic response of uniform beams on multiple
linear supports: a flexural wave analysis, Journal of Sound and Vibration, 141, pp. 465-
484.
Mead, D. J., 1971, Vibration response and wave propagation in periodic structures,
ASME Journal of Engineering for Industry, Aug., pp. 783-792.
Mead, D. J., 1994, Waves and modes in finite beams: application of the phase-closure
principle, Journal of Sound and Vibration, 171, pp. 695-702.
Mei, C. C., 1997, Mathematical Analysis in Engineering, Cambridge University Press,
New York
Meirovitch, L., 1967, Analytical Methods in Vibrations, Macmillan, New York.
Meirovitch, L., 1980, Computational Methods in Structural Dynamics, Sijthoff &
Noordhoff , Rockville, Maryland.
Meirovitch, L., 2001, Fundamentals of Vibrations, McGraw-Hill, New York.
Mioduchowski, A., 1986, Torsional waves and free vibrations of drive systems with
stepped shafts, Ingenieur Archiv, 56, pp. 314-320.
Moreles, M. A., Botello, S. and Salinas, R., 2005, A root-finding technique to compute
eigenfrequencies for elastic beams, Journal of Sound and Vibration, 284, pp. 1119-1129.
Morse P. M. and Feshbach H., 1953, Methods of Theoretical Physics, McGraw Hill, New
York.
Morse, P. M., 1948, Vibration and Sound, Macmillan, New York.
Munjal, M. L., Sreenath, A. V. and Narasimhan, M. V., 1973, An algebraic algorithm for
the design and analysis of linear dynamical systems, Journal of Sound and Vibration, 26,
pp. 193-208.
203
Munjal, M. L. and Prasad, M. G., 1986, On plane-wave propagation in a uniform pipe in
the presence of a mean flow and a temperature gradient, Journal of the Acoustical Society
of America, 80, pp. 1501-1506.
Nicholson, J. W. and Bergman, L. A., 1986, Free vibration of combined dynamical
systems, Journal of Engineering Mechanics, 112, pp. 1-13.
Peat, K. S., 1988, The transfer matrix of a uniform duct with a linear temperature gradient,
Journal of Sound and Vibration, 123, pp. 43-53.
Peat, K. S., 1997, Convected acoustic wave motion along a capillary duct with an axial
temperature gradient, Journal of Sound and Vibration, 203, pp. 855-866.
Perkins, N. C. and Mote, C. D. Jr., 1986, Comments on curve veering in eigenvalue
problems, Journal of Sound and Vibration, 106, pp. 451-463.
Polya, G., 1943, On the zeros of the derivatives of a function and its analytic character,
BULLETIN of the American Mathematical Society, 49, pp.178-191.
Pota, H. R. and Kelkar, A. G., 2000, Analysis of perfect noise cancelling controllers,
Proceedings of the 2000 IEEE International Conference on Control Applications, Sept,
pp. 25-27.
Pouyet, J. M. and Lataillade, J. L., 1981, Tosional vibrations of a shaft with non-uniform
cross section, Journal of Sound and Vibration, 76, pp. 13-22.
Prestel, E. C. and Leckie, F. A., 1963, Matrix methods in elastomechanics, McGraw-Hill,
New York.
Raspet, R., Bass, H. E. and Kordomenos, J., 1993, Thermoacoustics of traveling waves:
Theoretical analysis for an inviscid ideal gas, Journal of the Acoustical Society of
America, 94, pp. 2232-2239.
Rayleigh, J. W. S., 1945, The Theory of Sound, Dover, New York.
Ross, C. F., 1981, A demonstration of active control of broadband sound, Journal of
Sound and Vibration, 74(3), pp. 411-417.
Snowdon, J. C., 1968, Vibration and Shock in Damped Mechanical Systems, Wiley, New
York.
Subrahmanyam, P. B., Sujith, R. I. and Lieuwen, T. C., 2001, A family of exact transient
solution for acoustic wave propagation in inhomogeneous, non-uniform area ducts,
Journal of Sound and Vibration, 240, pp. 705-715.
204
Subrahmanyam, P. B., Sujith, R. I. and Lieuwen, T. C., 2003, Propagation of sound in
inhomogeneous media: Exact, transient solutions in curvilinear geometries, Journal of
Vibration and Acoustics, 125, 133-136.
Sujith, R. I., 1996, Transfer matrix of a uniform duct with an axial temperature gradient,
Journal of the Acoustical Society of America, 100, pp. 2540-2542.
Sujith, R. I., 2001, Exact solutions for modeling sound propagation through a combustion
zone, Journal of the Acoustical Society of America, 110, pp. 1839-1844.
Sujith, R. I., Waldherr, G. A. and Zinn, B. T., 1995, An exact solution for one-
dimensional acoustic fields in ducts with an axial temperature gradient, Journal of Sound
and Vibration, 184, pp. 389-402.
Swinbanks, M. A., 1973, The active control of sound propagation in long ducts, Journal
of Sound and Vibration, 27(3), pp. 411-436.
Tan, C. A. and Kuang, W., 1995, Vibration of a rotating discontinuous shaft by the
distributed transfer function method, Journal of Sound and Vibration, 183(3), pp. 451-
474.
Tang, Y., 2003, Numerical evaluation of uniform beam modes, Journal of Engineering
Mechanics, 129, pp. 1475-1477.
Timoshenko, S., 1981, Theory of Vibration with Applications, Prentice-Hall, Englewood
Cliffs, New Jersey.
Trindler, M. C. J. and Nelson, P. A., 1983, Active noise control in finite length ducts,
Journal of Sound and Vibration, 89, pp. 95-105.
Wei, G. W., Zhao, Y. B. and Xiang, Y., 2002, A novel approach for the analysis of high-
frequency vibrations, Journal of Sound and Vibration, 257, pp. 207-246.
Wittrick W.H. and Williams F.W., 1971, A general algorithm for computing natural
frequencies of elastic structures, Quarterly Journal of Mechanics and Applied
Mathematics, 8, pp. 263-284.
Wu, S. R., 2005, A classical solutions of forced vibration of rod and beam driven by
displacement boundary conditions, Journal of Sound and Vibration, 279, pp. 481-486.
Wu, X., 1997, Closed-form transient analysis of one-dimensional distributed dynamic
systems, Ph.D. Dissertation, University of Southern California.
205
Yang, B. and Fang, H., 1994, Transfer function formulation of non-uniformly distributed
parameter systems, ASME Journal of Vibration and Acoustics, 116, pp. 426-432.
Yang, B. and Mote, C. D. Jr., 1990, Vibration of band saws: theory and experiment,
Wood Science and Technology, 25, pp. 355-373.
Yang, B. and Mote, C. D. Jr., 1991, Controllability and observability of distributed
gyroscopic systems, ASME Journal of Dynamic Systems, Measurement, and Control,
113, pp. 11-17.
Yang, B. and Mote, C. D., Jr., 1991, Active vibration control of the axially moving string
in the s-domain, ASME Journal of Applied Mechanics, 58, pp. 189-196.
Yang, B. and Mote, C. D., Jr., 1992, On time delay in non-collocated control of flexible
mechanical systems, ASME Journal of Dynamic Systems, Measurement, and Control,
114, pp. 409-415.
Yang, B. and Tan, C.A., 1992, Transfer functions of one-dimensional distributed
parameter systems, ASME Journal of Applied Mechanic, 59, pp. 1009-1014.
Yang B., 1994, Distributed transfer function analysis of complex distributed parameter
systems, ASME Journal of Applied Mechanics, 61, pp. 84-92.
Yang, B., 1995, Linear vibration of a coupled string-rigid bar system, Journal of Sound
and Vibration, 183, pp. 383-399.
Yang, B., 1996, Closed-form transient response of distributed damped systems, part ii:
energy formulation for constrained and combined system, ASME Journal of Applied
Mechanics, 63, pp.1004-1010.
Yang, B., 2005, Stress, Strain, and Structural Dynamics: An Interactive Handbook of
Formulas, Solutions, and MATLAB Toolboxes, Elsevier Science, Boston.
Yang, B., 2008, A distributed transfer function method for heat conduction problems in
multilayer composite, Numerical Heat Transfer, Part B Fundamentals, 54, pp. 314-337.
Yang, B., 2010, Exact transient vibration of stepped bars, shafts and strings carrying
lumped masses, Journal of Sound and Vibration, 329, pp. 1191-1207.
Yang, K., 2008, A unified solution for longitudinal wave propagation in an elastic rod,
Journal of Sound and Vibration, 314, pp. 307-329.
206
Zu, J. W. Z. and Han, R. P. S., 1992, Natural frequencies and normal modes of a spinning
Timoshenko beam with general boundary conditions, ASME Journal of Applied
Mechanics, 59, pp. S197-S204
Zu, J. W. Z. and Han, R. P. S., 1992, Natural frequencies and normal modes of a spinning
Timoshenko beam with general boundary conditions, Discussion/Authors closure, ASME
Journal of Applied Mechanics, 59, pp. 1046-1047
Zu, J. W. Z. and Han, R. P. S., 1994, Dynamic response of a spinning Timoshenko beam
with general boundary conditions and subjected to a moving load, ASME Journal of
Applied Mechanics, 61, pp. 152-160.
207
Appendix A: Derivation of solutions in Section 3.5
The solutions given in Section 3.5 are obtained through introduction of proper
coordinate transformations, which convert the original (3.20) into differential equations
with known solutions.
We first consider that
2
( ) ( ) A x x , where function ( ) 0 x for 0 xL . Let the
solution of (3.20) be of the form
( , ) ( ) ( , ) z x s x u x s (A.1)
By which, (3.20) is converted to
2
2
2
0
z
z
x
(A.2)
where
22
d dx . Based on (A.2), solutions in Types 4 to 5 of Section 3.5 can be
obtained.
If ( ) 0 x , () x ax b , and (A.2) has two independent solutions
x
e
,
x
e
, which
eventually leads to the result for Type 4. If
2
/0 a which a being a constant,
()
ax
xe or
ax
e
, and (A.2) becomes
2
22
2
0
z
az
x
(A.3)
Solution of (A.3) yields the result for Type 3. Note that () x can also be cosh( ) ax b ,
which is a linear combination of
ax
e and
ax
e
. This leads to the two independent
208
solutions for Type 6. If
2
/0 a
, ( ) cos( ) x ax b , and (A.2) yields the solutions
of Type 5.
For Type 7, by the coordinate transformation similar to that in (Li, 2000)
1
2
( , ) (1 ) ( , )
( 1/ )
b
u x s ax z s
xa
(A.4)
Eq. (3.20) is converted to the modified Bessel equation, with order number (1 ) / 2 nb .
Thus, two independent solutions of (3.20) are obtained in Eq. (3.29). Type 2 is a special
case of Type 7, with b = 1.
209
Appendix B: Proof of Eq. (4.36)
Modified Bessel functions have the following relations with Bessel functions
( ) ( )
n
nn
I j j J
and
1
( ) ( ) ( )
2
n
n n n
K j j J jY
(B.1)
and the asymptotic relations for large values of z is expressed by (Mei, 1997)
22
( ) ~ cos (2 1) , ( ) ~ sin (2 1)
44
nn
J z z n Y z z n
zz
(B.2)
These formulas further confirm that ()
n
Jz and ()
n
Yz are the mathematical counterparts of
cosx and sinx, respectively.
The state transition matrix of a duct with a linear temperature profile is computed from
the fundamental matrix in Section 4.5. Using the relations (B.1) and letting sj in Eq.
(4.16), the characteristic equation can be rewritten as
0 2 0 1 0 2 0 1
0 J T Y T Y T J T
a a a a
(B.3)
where
2
m
aR .
By Eq. (B.2), Eq. (B.3) yields
12
4
21
21
sin 0
a
TT
a TT
(B.4)
The asymptotic eigenfrequencies of a duct are roots of (B.4) are
12
/,
k
ak T T
1,2, k . By Eq. (4.15), the eigenfunctions of a duct have of a form
210
0 0 1 0 0 1
()
k k k k
kk
v x A J T Y T Y T J T
a a a a
(B.5)
Assume that a duct has a linear temperature profile,
0
() T x T mx and
0
/1 mT . By
Eq. (B.2), Eq. (B.5) yields
1
1
2
( ) ~ sin
k
kk
k
a
v x A T T
a T
(B.6)
By Eq. (B.6), it is easy to show that
*
1
*
2 1 1
*
1
RHS of (4.36) ~ sin
cos sin
sin ~ LHS of (4.36)
dk
dk
T k
ks
T kk
ks
k
ks
A T T e
a
A T T T T e
aa
A T T
a
(B.7)
where
*
1
2
kk
k
a
AA
T
.
211
Appendix C: Proof of Eq. (4.37)
The state transition matrix of a duct with a conical temperature profile is computed
from the fundamental matrix in Section 4.5. By Eq. (4.16) and the given state transition
matrix, the characteristic equation of a duct with a conical temperature profile is given by
1 2 2 1
2 1 2 1
0 T T T T
(C.1)
where
2
2
10
a
,
12
and are the roots of the equation, i.e.
2 2 2
/0 m , and
2
1,2 2
11
1 with
2 2 2
m
aR
a
(C.2)
The eigenfrequencies
k
of a duct are roots of Eq. (C.1). If
2
2
10
a
, the roots
1,2
are
real and there is no eigenfrequencies satisfying Eq. (C.1).
If
2
2
10
a
, the roots
1,2
are complex. When 0 z and the exponent c is any
complex number, the function
c
z is defined by
ln c c z
ze (Brown and Churchill, 2003). Eq.
(C.1) can be rewritten as
2
1
1
2
12 2
2
2 sin 1ln 0
T
j TT
aT
(C.3)
The eigenfrequencies
k
of a duct are non-zero positive roots of Eq. (C.3) which are
given by
212
2
12
1 2 1 2
ln /
1 ~ , 1,2,
ln / ln /
k
TT
ak ak
k
k T T T T
(C.4)
A duct has a conical temperature profile (Type 4). Assume that
0
/1 mT . By Eq.
(4.15), the eigenfunction of a duct have of a form
1
11
( ) ~ sin (ln ln )
k k k
v x A T q T T
(C.5)
where
2
2
1
k
k
q
a
.
By Eq. (C.5), it is easy to show that
1
11
1
1 1 2 1
1
1 1 2 1
1
11
RHS of (4.37) ~ sin ln ln
sin (ln ln ) (ln ln )
cos ln ln sin (ln ln )
sin (ln ln ) ~ LHS of (4.37)
dk
dk
dk
T
k k s
T
k k a
T
k k k a
k k a
A T q T T e
A T q T T T T e
A T q T T q T T e
A T q T T
(C.6)
213
Appendix D: Proof of Eq. (4.39)
By Eq. (4.15), the eigenfunctions of a duct have of a form
0 1 1 0 1 1
()
k k k k
kk
v x A Y T J T J T Y T
a a a a
(D.1)
Using Eq. (B.2), Eq. (D.1) yields
1
1
2
( ) ~ cos
k
kk
k
a
v x A T T
a T
(D.2)
By Eq. (D.2), it is easy to show that
3 2 1
3 1 2
21
*
( ) ( ) ~ ( ) ( )
~ sin
k k a k k
j T j T j T j T
k a k s k s k s
k
k s s
v x e v x e v x e v x e
jA T T
a
(D.3)
where
*
1
2
kk
k
a
AA
T
and
3 12
1 11 1
1 2 3
, , ,
s ss a
a
TT T T T T TT
T T T T
a a a a
(D.4)
Rearranging Eq. (D.3) and using Eq. (D.4) leads to Eq. (4.39).
214
Appendix E: Proof of Eq. (4.40)
By Eq. (4.15), the eigenfunctions of a duct have of a form
11
1
( ) sin (ln ln ) cos (ln ln )
2
k k k k k
v x A q T T q q T T
(E.1)
where
2
2
1
k
k
q
a
.
By Eq. (E.1), it is easy to show that
3 2 1
3 1 2
12
( ) ( ) ~ ( ) ( )
1
~ sin ln ln
2
k k a k k
j T j T j T j T
k a k s k s k s
k k k s s
v x e v x e v x e v x e
A jq q T T
(E.2)
where
1 2 3
1 1 1 1
1 2 3
ln , ln , ln , ln
k k k k
a
k s k s k s k a
q q q q T T T T
T T T T
T T T T
(E.3)
By asymptotic formula, Eq. (E.3) can be approximated as
3 12
1 11
1
1 2 3
ln( / ) ln( / ) ln( / )
ln( / )
, , and
s ss
a
a
TT T T T T
TT
T T T T
a a a a
The accuracy increases as the mode goes higher.
215
Appendix F: Derivation of Eq. (4.62)
With the fundamental matrix in Section 4.5, the state transition matrix is obtained by
1
( , , ) ( , ) ( , ) x s x s s
Φ UU . For convenience of analysis, the state transition matrix is
defined as
12
34
( , , ) ( , , )
( , , )
( , , ) ( , , )
x s x s
xs
x s x s
Φ (F.1)
For a duct with a linear temperature distribution
0
() T x T mx , the open-loop transfer
function ˆ ( , , )
o
g x s is the 1-row and 2-column element of the distributed transfer function
Ĝ, and is given by
21
12
2
22
2
( ,0, ) ( ,0, )
( ,0, ) (0, , ),
( ,0, )
ˆ ( , , )
( ,0, ) ( , , )
,
( ,0, )
o
x s L s
x s s x
Ls
g x s
x s L s
x
Ls
(F.2)
where
i
is the element of the state transition matrix.
Let the open-loop poles be
k
j . By Eqs. (B.1), (B.2) and (4.42), it is found that
3/4
12
1/4
0
21
21
3/4
2 2
1/4
21
sin sin
()
()
cos
()
sin
()
kd
jT
kk
sa
k a k
k
kk
s
a ak
k
kk
s
T T T T e
ds aT N j j
aa
d D j
T T T
TT
a
TT aT N j j
D j a
T T T
(F.4)
where
12
d
TT
T
a
.
Abstract (if available)
Abstract
An analytical method, based on the distributed transfer function method (DTFM), is generalized and developed for modeling, analysis and control of one-dimensional distributed parameter systems. Distributed parameter systems encompass a variety of applications in wide engineering fields. Aerospace structures, building, bridges, turbines, rotating machines and active sound suppression are examples of the application of structural dynamics and control. ❧ Distributed parameter systems studied in this work are classified as a single-body distributed systems and a multi-body (complex) distributed system. A single-body distributed system contains one continuum with boundary constraints
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Static and dynamic analysis of two-dimensional elastic continua by combined finite difference with distributed transfer function method
PDF
Medium and high-frequency vibration analysis of flexible distributed parameter systems
PDF
Modeling, analysis and experimental validation of flexible rotor systems with water-lubricated rubber bearings
PDF
An approach to dynamic modeling of percussive mechanisms
PDF
Control of spacecraft with flexible structures using pulse-modulated thrusters
PDF
Dynamic modeling and simulations of rigid-flexible coupled systems using quaternion dynamics
PDF
Optimal design, nonlinear analysis and shape control of deployable mesh reflectors
PDF
Mathematical model and numerical analysis of a shape memory polymer composite cantilever beam for space applications
PDF
Control of two-wheel mobile platform with application to power wheelchairs
PDF
Modeling and dynamic analysis of coupled structure-moving subsystem problem
PDF
Modeling and vibration analysis of wheelchair-users
PDF
On the synthesis of dynamics and control for complex multibody systems
PDF
An approximation of frequency response of vibration system in mid-frequency area
PDF
Transient modeling, dynamic analysis, and feedback control of the Inductrack Maglev system
PDF
Periodic solutions of flexible systems under discontinuous control
PDF
Nonlinear control of flexible rotating system with varying velocity
PDF
Dynamic social structuring in cellular self-organizing systems
PDF
Design, modeling and analysis of piezoelectric forceps actuator
PDF
Active delay output feedback control for high performance flexible servo systems
PDF
Robust control of periodically time-varying systems
Asset Metadata
Creator
Noh, Kyoungrae
(author)
Core Title
Dynamic analysis and control of one-dimensional distributed parameter systems
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Publication Date
02/01/2013
Defense Date
06/12/2012
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
Control,distributed parameter systems,distributed transfer function method,dynamic analysis,OAI-PMH Harvest
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Yang, Bingen (Ben) (
committee chair
), Flashner, Henryk (
committee member
), Shiflett, Geoffrey R. (
committee member
), Wellford, L. Carter (
committee member
)
Creator Email
knoh@usc.edu,kyoungrae.noh@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-82178
Unique identifier
UC11289883
Identifier
usctheses-c3-82178 (legacy record id)
Legacy Identifier
etd-NohKyoungr-1106.pdf
Dmrecord
82178
Document Type
Dissertation
Rights
Noh, Kyoungrae
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
distributed parameter systems
distributed transfer function method
dynamic analysis