+D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 157

+D=FJAH

Dynamics of Multi-Degree-
of-Freedom Systems

11.1 INTRODUCTION
A Multi-Degree-of-Freedom (MDOF) system, as the name suggests, is one that requires two or
more independent coordinates to describe its motion. The coordinates normally used to describe
the motion of a structural system, may be related to each other via some constraints, which could
either be simple kinematic relations between various coordinates, or they could arise from the
consideration of equilibrium of forces. The number of generalised (independent) coordinates is
given by the difference between the total coordinates describing the motion of a system and the
number of constraint relations. For example, consider
the case of a double pendulum, which is constrained to Y
move in XY plane as shown in Figure 11.1. In Cartesian
l1
coordinate system the positions of two masses m 1 and m2 q1
are described by two pairs of Cartesian coordinates m1(x1, y1)
(x1, y1) and (x2, y2) respectively. These four coordinates
l
viz. x1, y1, x2, y2, however, are related to each other q2 2
through two constraint relations: m2(x2, y2)
X
x21 + y21 = l21, (11.1)
and (x2 – x1)2 + (y2 – y1)2 = l 22 (11.2) FIGURE 11.1 A double pendulum.

Thus the number of degrees of freedom (or, generalised coordinates) of the structural system
is 4 – 2 = 2. The angles q1 and q2 can be taken as the two independent generalised coordinates
to describe the motion of the masses. In certain systems it is possible to eliminate dependent
variables by using constraint relations and derive a set of generalised coordinates, are said
to possess holonomic constraints. On the other hand, there may exist some constraints, called
nonholonomic for which it is not possible to derive a set of independent coordinates.
Nonholonomic constraints are rarely encountered in practice, so it will be assumed in the
following that the equations of the dynamic equilibrium of the system are specified in the
unconstrained coordinate system.
157
 158 Earthquake Resistant Design of Structures

11.2 SYSTEM PROPERTY MATRICES

As mentioned earlier, every dynamical system comprises
(i) a mechanism for storing strain energy due to deformations,
(ii) some means of storing kinetic energy of the system in motion, and
(iii) an energy dissipation mechanism.
In MDOF system, described by a set of N generalised coordinates (say, v1, v2, …, vN), these
energy functionals depend on the motion of the system described by the generalised coordinates.
N
1 ∂ 2U
Potential energy = U(v1, v2, … vN, t) =
2 Â vv
∂vi ∂v j i j
i , j =1

N 2
Kinetic energy = T ( v
1 , v
2 , º, v
N , t ) =
1
2 Â ∂∂v
∂Tv
i j
v
i v
j
i , j =1
N 2
Rayleigh dissipation function = R ( v
1 , v
2 , º, v
N , t ) =
1
2 Â ∂∂v
∂Rv
i j
v
i v
j
i , j =1

where Rayleigh dissipation function represents the energy loss through velocity proportional
viscous damping force.
∂2 R
The element cij of damping matrix (C) is given by the coefficient and represents
∂v
i ∂v
j
the damping force at ith DOF corresponding to the unit velocity at j th DOF with the velocities
∂ 2U
at all other DOFs remaining zero. The coefficient is the element kij of stiffness matrix
∂vi ∂v j
(K) and represents restoring force at the ith DOF corresponding to the unit displacement at jth
DOF with displacements at all other DOFs being constrained to zero. Similarly, the coefficient
∂2T
is the element mij of inertia matrix (M) and represents the inertia force at the ith DOF
∂v
i ∂v
j
corresponding to the unit acceleration at the jth DOF with accelerations at all other DOFs being
constrained to zero. The governing differential equation of motion for an MDOF system can
be derived from the same principles as used in the case of SDOF systems. As an illustration,
let us consider the 2-DOF system as shown in Figure 11.2. By invoking the d’Alembert’s
principle, the equations of motion for free vibration of this system can be written as,
m v

1 + 2c v
1 – c v
2 + 3kv1 – kv2 = 0
m


v2 – c v
1 + 2c v
2 – kv1 + 3kv2 = 0
which can be arranged in the matrix form as,
LMm 0 OP FG v

IJ + LM 2c
1 -c OP FG v
IJ + LM3k
1 -k OP FG v IJ = FG 0IJ ,
1
or M


v + C v
+ Kv = 0 (11.3)
N 0 mQ H v

K N- c
2 2 c Q H v
K N - k
2 3k Q H v K H 0K
2
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 159
v1 v2
2k k 2k
m m

c c c

2kv1 .. k(v2 – v1) .. 2kv2

. mv1 . . mv2 .
cv1 c(v2 – v1) cv2

Free body diagram

FIGURE 11.2 A 2-DOF system.

The nature of damping forces is assumed to be of the viscous type primarily as an approximate
representation of the combined action of all energy dissipation mechanisms present in a vibrating
system. Since the extent of damping in structural systems is usually very small, precise nature
of the damping force is not very important for dynamic response computations.

F0 sin wt
11.3 DYNAMICS OF TWO DEGREE
OF FREEDOM SYSTEMS m1 m2
k1 k2
Let us consider the response of a harmonically excited 2-DOF
FIGURE 11.3 Harmonic exci-
(undamped) system as shown in Figure 11.3. The governing
tation of a 2-DOF system.
equation of motion for this system can be given as:

LMm 1 OP FG v

IJ + LMk + k

0 1 1 2 - k2 OP FG v IJ = FG F sinwtIJ
1 0
(11.4)
N0 m Q H v

K N - k
2 2 2 k2 Q Hv K H 0 K
2

This system is the characteristic of an industrial building with a reciprocating machine installed
at one of the floors. Since such machines typically operate at a fixed speed, the force exerted
by these machines on the building floor will be harmonic and the steady-state response of the
system to this harmonic excitation will also be harmonic of the same frequency.
Thus assuming the harmonic response as v = [v1, v2 ]T = sin wt[V1, V2]T, Equation (11.4)
can be solved for the response amplitudes V1 and V2 of the two masses as:

