Numerical solution methods for natural frequencies and mode shapes in relation to structural dynamics during earthquakes: Abstract: This chapter discusses basic solution schemes as well as approximate methods for finding natural frequencies and mode shapes. The methods include Vianelloâ€'Stoodala power method, transfer matrix method, Jacobi method, Holzer method, Rayleighâ€™s approximation and Dunkerleyâ€™s approximation. Some of the approximate methods will lead to upper bound solutions and some lower bound solutions. A relevant program in MATHEMATICA is also given.
Numerical
solution methods for natural frequencies and mode shapes in relation to
structural dynamics during earthquakes
Abstract: This
chapter discusses basic solution schemes as well as approximate methods
for finding natural frequencies and mode shapes. The methods include
Vianello�'Stoodala power method, transfer matrix method, Jacobi method, Holzer
method, Rayleigh’s approximation and Dunkerley’s approximation. Some of the
approximate methods will lead to upper bound solutions and some lower bound
solutions. A relevant program in MATHEMATICA is also given.
Key words: banded matrix, sweeping
technique, deflation, transfer matrix, Holzer method.
Introduction
The important step in the dynamic
analysis of a multiple-degrees-of-freedom (MDOF) system is quite often the
solution of the eigenvalue problem, or the determination of the system natural
frequencies and the corresponding normal vibration modes. This is particularly
true if a mode superposition analysis is to be conducted. Several
procedures for solving the eigenvalue problem have been discussed in many
books. Predicting or finding the roots of the characteristic polynomial as
discussed in Chapter 10 is satisfactory only for systems having few degrees of
freedom. For large MDOF systems, extracting the roots of the characteristic
polynomial requires computational effort and is quite often an indeterminable
task. This chapter discusses the basic solution schemes as well as approximate methods
for finding frequencies.
General solution methods for
eigen problems
In structural dynamics, the basic eigen problem for an MDOF
system having ‘n’ degrees of freedom is represented as
[k]{φ} = λ[m]{φ} … . . . . 11.1
Let [k] be the stiffness matrix of order ‘n’ and
[m], the mass matrix, also of order ‘n’. For most structural
systems, [k] is normally banded matrix and [m] is a
diagonal matrix for a lumped mass formulation (without rotary inertia) coupling
or a narrowly banded matrix for a consistent mass formulation.
There are ‘n’ eigenvalues and ‘n’ eigenvectors
satisfying the above equation. The rth eigen pair is determined by (λ_{r},{φ}_{r}). In dynamic problems,
the eigenvalues are the square of the natural frequencies ω _{n}^{2}
such that
0 < λ_{1} < λ_{2} <…< λ_{n ……….
.. .} 11.2
The dynamic response of MDOF systems having a large number of
degrees of freedom is generally confined to a relatively small subset of the
lowest vibration modes of the system. Therefore for such systems, only ‘p’
eigen pairs need to be solved for where p << n. The
solution of p eigen values and the corresponding eigenvectors can be
written as
The
majority of the eigen problem solution techniques can be classified as
• vector
iteration methods;
• transformation
methods;
• polynomial
iteration methods.
Clearly all the methods are
iterative in nature because the solution of the eigen problem as defined in Eq.
11.1 is tantamount to solving the characteristic polynomial of order ‘n’.
Since explicit formulas for the determination of roots to the characteristic
polynomial having an order higher than 4 do not exist, an iterative solution is
mandatory.
The main essence of each method is very distinctive. The
vector iteration methods are based on the property
Let [K] and [M] be modal stiffness and modal
mass respectively. The polynomial iteration is based on the property that the
characteristic polynomial is a function of λ_{r}, and
The characteristic polynomial is of order ‘n’. Solution
of characteristic polynomial has been discussed in Previous Pages.
Vector iteration
technique
1 Vianello and
Stoodala method (power method)
This iterative method can be applied to extract the highest
eigenvalue of either symmetric or unsymmetric matrix of any order. Consider the
homogeneous equation written in the form of
Assume any vector {x}_{0} and multiplying with
[A] matrix gives {y}_{1} can be written in terms of λ_{1}{x}_{1}
by taking highest element (associated vector) outside and this can be used in
the next iteration as
It can be seen that the
eigenvalue λ as well
as eigenvector will converge. The iteration can be stopped when |λ_{n}_{+1} �' λ_{n}| < ε
(tolerance) and λ is the
highest eigenvalue in the absolute sense.
Consider the problem in dynamics
[D]_{1} is called the dynamic matrix.
Application of the power method will converge to (1/λ)_{max} or λ_{min} = λ_{1} · λ_{1} is
called the least dominant eigenvalue and {ϕ}_{1}
is called least dominant eigenvector. However the method may be modified to
calculate eigenvalues and the corresponding eigenvectors for higher modes by matrix
deflation or deflation of the iteration vectors.
Assume using the power method we
get the (1/λ)_{max}
or λ_{min} = λ_{1} and {φ}_{1} pair. The basic
premise for vector deflation is that, for an iteration vector to converge to a
required eigenvector, the iteration vector must be orthogonal to the
eigenvector. Therefore, for this case at hand, this can be interpreted as
meaning that if the iterative vector is orthogonalized to the eigenvectors
already calculated (for example {φ}_{1}) the vector is precluded and occurs instead to
another (higher) eigenvector. More succinctly, an eigen pair other than (λ_{1}, {φ}_{1}) becomes the least
dominant eigen pair.
2 Method 1
sweeping technique
To find the second eigenvalue assume a trial vector that is a
linear combination of all eigenvectors.
In the trial vector assumed in
Eq. 11.13, let us sweep out the effect of mode 1 as
If {φ}_{1}
is the normalized eigenvector, as proved in Chapter 10, {φ }_{1}^{T} [ m]{φ}_{1} = 1 because of the normalization principle. In that case
dynamic matrix is written as [D]_{2} = [D]_{1}[S]_{1}
where [S]_{1} is given as in Eq. 11.18 as
In Eq. 11.19, [S]_{1}
is called sweeping matrix for the first mode. The purpose of [S]_{1}
is to eliminate or sweep out the effect of mode 1. Similarly [S]_{n}
is to sweep out the effect of modes for 1 to n and allow the mode (n
+ 1) to become least dominant.
3 Method 2 �'
deflation method
Assume we find (λ_{1}, {φ}_{1}) pair to obtain
some eigenvalue using the ortho-normal principle.
The power method will converge to
the second eigenvalue and the corresponding normalized eigenvector {ϕ}_{2} can be
obtained.
To extract the third eigenvalue,
Using the power method, the
iteration will converge to third eigenvalue and the corresponding normalized
eigenvector can be obtained.
Let us use the same example and assume we obtain the highest
eigenvalue as 2.8453 and the corresponding normalized eigenvector as given in
method 1. Substituting the corresponding values in Eq. 11.25, we get
Using the first normalized eigenvector we obtain the equation
x_{1} = �'1.611x_{2}
�' 1.656x_{3 }(a)
Using the second normalized eigenvector we obtain the equation
The procedure can be programmed very easily in the EXCEL
package.
Jacobi’s method
The matrix iteration method
discussed in Section 11.3 produces the eigenvalues and eigenvectors of matrix [D]_{1}
one at a time. Jacobi’s method is also an iterative method but produces all the
eigenvalues and eigenvectors of the matrix [D]_{1}
simultaneously. [D]_{1} is a real symmetrix matrix and has only
real eigenvalues. There is an orthogonal matrix [R] such that [R]^{T}[D]_{1}[R]
is a diagonal matrix. The diagonal elements are the eigenvalues and the columns
of [R] is generated as a product of several rotation matrices as
[R] = [R_{1}][R_{2}][R_{3}]… ---- 11.36
Consider the highest off-diagonal term of matrix [D]_{1}.
Let it be d_{ij}. Find
Again find the highest off-diagonal term and form [R_{2}]
matrix. While making this off-diagonal term as zero, it introduces non-zero
contributions to formerly zero positive. However, successive matrices of the
form
[R_{2}]^{T}[R_{1}]^{T}[D]_{1}[R_{1}][R_{2}] - - - - - - 11.41
[R_{3}]^{T}[R_{2}]^{T}[R_{1}]^{T}[D]_{1}[R_{1}][R_{2}][R_{3}] - - - - - - 11.42
converges to the required diagonal form. Find matrix [R]
such that
[R] = [R_{1}][R_{2}][R_{3}]…[R_{n}] - - - - - - 11.43
is the required eigenvector.
Consider Example 11.1 with mass matrix as [I]. The
dynamic matrix is written as
Transfer matrix
method to find the fundamental frequency of a multi-storeyed building (shear
frame)
Consider a shear frame shown in Fig. 11.1 consisting of ‘n’
storeys. Consider the free body diagram of the ith storey as shown in
Fig. 11.2. Here ‘m’ is the mass, ‘k’ is the stiffness, ‘V’
is the shear, ‘ω’ is the
natural frequency and ‘v’ is the displacement.
Considering inertia force, the shear in the ith storey
can be written in terms of shear of (i+1)th storey.
Initially the displacement at the
top storey level and the natural frequency are both assumed. By using the
transfer matrix method, it is possible to find the displacement at the base as
well as shear in the base storey. If the support displacement is not zero, a
new value for the natural frequency is assumed and the procedure is repeated
till we get the value of the base displacement as zero.
Example 11.2
Find the natural frequency of the
three storeyed shear building as shown in Fig. 11.3, given the mass and the
stiffness.
Solution
Since ‘MATHEMATICA’ can solve the problem using symbolic
processing, we will write in symbolic form. Transfer matrix at top frame (
assume p = ω _{n}^{2}
)
One we get natural frequency one can also get the mode shape
as shown below.
Holzer method for
torsional vibrations
The Holzer method falls under the
determinant search technique. The method can be applied to rectilinear or
angular motions for damped as well as undamped systems. The method is best
suited for systems where the components are arranged along the basic axis. Let
us consider the torsional vibration for shafts as shown in Fig. 11.4.
The torsional moment equilibrium at node 1 can be written as
Approximate methods
for finding the natural frequencies
1 Rayleigh’s
quotient
In many practical situations
involving MDOF systems, only the accurate estimation of the fundamental
frequency is required. In such cases, laborious calculations to extract all the
normal vibration modes of the system are not warranted and the approximate methods
are desirable. This section discusses two approximate methods for estimating
the fundamental frequency of MDOF systems.
The first method, Rayleigh’s
method, is an upper bound method based on energy principles and stiffness
approach. The second method, Dunkerley’s approximation, is based
on the flexibility of the system eigenvalue problem and therefore
provides lower bound estimation of the fundamental frequency. Thus the upper
bound estimation of the fundamental frequency provided by Rayleigh’s method can
be complemented by the lower bound estimation afforded by Dunkerley’s
approximations to envelope true fundamental frequency.
Consider that an undamped single-degree-of-freedom (SDOF) mass
spring system is in free harmonic motion given by
2 Rayleigh’s
quotient method to MDOF
In the above we have discussed
Rayleigh’s method to determine the fundamental frequency of an SDOF system.
Application of the Rayleigh’s method to determine the fundamental frequency of
an MDOF system is presented in this section.
Consider the eigenvalue problem of the MDOF system represented
by equation
In Eq. 11.70b, the denominator is related to the kinetic
energy for the rth mode and the numerator is related to the potential
energy, or the strain energy of the rth mode. If the modal vector {φ}_{r} is replaced
with any arbitrary vector {A}, Eq. 11.70b is written as
where R({A}) is a scalar quantity referred to as
Rayleigh’s quotient. It is evolved from Eq. 7.51 that Rayleigh’s
quotient is dependent upon the known matrix [m] and [k] and the
unknown arbitrary vector {A}. Obviously if {A}
coincides with one of the systems normal modes then λ_{R} is the corresponding
eigenvalue or normal frequency of the system. A very important property of
Rayleigh’s quotient is
and it also follows that for any vector {A}, if [K]
is positive definite, R({A}) > 0. If [K] is positive
semi-definite,
Equation 11.65 thus indicates that the Rayleigh’s quotient is
never lower than the fundamental eigenvalue, and furthermore the minimum value
the Rayleigh’s quotient can assume is that of the fundamental eigenvalue
itself. Therefore, Rayleigh’s quotient is very good technique to estimate the
fundamental frequency of MDOF systems. A reasonable estimate for the vector {A}
corresponding to the fundamental mode is the vector of static displacement
resulting from subjecting the masses in the system to forces proportional to
their weights. Many seismic design code present expressions to estimate the
fundamental frequency of high-rise building based on this concept. The natural
frequency thus obtained is called the Rayleigh frequency ω_{R}
expressed as
The accuracy of the Rayleigh
frequency ω_{R} depends
entirely on the displacement vector {A} used to represent the vibration
mode shape. In principle, any vector {A} may be selected which satisfies
the geometric boundary conditions. However, any vector other than the true
modal vector requires the action of additional external constraints to maintain
equilibrium, which would in turn stiffen the structure, resulting in increased
computed frequency. Therefore the true vibration mode will yield the lowest
frequency
obtained by Rayleigh’s method. Hence the approximation
yielding the lowest frequency for a particular case is the best result.
Steps
1. Estimate
the fundamental mode of vibration. This may be done either by assuming the
displacement of the modes directly or computing the displacement from the
associated forces.
If the frequency is computed for
several displacement-assumed configurations, the smallest of the computed
values will be close to the exact value of ω_{n} and the associated configuration
is closest to the actual configuration.
Example 11.6
Determine the fundamental frequency of the shear frame shown
in Fig. 11.3 by the improved Rayleigh method.
Solution
R00 method
Assume the mode shape as shown in Fig. 11.9. Maximum potential
energy:
R_{01} method
By calculating the shear in each floor let us improve the mode
shape. The shear in each floor can be calculated as (see Fig. 11.10)
Equating this improved potential energy to previously
calculated kinetic energy
0.009 06 ω _{n}^{4} = 2.25ω _{n}^{2}
Or
ω _{n}^{2} = 234.375;
ω _{n} = 15.309
R_{11} method
Improve the kinetic energy. Velocity at each storey level:
v_{1} = 0.0025 ω _{n}^{3} ; v _{2} = 0.004 58 ω _{n}^{3} ; v_{3} = 0.006 25 ω _{n}^{3}
T_{max(improved)} = ω _{n}^{6} /2[2 �- 0.0025^{2} + 1.5 �- 0.004 58^{2} + 1 �- 0.006 25^{2} ]
= 4.15135 �- 10 ^{�'}^{5} ω _{n}^{6}
Equating maximum kinetic energy to maximum potential energy,
we get
0.009 06 ω _{n}^{4} = 4.15135 �- 10 ^{�'}^{5} ω _{n}^{6}
ω_{n} =
14.77rad/s
Dunkerley’s approximation
It is another approximate method for estimating the
fundamental frequency for MDOF systems. The method yields accurate results for
systems for which
damping is negligible and the
natural frequencies are well separated. Dunkerley’s equation provides a ‘lower
bound’ estimate with fundamental frequency and is therefore complementary with
the Rayleigh method that provides an ‘upper bound’ estimate with fundamental
frequency.
To derive Dunkerley’s equation consider the equation
where [a] is the flexibility matrix.
The frequency equation is
obtained by expanding the determinant of the characteristic matrix in Eq.
11.101.
Let us consider a two-degrees-of-freedom system with lumped
mass diagonal matrix. Thus the resulting characteristic determinant becomes
Expanding Eq. 11.102 results in the system frequency equation,
i.e. second order equation in λ = 1/ω^{2} given
by,
If the roots are 1/ ω_{1}^{2} , 1/ ω _{2}^{2}
The relationship represented by Eq. 11.104 also holds true for
systems having ‘n’ degrees of freedom. Extending this ‘n’ degrees
of freedom system
Dunkerley’s approximation to the fundamental frequency is made
on the assumption that if the fundamental frequency ‘ω’ is much lower than the higher
harmonics (ω_{2},… ω_{n}) then
the terms on the left hand side 1/ ω _{2}^{2} … 1/ ω _{n}^{2}
can be calculated. The elimination of these terms yields an estimate of 1/ ω_{1}^{2} which is
higher than the true value thereby making the estimate of ‘ω_{1}’ lower
than the exact value of fundamental frequency. The Dunkerley’s lower bound
estimate of ‘ω_{1}’ is
approximated to
In Eq. 11.106 the term ‘a_{ii}m_{i}’
represents the contribution of each mass to 1/ω_{1}^{2} in the
absence of all other masses. Thus
where ω _{ii}^{2}
is the natural frequency of an SDOF system with mass ‘m_{i}’
acting alone at state i. Hence Dunkerley’s equation is given by
Summary
In this chapter the sweeping
technique combined with power method and transfer matrix methods have been
discussed to find the natural frequencies of n-degrees-of-freedom
system. In addition, Rayleigh’s coefficient method and Dunkerley’s approximate
methods are also discussed to find the approximate fundamental frequency of an n-degrees-of-freedom
system.