LMk1 + k 2 - w 2 m1 - k2 OP FG V IJ sin wt = FG F sinwtIJ

1 0

N - k2 k 2 - w 2 m2 Q HV K
2 H 0 K
or,
FG V IJ = 1 LMk
1 2 - w 2 m2 k 2 OP FG F IJ 0
(11.5)
HV K D N
2 k2 k + k -w m QH 0 K
2 2
2
1

where D = mm 21 [w4 – w 2 {w *1 + w*2 (1 + m)} + (w *1w *2 )2 ], w *1 = k1 / m1 , w *2 = k2 /m2 , and

 160 Earthquake Resistant Design of Structures

m = m2 /m1. Thus the system response can be given by,

F0 ( k2 - w 2 mm1)
v1(t) = sin wt
mm12 [w 4 - w 2 {w *1 + w *2 (1 + m )} + (w *1 w *2 ) 2 ]

F0 k2
v2(t) = sin wt (11.6)
mm12 [w 4 -w 2
{w *1 + w *2 (1 + m )} + (w *1 w *2 )2 ]

It may be noted from Equations 11.6 that it is possible to force the amplitude of response of
the first mass (m1) to vanish by a suitable choice of parameters (also refered to as tuning)
k2 and m2 (or, m). This concept can be exploited in designing vibration absorbers for industrial
structures and can be achieved by attaching an auxiliary/secondary mass to the primary structure
which is subjected to a harmonic excitation. This can be quite effective, when the operating
frequency of the reciprocating machine and a natural frequency of the supporting structure are
nearly equal, causing large amplitude vibrations of the supporting structure due to resonance.
Let us assume for simplicity, that the second mass m2 in Figure 11.3 is an auxiliary mass
attached to the primary structure, which is excited by a reciprocating machine installed on the
floor. The amplitude of displacement response (normalized with respect to the static deflection
F0 /k1) of the primary mass is shown in Figure 11.4(a) for a range of operating frequencies. The
unbounded amplitude for w/w *1 = 1.0, corresponds to the condition of resonance and it should
be avoided (w *1 = k1 / m1 represents the natural frequency of primary structure alone). Let us
now consider the use of a tuned mass damper (or, vibration absorber) in altering the dynamic
response of primary structure. Figures 11.4(b-f) show the response characteristics of primary
structure with an auxiliary/secondary structure for different parametric variations. The objective
is to find a suitable set of parameters for secondary structure, so as to limit the normalized
response of primary structure to be less than unity in the neighbourhood of w/w *1 = 1. The effect
of adding a secondary structure is to split one resonant peak into two resonant peaks,
corresponding to two natural frequencies of the combined 2-DOF system. A light (m < 1)
secondary system attached to the primary structure by means of a relatively flexible attachment,
does take away a significant part of the vibration energy from the primary structure at w/w *1
= 1, as can be seen in Figure 11.4(b). However, such an arrangement is only effective (in
reducing the vibration amplitude of primary mass) for a very small range of operating
frequencies, by means of increasing the deformations in the secondary structure. Attaching a
heavy secondary mass with a flexible spring is not at all effective, as can be infered from Figure
11.4(c). Figure 11.4(d) shows that the use of a light mass and relatively stiff attachment does
not lead to any significant change in the resonant frequency of the original structure and hence
is not an effective solution for the vibration problem. Figures 11.4(e) and (f) show the perform-
ance for secondary structures comprising heavy mass with stiff attachments. By comparing these
two figures, it may be noted that the two resonant peaks are well separated and are sufficiently
away from the original location of resonant peak (for the primary structure alone) at w/w *1 =
1. Further, the amplitude of displacement of primary structure vanishes at w/w *1 = 1 and is less
than the static deformation (F0 /k1) for a wide range of operating frequencies. This effective
range of operating frequencies increases with increase in the mass ratio (m = m2/m1).
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 161

(a) Without auxiliary mass (b) k2/k1 = 0.5, m2/m1 = 0.5

4.0 4.0
Primary mass Primary mass
3.0 3.0
Normalized displacement

Normalized displacement
Secondary mass
2.0 2.0
1.0 1.0
0.0 0.0
–1.0 –1.0
–2.0 –2.0
–3.0 –3.0
–4.0 –4.0
0.0 1.0 2.0 3.0 4.0 5.0 0.0 1.0 2.0 3.0 4.0 5.0
w/w *1 w/w *1
(c) k2/k1 = 0.5, m2/m1 = 2.0 (d) k2/k1 = 2.0, m2/m1 = 0.5
4.0 4.0 Primary mass
Primary mass
3.0 3.0 Secondary
Normalized displacement

Normalized displacement
Secondary mass
mass
2.0 2.0
1.0 1.0
0.0 0.0
–1.0 –1.0
–2.0 –2.0
–3.0 –3.0
–4.0 –4.0
0.0 1.0 2.0 3.0 4.0 5.0 0.0 1.0 2.0 3.0 4.0 5.0
w/w *1 w/w *1
(e) k2/k1 = 2.0, m2/m1 = 2.0 (f) k2/k1 = 4.0, m2/m1 = 4.0
4.0 4.0
Primary mass Primary mass
3.0 3.0
Normalized displacement

Normalized displacement

Secondary mass Secondary

mass
2.0 2.0
1.0 1.0
0.0 0.0
–1.0 –1.0
–2.0 –2.0
–3.0 –3.0
–4.0 –4.0
0.0 1.0 2.0 3.0 4.0 5.0 0.0 1.0 2.0 3.0 4.0 5.0
w/w *1 w/w *1
FIGURE 11.4 Performance of vibration absorber/tuned mass damper.

Thus a secondary structure with heavy absorber mass can be very effective in controlling
machine induced vibrations in a building/structure. However, the large size associated with a
heavy mass can sometimes impose a practical limitation on the usable range of operating
frequencies. The required stiffness (k 2) of the secondary attachement can be calculated from
the maximum allowable displacement of the secondary system (at w = w *2 ) V2,max =
F0 /k 2.
 162 Earthquake Resistant Design of Structures

11.4 FREE VIBRATION ANALYSIS OF MDOF

SYSTEMS
By considering the fact that the damping levels are usually very small in structural systems, let
us consider the response of an undamped MDOF system. The effect of damping will be dealt
with at a later stage. The equation of free vibration then reduces to,
M


v + Kv = 0 (11.7)
We look for a solution in the form vi = q(t)fi, i = 1, 2, …, N, where the dependence on time
and that on space variables can be separated. This implies that the ratio of amplitudes of any
two coordinates is independent of time. Physically, it implies that all degrees of freedom
perform synchronous motion and the system configuration does not change its shape during
motion but only its amplitude changes.
Substituting for v, the equation of motion may be written as,
M{f} q

(t) + K{f}q(t) = 0 (11.8)
which is a set of N simultaneous equations of the type
N N

 mij f j q

(t ) + Â k f q(t ) = 0; i = 1, 2, …, N

ij j (11.9)
j =1 j =1

where the separation of variables leads to

k f
N

-
q

(t )
=
 j =1 ij j
; i = 1, 2, …, N (11.10)
m f
q(t ) N
 j =1 ij j

Since the terms on either side of the equality sign are independent of each other, this equality
can hold only when each of these terms are equal to a positive constant (say, w 2).1 Thus we have,
q

(t) + w2q(t) = 0 (11.11)
N

 (k ij - w 2 mij )f j = 0; i = 1, 2, …, N. (11.12)
j =1

The solution of Equation (11.11) is q(t) = sin(wt – a), a harmonic of frequency w. Thus we
may conclude that the motion of all coordinates is harmonic with same frequency w and same
phase difference a. However, it still needs to be established if such a synchronous, harmonic
motion is possible for all frequencies. To investigate this issue let us consider the Equation

1
The choice of the sign is dictated by physical considerations. For a conservative system the displacements must
remain finite at all instances. If we had chosen a negative constant then the solution would involve exponential
functions which would grow without bounds with time t. The choice of positive constant, on the other hand,
provides a harmonic solution which has finite energy at all times.
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 163

(11.12), which is a set of N simultaneous linear homogeneous equations in unknowns fj . The

problem of determining constant (w2) for which the Equation (11.12) has a non-trivial solution
is known as the characteristic value or eigenvalue problem. The eigenvalue problem may be
rewritten, in matrix notation as,
(K – w 2 M){f} = 0 (11.13)
A non-trivial solution for Equation (11.13) is feasible only when the determinant of the
coefficient matrix vanishes, i.e.,
|K – w 2M| = 0 (11.14)
The expansion of the determinant in Equation (11.14) yields an algebraic equation of Nth order
in w 2, which is known as the characteristic equation. The roots of characteristic equation are
known as the eigenvalues and the positive square root of these eigenvalues are known as the
natural frequencies (wi) of the MDOF system. It is only at these N frequencies that the system
admits synchronous motion at all coordinates. For stable structural systems with symmetric and
positive definite stiffness and mass matrices the eigenvalues will always be real and positive. For
each eigenvalue the resulting synchronous motion has a distinct shape and is known as natural/
normal mode shape or eigenvector. The normal modes are as much a characteristic of the system
as the eigenvalues are. They depend on the inertia and stiffness, as reflected by the coefficients
mij and kij. These shapes correspond to those structural configurations, in which the inertia forces
imposed on the structure due to synchronous harmonic vibrations are exactly balanced by the
elastic restoring forces within the structural system. These eigenvectors are determined as the
non-trivial solution of Equation (11.13). Since the determinant of the coefficient matrix
evaluated at one of the natural frequencies is singular, a unique solution for eigenvectors can
not be found. It is, however, possible to compute the amplitudes of the synchronous motion at
N – 1 coordinates relative to the amplitude of motion at the remaining coordinate, which may
be selected arbitrarily. Thus an additional constraint—known as normalisation condition—must
be supplied in addition to Equation (11.13) to completely determine an eigenvector. Two of the
most commonly used normalisation procedures are:
(i) assume the amplitude of synchronous motion at the first degree of freedom as unity,
(ii) constrain a length measure of the eigenvector to be unity. For example, for any eigen-
vector {f (i)} it is possible to determine elements of {f (i)} such that {f i}TM{f (i)} = 1.
Such a normalisation, using mass/inertia matrix (M) is known as mass renormalisation
and the resulting mode shape is known as mass orthonormal mode shape.
It can be shown that the N eigenvectors of an N-DOF system completely span the N-dimensional
vector space, and therefore, can be used as basis vectors for representing any Nth order vector.
Since the condition of orthogonality is a necessary condition for any set of base vectors, it will
now be shown that the eigenvectors also satisfy this condition.

11.4.1 Orthogonality Conditions

An important property of the mode shapes or eigenvectors is that they are mutually orthogonal
with respect to the mass and stiffness matrices. More precisely, the product involving
 164 Earthquake Resistant Design of Structures

multiplication of mode shapes corresponding to two different modes vanishes.

{f (j)}T [M]{f (i)} = Midij and {f(j)}T [K]{f (i)} = Kidij (11.15)
where Mi and K i are called the generalised mass and stiffness respectively for the i th mode and
dij is the Kronecker delta. For the case when mode shapes have been orthonormalized with
respect to the mass, Mi = 1 and Ki = w 2i .
In order to prove this proposition, let us assume that w 2i and {f (i)} denote the eigenvalue
and corresponding eigenvector for i th mode and w 2j and {f ( j)} correspond to the j(π i)th mode.
It follows that both of these eigenpairs satisfy Equation (11.13). Thus,
K{f(i)} = w 2i M{f(i)} (11.16)
K{f } =(j)
w 2j M{f (j)} (11.17)
Pre-multiplying Equation (11.16) by {f(j)}T and Equation (11.17) by {f (i)}T, we get,
{f(j)}TK{f(i)} = w 2i {f(j)}TM{f(i)} (11.18)
(i) T
{f } K{f } = (j)
w 2j (i) T
{f } M{f } (j)
(11.19)
Subtracting Equation (11.19) from Equation (11.18) and noting {f(j)}TK{f(i)} = ({f(i)}TKT{f(j)})T
and the fact that K and M are symmetric matrices we have,
(w 2i – w 2j ){f(j)}TM{f(i)} = 0 (11.20)
For all modes i π j with distinct eigenvalues (w i π wj), Equation (11.20) can be specified only
if the matrix inner product {f(j)}TM{f(i)} vanishes. This proves the first half of the proposition
stated in Equation (11.15). The other half follows by substituting this result in either of
Equation (11.18) or (11.19).
Since the computed mode shapes of a N-DOF system form a set of orthogonal vectors, they
span the N-dimensional space completely. In other words, these mode shapes can be used as a
set of basis vectors in the N-dimensional space and any vector in this space can be represented
as a linear combination of these mode-shapes.
The orthogonality property of mode shapes leads to a very powerful theorem, modal
expansion theorem, which states that any vector x in N-dimensional vector space can be
represented as a linear combination of mode-shape vectors,
N
x= Â q {f i
(i )
} (11.21)
i =1

where {f(i)} represents the ith mode shape and qi denotes the corresponding modal coordinate.
For a given vector x, the modal coordinates qi may be computed by using the property of
orthogonality of mode-shapes as,

{f ( i )}T Mx
qi = . (11.22)
{f (i )}T M{f (i )}
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 165

11.5 DETERMINATION OF FUNDAMENTAL

FREQUENCY
The determination of the eigenspectrum of a system is an important part of the dynamic analysis
of the system. Since the response of MDOF system is usually contained in the lower modes of
vibration, determination of the characteristics of the fundamental mode is of primary interest.

11.5.1 Rayleigh Quotient

For any arbitrary vector, {u}, representing a displacement configuration of a N-DOF system,
the Rayleigh quotient is defined as the ratio

u T Ku
r= (11.23)
u T Mu
For a particular case when the vector u represents the amplitudes of the harmonic oscillations
of the N-DOF system or the Rayleigh quotient, r, corresponds to square of the frequency of
harmonic oscillations. This result follows from the principle of conservation of energy by
equating the maximum potential energy stored in the system to the maximum kinetic energy.
Further, the Rayleigh quotient has the property of being stationary in the neighbourhood of the
natural modes of the system. It is a global minimum for the fundamental mode and global
maximum for the highest mode of vibration—also known as the minimax property of Rayleigh
quotient.

11.5.2 Stodola Method

By transforming the generalised eigenvalue problem to the standard eigenvalue problem,
D{f} = l{f} (11.24)

where D = K–1M is known as the dynamical matrix of the system and l =

1
. Stodola method
w2
starts with the choice of a trial vector, say, {f ( 0)} . Pre-multiplying {f ( 0)} by the dynamical
matrix, D yields another vector {f (1)}, which is an improved estimate of the eigenvector. An
estimate of the eigenvalue is obtained by taking the ratio of any element of new vector {f (1)}
to the corresponding element of the trial vector, i.e.,

f (j1)
l (1) = (11.25)
f ( 0 )
j

If {f (1)} were a true eigenvector, this ratio would be constant for any choice of the element of
these vectors. In general, however, this ratio will be different for different choice of elements
of these vectors. In the special case of symmetric coefficient matrices, the minimum and the
maximum values of this ratio provide the upper and lower bounds on the eigenvalue. The
 166 Earthquake Resistant Design of Structures

iteration resumes with the new trial vector of {f (1 )} = {f (1)} . Thus the equation for ith
1
l (1)
iteration is given as,
D{f (i )} = l (i +1) {f (i +1)} (11.26)
The Stodola method can be viewed as an iterative solution of a system of simultaneous equations
to arrive at that configuration of generalised displacements for which the inertia forces are
exactly balanced by the elastic forces in the structural members.

Why should iterative procedure converge to the first mode always?

To answer this natural query, let us take recourse to the modal expansion theorem and expand
an arbitrary trial vector {f} as,

{f} = q1{f (1)} + q2{f(2)} + ◊ ◊ ◊ + qN{f(N)} (11.27)

where {f(1)}, {f(2)}, …, {f(N)} denote the eigenvectors of the dynamical system. The first
iteration results in,
q1 (1) q q
D{f} = {f } + 22 {f ( 2 )} +  + N2 {f ( N )} (11.28)
w1
2
w2 wN

Thus each iteration results in amplification of the ith term in the modal expansion by a factor
1/w 2i . So that after p successive iterations,
q q q
Dp{f} = 21p {f (1)} + 22p {f (2 )} +  + 2Np {f ( N )} (11.29)
w1 w2 wN
Assuming that the natural frequencies (wi) are all distinct and are numbered in the ascend-
ing order i.e. w1 < w2 < ◊ ◊ ◊ < wN, it follows that after sufficient number of cycles
1 1 1
>> 2 p >>  >> 2 p . Therefore the first term in the modal expansion becomes progres-
w 12 p w2 wN
sively more dominant with each iteration and eventually converges to the first mode {f (1)}.

11.5.3 Converging to Higher Modes

The iteration method described earlier will always converge to the lowest mode, unless the
chosen trial vector exactly resembles a higher natural mode. Therefore to determine the higher
modes using iteration procedure, it is necessary to sweep out all the lower modes. For example,
let us assume that the first mode shape has already been determined and has been mass
orthonormalized (such that {f(1)}TM{f(1)} = 1.0). Considering any arbitrary trial vector {f} and
pre-multiplying it by {f(1)}TM and invoking the orthogonality of mode shapes,

{f(1)}TM{f} = {f(1)}TM{f(1)}q1 + {f(1)}TM{f(2)}q2 + ◊◊◊ + {f(1)}TM{f(N)}qN

= q1{f(1)}TM{f(1)} (11.30)
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 167

where, use has been made of the modal expansion theorem (Equation (11.27)). The coefficient
q1 so determined, gives the extent of representation of the mode-shape {f(1)} in the trial vector
{f} . Let us define a new trial vector {y } = {f} – q1{f(1)} by sweeping out the traces of known
eigen vector {f(1)}. Using this purified trial vector in the iteration procedure we would converge
to the next lowest mode i.e. 2nd mode, {f(2)}. This process can be repeated to compute any
desired eigenvector by sweeping out the traces of all the previous lower mode eigenvectors from
the trial vector. A geometrical interpretation of the process of sweeping is to determine a trial
vector which is orthogonal to all the previously determined eigenvectors and this approach is
known as vector purification/deflation.
In theory, though it is possible to sweep out completely the traces of a known eigenvector
from an assumed trial vector, in practice, however, it is necessary to sweep out the known
eigenvectors from trial function before the beginning of each iteration. This precaution is
necessary because the round off errors due to finite precision arithmetic, on a computer always
leave some small traces of swept out eigenvector(s) in the trial vector at the end of the iteration.
It is possible to automate the process of sweeping in each iteration by sweeping out the traces
of known modes from the coefficient matrix. Let us consider that first n(< N) modes are known
and it is required to converge to the n+1th mode via iteration. Since the need for sweeping the
traces of known modes from trial vector at each iteration may be computationally expensive,
it is worthwhile to look for the possibility of a more elegant formulation for this procedure. Let
us consider that {y ~} be the trial vector from which the traces of first n(< N) modes are to be
removed. We have, by modal expansion theorem,
N
{y
~} =
 {f ( j )} q j
j =1

For first n modes, which are known, the coefficients qj can be computed by using the orthogo-
nality property. The purified trial vector {y } can then be given as,

N
{y } = {y~} - Â {f ( j )} q j
j =1

F n
1 I
{f ( j )}{f ( j )}T M {y
~}
GH Â {f
= I-
j =1
} M{f }
( j) T ( j) JK (11.31)

= S{y
~}

where S is known as the sweeping matrix and the entire process of sweeping out the known
modes from a trial vector has been reduced to a simple matrix multiplication. In practice, the
coefficient matrix of the eigenvalue problem is post-multiplied by the sweeping matrix and the
resulting updated coefficient matrix is used in the iteration procedure to converge to the n + 1th
mode. The sweeping matrix is then updated to sweep out the first n + 1 modes by extending
the summation in Equation (11.31) to include the n + 1th mode.
 168 Earthquake Resistant Design of Structures

v1
For the demonstration of the procedure, let us consider m
a three-storey shear building shown in Figure 11.5. The
system parameters are given as m = 3500 kg, k1 = k = 1500 k1
v2
kN/m, k2 = 1.5k, and k3 = 2.0k. The mass and stiffness m
matrices can be written as,
k2
L1 0 0 O L1 -1 0 O m
v3
M = m M0 0P K = k M -1 -1.5P
MM0 1
P
1 PQ
MM 0 2.5
P
3.5PQ k3
N 0 N -1.5

For inverse iteration, the system coefficient matrix for the

standard form of eigenvalue problem is given as, FIGURE 11.5 A 3-storey shear
building.

m
LM
2.167 1167
. 0.5 OP
D = K–1M =
k
.
MM
1167 .
1167 0.5
PP
0.5N 0.5 0.5 Q
TABLE 11.1 Iteration for the first mode
k
D {y (0)} {y (1)} {y(1)} {y (2)} {y (2)} {y (3)} {y (3)} {y (4)} {y (4)}
m
2.167 1.167 0.500 1.00 2.875 1.00 3.075 1.00 3.109 1.00 3.12 1.00
1.167 1.167 0.500 0.50 1.875 0.65 2.075 0.67 2.109 0.68 2.12 0.68
0.500 0.500 0.500 0.25 0.875 0.30 0.975 0.32 0.995 0.32 1.00 0.32

Thus, the iteration procedure converges to the first eigenvalue l1 = 3.12m/k and the
corresponding eigenvector is {f(1)} = [1.00, 0.68, 0.32]T. Table 11.1 shows the iteration for the
first mode and accordingly, the natural frequency of the first mode of the structural system is
k
given by w 21 = 0.32 .
m

Sweeping
The sweeping matrix for removing the first mode is given by,

1
S1 = I – {f (1)}{f (1)}T M
{f (1)}T M{f (1)}

LM 0.361 -0.435 - 0.204 OP

= -0.435 0.705 - 0.139
MM-0.204 PP (11.32)
N -0.139 0.935 Q
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 169

Table 11.2 shows iteration for second mode and the coefficient matrix for the iteration for
second mode is given by,

m
LM
0.173 -0.189 -0.137 OP
-0.189 0.246 0.067
D1 = DS1 =
k MM PP
N
-0.137 0.067 0.296 Q
TABLE 11.2 Iterations for second mode
k
D {y(0)} {y (1)} {y(1)} {y ( 2)} {y(2)} {y (3)} {y(3)} {y ( 4)} {y(4)} {y ( 5)}
m 1
0.173 –0.189 –0.137 1.00 0.302 1.00 0.49 1.00 0.496 1.0 0.498 1.00 0.499
–0.189 0.246 0.067 –0.50 –0.329 –1.09 –0.51 –1.04 –0.506 –1.02 –0.505 –1.01 –0.503
–0.137 0.067 0.296 –0.25 –0.245 –0.81 –0.45 –0.92 –0.479 –0.966 –0.49 –0.986 –0.497
The approximation to eigenvector after 5th iteration is {y(5)}T = [1.000, –1.008, –0.996].

Thus, as the iterations proceed, the iteration procedure converges to the second eigenvalue
l2 = 0.5m/k and the corresponding eigenvector is {f(2)} = [1.00, –1.00, –1.00]T. Accordingly,
k
the natural frequency of the second mode of the structural system is given by w 22 = 2.0 . From
m
the elementary linear algebra, it is known that the trace of a square matrix is equal to the sum
of its eigenvalues. Thus, it follows that,
3
Tr(D) = Â lj
j =1

m k
and l3 = 0.214 , or w 23 = 4.68 . The corresponding eigenvector can be computed as {f(3)}
k m
= [1.0, –3.68, 4.68]T. Alternatively, the third eigenpair could have been computed by first
constructing a new sweeping matrix as
2
1
S2 = I Â {f ( j )}{f ( j )}T (11.33)
j =1
{f ( j )}T M{f ( j )}

and then deriving the coefficient matrix for iteration for the third mode (D2) by pre-multiplying
the new sweeping matrix by D, i.e.,
D2 = DS2
The remaining eigenpair may then be computed via iterations.

11.6 FORCED VIBRATION ANALYSIS

The forced vibrations of an MDOF system are described as a set of N coupled, non-homogeneous
differential equations in v as,
M


v + C v
+ Kv = f (11.34)1
 170 Earthquake Resistant Design of Structures

These equations, in coupled form, are extremely cumbersome, and we shall look for some
suitable transformation of the unknowns v to reduce the system of N coupled differential
equations to a set of N uncoupled differential equations. This method of solution is known as
mode-superposition method.

11.6.1 Mode-superposition Method

We know by modal expansion theorem that any arbitrary vector v in an N-dimensional space
can be represented as a linear combination of mode-shapes. Thus,
N
v= Â q (t) {f (r)} = F q
r (11.35)
r =1

where, F is the modal matrix with each of its columns representing a mode-shape of the MDOF
system and q is a vector of modal coordinates related to the system coordinate vector v through
a linear transformation given by Equation (11.35). Substituting this transformation in Equation
(11.34) and pre-multiplying it by F T, we have,
F TMF
Fq

+ FTCF
F q
+ FTKF
Fq = FT f (11.36)
Since the mode shapes F are orthogonal with respect to M and K matrices, the matrix triple
products involving M and K in Equation (11.36) both yield diagonal matrices. The damping
matrix C is not in general amenable to such diagonalization procedure. However, for a specific
class of damping matrices—called classical damping—such a diagonalization using (undamped)
mode shapes is indeed possible. A sufficient condition for a damping matrix C to be diagonalized
using undamped mode shapes is to have the following series expansion:

C=M Â a [M–1K]n
n (11.37)
n

A special case of Equation (11.37) obtained by retaining only two terms of the series for
n = 0 and n = 1 and is known as Rayleigh damping and is very widely used in structural dynamics
applications. It is also known as proportional damping as the damping matrix is proportional
to stiffness and mass matrices in this case. Since the matrices M and K are known for a given
structural system, a classical damping matrix C can be completely specified if the coefficients
an in the series of Equation (11.37) are specified. These coefficients can also determined so as
to have desired damping values in different modes of vibrations. If we assume that the C in
Equation (11.36) is a classical damping matrix, then the system of coupled equations reduces
to a set of N uncoupled differential equations in qr, r = 1, 2, …, N, as:
m r* q

r + c*r q
r + k*r qr = fr*, " r = 1, 2, …, N (11.38)
where, {f(r)}TM{f(r)} = m *r represents the modal mass for mode r, {f(r)}TC{f(r)} = cr* is the
coefficient of viscous damping in rth mode, {f(r)}TK{f(r)} = kr* denotes the modal stiffness for
rth mode, and {f(r)}T f = fr*, the modal force in mode r. It may be noted that if the mode shapes
have been mass-orthonormalized, then these modal parameters reduce to m *r = 1.0, cr* = 2z rwr,
and kr* = w 2r (note the similarity of form with the equation of motion for SDOF system in
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 171

Equation (7.5)). The Equations (11.38) can now be solved for unknowns qr, independent of each
other, by using the solution methods developed for single degree of freedom systems. Once the
solution for modal coordinates qr is available, the response in system coordinates can be obtained
by using the linear transformation of Equation (11.35).

11.6.2 Excitation by Support Motion

The forced vibration of MDOF system excited by support motions is described by the coupled
system of differential equations as
M


v + C v
+ Kv = –Mr


vg (11.39)
where v

g denotes ground acceleration, v is the vector of structural displacements relative
to the ground displacements, and r is a vector of influence coefficients. The ith element of vector
r represents the displacement of ith degree of freedom due to a unit displacement of the base.
The nature of this equation is similar to that of standard forced vibration problem as given by
Equation (11.34) and hence the method of solution (using mode-superposition) is also similar.
Thus the equation can be decoupled as
q

r + 2z rwr q
r + w 2r qr = –Gr v

g , " r = 1, 2, …, N (11.40)

{f (r )}T Mr
where, Gr = is known as the mode-participation factor for the rth mode.
{f ( r )}T M{f ( r )}
Note, however, that the Equation (11.40) differs from the equation of motion of a SDOF
system excited by support acceleration v

g by a scaling factor for the excitation. Since the
maximum response of SDOF system to ground acceleration is generally available in the form
of response spectra, it follows that the maximum value of the rth modal coordinate qr can be
determined directly from the response spectra without solving the differential equation of motion.
Therefore, assuming that the spectral displacement ordinate for frequency wr and damping zr
is given as Sd (wr, z r), the maximum response for rth modal coordinate qr,max is given as,
qr, max = Gr Sd(wr, zr); " r = 1, 2, …, N
This information about the maximum response in modal coordinates is, however, not very useful
for structural design, which is concerned with the maximum response in physical coordinates
v. It is possible to estimate probable maximum response values in physical coordinates from the
knowledge of maximum response in modal coordinates by using modal combination rules. Two
of the most commonly used modal combination rules are:

Absolute sum method

Assuming that the maximum of each modal coordinate occurs at the same instant of time, the
maximum response in physical coordinates at ith DOF is given by,
N

Âq r ,max f i
(r)
vi, max ª (11.41)
r =1
 172 Earthquake Resistant Design of Structures

The absolute sum method of modal combination provides a very conservative estimate of the
maximum response in physical coordinates since the time of occurrence of maxima in each mode
in general, is different.

Square root of sum of squares (SRSS) method

If we relax the assumption regarding simultaneous occurrence of peak response in all modes,
and assuming that the natural frequencies are not very closely spaced then the maximum
response in physical coordinate system can be estimated as,

L N
i OPP
(r) 2
1/ 2

vi,max ª MÂ d q r ,max f i (11.42)

MN r =1 Q
It must be emphasized here that these modal combination rules are approximate procedures for
combining the maximum modal responses to get a probable estimate of the maximum response
of physical system. These modal combination rules may be used to estimate the probable
maximum value for any response quantity of interest such as, shear force, bending moment,
drifts, etc. However, care should be taken to ensure that the maximum of each desired response
parameter is first calculated for each mode and then these modal maxima are combined
according to a modal combination rule. An example which illustrates this procedure is given
below.

Example 1 Consider a 3-storey shear building shown in Figure 11.5 with the following
properties:

LM30.0 0.0 0.0 OP LM1.000 1.000 1.000O F 3.88 I

0.0 tonnes, F = 0.548 -1.522 -6.260 P , w = G 9.15 J rad/s
MM 0.0
M = 0.0 30.0
30.0 PQ
P MM0.198 P
12.10 PQ
n
GH15.31JK
N 0.0 N -0.872

Compute the floor displacements, inter-storey drifts, storey shears and overturning moments of
this building when excited by an earthquake. The pseudo-spectral acceleration ordinates of the
earthquake ground acceleration for the three modes are given as Sa = 2.94, 1.57, and 3.93 m/s2.
Assume the storey heights to be 3.0 m and use SRSS rule for combining modal responses.
Solution The modal mass in the r th mode of vibration can be computed as m *r =
T
{f(r)} M{f(r)}. For this problem, the modal masses are m *1 = 40.185, m *2 = 122.306, and m *3
F 1 (r ) T I
{f } Mr can be computed as G1 =
= 5597.928. The mode participation factors Gr = GH m*r
JK
1.303, G2 = –0.342, and G3 = 0.037.
The maximum floor displacements are given by,

vi,max
L 3
ª MÂ {f (r )
Gr ( Sar /w r2 )}2
OP 0.5

MN
r =1
j
PQ
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 173

which, for the current problem are given by

F [{1.303 ¥ (2.94 / 3.88 )} + {-0.342 ¥ (1.57 / 9.15 )} + {0.037 ¥ (3.93 / 15.31 )} ]

2 2 2 2 2 2 0 .5
I
ª G [{0.714 ¥ ( 2.94 / 3.88 )} + {0.520 ¥ (1.57 / 9.15 )} + {-0.232 ¥ (3.93 / 15.31 )} ]
2 2 2 2 2 2 0 .5 JJ
vmax
GG [{0.258 ¥ (2.94 / 3.88 )} + {0.298 ¥ (1.57 / 9.15 )} + {-0.448 ¥ (3.93 / 15.31 )} ] JK
H 2 2 2 2 2 2 0 .5

F 0.255I
= G 0.139J m
GH 0.159JK
The maximum inter-storey drifts are given by,

Dij,max
L 3
ª MÂ n(f (r)
- f (jr ) ) Gr Sar /w r2
2
s OPP
0 .5

MN r =1
i
Q
which, in case of current problem leads to,

F LR(1 - 0.548) ¥ 1.303 ¥ 2.940 U + R(1 + 1.522) ¥ - 0.342 ¥ 1.57 U I

2 2

GG MMNST V ST
3.88 W 2
9.15 W
V J
JJ 2

GG R 3.93 U O JJ 2 0 .5

GG + S(1 + 6.260) ¥ 0.037 ¥

T VP
15.31 W PQ 2

GG LMRS(0.548 - 0.198) ¥ 1.303 ¥ 2.940 UV + RS(- 1.522 + 0.872) ¥ 0.342 ¥ 1.57 UV JJJ
2 2

MNT 3.88 W T2
9.15 W 2
D max ªG JJ
GG R 3.93 U O JJ
2 0.5

GG
+ S(- 6.260 - 12.100) ¥ 0.037 ¥
T 15
VP
.31 W PQ 2

GG LMR(0.198) ¥ 1.303 ¥ 2.940 U + R(-0.872) ¥ - 0.342 ¥ 1.57 U

2 JJ 2

M
N ST 3.88 W
V ST
2
9.15 W
V JJ 2
GG
3.93 U O
0. 5

GH R
+ S(12.100) ¥ 0.037 ¥ V P
2
JJ
T 15.31 W PQ 2
K
F 0.116I
= G 0.089J m
GH 0.051JK
The maximum storey shears are given by,

L R| 3 j
U| OP
2 0.5

Vj,max ª MÂ SG S Â m f (r)
V| P
MN |T r =1
r ar
i =1
ii i
WQ
 174 Earthquake Resistant Design of Structures

or, in vector form

F[{30.0 ¥ 1.0 ¥ 1.303 ¥ 2.94} + {30.0 ¥ 1.0 ¥ - 0.342 ¥ 1.57} I

2 2

GG + {30.0 ¥ 1.0 ¥ 0.037 ¥ 3.93} ] JJ

2 0.5

GG [{30.0 ¥ (1 + 0.548) ¥ 1.303 ¥ 2.94} JJ 2

2
ªG JJ
+ {30.0 ¥ (1 - 1.522) ¥ - 0.342 ¥ 1.57}
Vmax
GG + {30.0 ¥ (1 - 6.260) ¥ 0.037 ¥ 3.93} ]
JJ
2 0.5

GG [{30.0 ¥ (1 + 0.548 + 0.198) ¥ 1.303 ¥ 2.94}

JJ
2

2
GH ++{30{30.0.0¥¥(1(1--6.1260
.522 - 0.872) ¥ - 0.342 ¥ 1.57}
+ 12.100) ¥ 0.037 ¥ 3.93} ] K 2 0 .5

F 116.13I
= G 179.57J kN
GH 204.10JK
The maximum overturning moments are given by,

L |R3 j
|UV OP
2 0 .5

Mj,max ª MÂ SG S Â (h - h ) m f (r )
MN |T
r =1
r ar
i =1
i j ii i
|W PQ
or, in vector form,

F [{30.0 ¥ 0.0 ¥ 1.0 ¥ 1.303 ¥ 2.94} + {30.0 ¥ 0.0 ¥ 1.0 ¥ - 0.342 ¥ 1.57}
2 2 I
GG + {30.0 ¥ 0.0 ¥ 1.0 ¥ 0.037 ¥ 3.93} ] 2 0 .5 JJ
GG [{30.0 ¥ (1 ¥ (9.0 - 6.0) + 0.548 ¥ 0.0) ¥ 1.303 ¥ 2.94} 2
JJ
2
ªG JJ
+ {30.0 ¥ (1 ¥ (9.0 - 6.0) - 1.522 ¥ 0.0) ¥ - 0.342 ¥ 1.57}
Mmax
GG + {30.0 ¥ (1 ¥ (9.0 - 6.0) - 6.260 ¥ 0.0) ¥ 0.037 ¥ 3.93} ] 2 0.5
JJ
GG [{30.0 ¥ (1 ¥ (9.0 - 3.0) + 0.548 ¥ (6.0 - 3.0) + 0.198 ¥ 0.0) ¥ 1.303 ¥ 2.94} 2

2 JJ
GH ++{30{30.0.0¥¥(1(1¥¥(9(.90.0--3.30.)0-) -6.1260
.522 ¥ (6.0 - 3.0) - 0.872 ¥ 0.0) ¥ - 0.342 ¥ 1.57}
¥ (6.0 - 3.0) + 12.100 ¥ 0.0) ¥ 0.037 ¥ 3.93} ] 2 0 .5
K
F 0.0 I
= G 348.39J kN.m
GH 880.55JK
Similarly, the calculations for maximum overturning moment at the base can also be
performed.
In the above-mentioned example, only one component of ground acceleration was consid-
ered for excitation. In general, the structure would be subjected to three mutually orthogonal
 +D=FJAH  Dynamics of Multi-Degree-of-Freedom Systems 175

translational components of ground motions simultaneously at a support point. Computations

for any response quantity of interest for simultaneous excitation by multiple components would,
in general, yield different estimates than for any single ground motion component acting alone.
It is not adequate, for design purpose, to consider the maximum response out of the three
estimates obtained for different ground motion components independently. To find the response
parameters for use in design, the response estimates for excitation by individual components
may be combined together by SRSS rule. For any generic response quantity of interest, say, R,
the value to be adopted for design calculations Rdes can be obtained as,

Rdes = Rx2 + Ry2 + Rz2

where, Rx, Ry, and Rz represent the estimate of response R due to excitation by ground motion
in x, y, and z directions, respectively.

11.6.3 Mode Truncation

Generally, the mathematical models for real civil engineering structural systems may involve
millions of degrees of freedoms, implying that the total number of equations to be solved for
modal coordinates (Equations (11.38)) could be of the same order—a formidable task even for
the powerful desktop computers available today. Fortunately, it is not necessary to include
response in all the modes to get a rational estimate of the total response. Since most of the energy
of the dynamic loads of civil engineering structures (such as earthquake ground motions, wind
forces, ocean waves, etc.) is concentrated in low frequencies (typically < 35 Hz for earthquakes)
the higher modes (with larger natural frequencies) are not excited by these low-frequency
forces.2 Thus it is possible to truncate the modal summation in Equation (11.35) to the sum of
only a few of the lower modes. The total number of terms in such truncated modal summation
seldom exceeds a few hundreds, even in very complex structural systems. Thus the response
vector v can be approximately determined as,
N
v= Â q (t ) {f(r)}
r (11.43)
r =1

where, N << N. The decision about the number of modes to be included in the response
computations may be based on the following two criteria:
(i) All modes having natural frequencies less than or equal to the highest frequency in the
excitation should be included in the modal summation.
(ii) At least 90% of the total mass of the structural system should be included in the dynamic
response computation. This criterion in assessed by considering the cumulative effective
modal mass (= Sr m *r G 2r ) for all modes included in the summation, which should be
more than 90% of the total mass of the system.

2
This inference can be drawn by considering the nature of response of SDOF systems to harmonic excitations. The
dynamic amplification factor for the oscillator response approaches unity as the ratio of frequency of excitation to
natural frequency (w/wn) decreases and the response approaches that for a static case (see Figure 7.5).