Modal analysis

Modal analysis is a set of mathematical techniques for determining the modal characteristics of a dynamical system. Modal analysis can be used to determine the resonant frequencies and vibration responses of a structure and is closely related to failure analysis. The technique often assumes a system to behave linearly as an approximation so that matrix algebra can be used, for which computers are especially well suited. This assumption turns out to be pretty good, in fact, as unperson points out, since stable equilibrium configurations correspond to local minima in potential energy, which can be approximated locally by parabolas, the hallmarks of the simple harmonic oscillator.

Consider three masses m₁ = m, m₂ = M, and m₃ = m constrained to move along a single direction, say, the x axis, with m₁ and m₃ on either side of m₂. Now join m₁ and m₃ separately to m₂ with two springs having spring constants equal to k. Let the three generalized coordinates of the masses be x₁, x₂ and x₃, each of which are zero at the equilibrium positions of the masses.

The total kinetic energy of the system is

T = ½m(dx₁/dt)² + ½M(dx₂/dt)² + ½m(dx₃/dt)²

From Hooke's law, the total potential energy of the system is

V = ½k(x₂-x₁)² + ½k(x₃-x₂)²

V = ½kx₁² + kx₂² + ½kx₃² - kx₁x₂ - kx₂x₃

So the Lagrangian for the system is L = T - V:

L = ½m(dx₁/dt)² + ½M(dx₂/dt)² + ½m(dx₃/dt)² - ½kx₁² - kx₂² - ½kx₃² + kx₁x₂ + kx₂x₃

Now, assuming there aren't any damping forces present, we can obtain the equations of motion from

d/dt (∂L/∂x'_i) - ∂L/∂x_i = 0

Where x'_i is taken to be dx_i/dt. So,

∂L/∂x'₁ = mx₁
∂L/∂x₁ = kx₁ - kx₂

∂L/∂x'₂ = Mx₂
∂L/∂x₂ = 2kx₂ - kx₁ - kx₃

∂L/∂x'₃ = mx₃
∂L/∂x₃ = kx₃ - kx₂

This results in a set of three equations:

md²x₁/dt² - kx₁ + kx₂ = 0
Md²x₂/dt² - 2kx₂ + kx₁ + kx₃ = 0
md²x₃/dt² - kx₃ + kx₂ = 0

Now, we let x = [x₁, x₂, x₃]^T and x" = [d²x₁/dt², d²x₂/dt², d²x₃/dt²]^T, and define the following matrices:

    /m 0 0\
M = |0 M 0|
    \0 0 m/

    / k  -k  0\
K = |-k  2k -k|
    \ 0  -k  k/

Here M is the mass matrix and K is the stiffness matrix, although in general these matrices may be in units other than that of mass or stiffness. With these definitions, we may write our equations of motion as:

Mx" - Kx = 0

Now we may make a series of substitutions in order to simplify calculations. First, we define the matrix square root. The square root of a matrix A, if defined, is A^½ such that A^½ A^½ = A. A convenient way of finding the square root of a matrix is by noting that if A is diagonalizable by an orthogonal matrix P, such that P^TAP = D where D is a diagonal matrix, then A^½ = PD^½P^T. This is because PP^T = P^TP = I and PD^½P^TPD^½P^T = PD^½D^½P^T = PDP^T = A.

For those unfamiliar with the diagonalization process, D is a diagonal matrix with diagonal entries being the eigenvalues of A. The columns of P are the eigenvectors of A corresponding to the eigenvalues of A. The columns must be in the same order as their corresponding eigenvectors in D.

Now D^½ is easy to calculate because only the diagonal terms in D are present, so we take their square roots individually. Finally, let D^-½ be the inverse of D^½.

With these in place, let x = M^-½q, or equivalently let q = M^½x. In our example. We do not need to use diagonalization to calculate M^-½ because M is already a diagonal matrix. Therefore substitute into the equation of motion to obtain

MM^-½q" - KM^-½q = 0

Left multiply by M^-½ to obtain

M^-½MM^-½q" - M^-½KM^-½q = 0

q" - M^-½KM^-½q = 0

Now, let K = M^-½KM^-½, which is the modified stiffness matrix. At this point we obtain

q" - Kq = 0

This is an eigenvalue problem and can be solved directly, but to further simplify the problem, we can turn use some more substitution. First, we must diagonalize K, i.e. find a matrix P, such that P^TKP = D where D is a diagonal matrix. Then let q = Pr so that

Pr" - KPr = 0

This time, left multiply by P^T so that:

P^TPr" - P^TKPr = 0

r" - Dr = 0

With the problem simplified thusly, we can easily express the solutions to the equations of motion, since by changing coordinates from x_i to q_i and then to r_i, we've managed to decouple the systems completely (D is diagonal), so that we may write:

r"_i - λ_i r_i = 0

Writing λ_i = ω_i², where λ_i is the i^th eigenvalue of K and diagonal element of D, and where { ω_i } is the set of natural or resonant frequencies of the system. We can solve each of these equations independently, so that r_i = A_i cos(ω_it + φ_i). Any general set of coordinates in which the system oscillates at a single frequency in each coordinate is called a set of normal coordinates. Thus {r_i} represent the normal coordinates for the system. We can get back our original coordinates by applying the transformations in reverse:

x = M^-½q
q = Pr

So,

x = M^-½Pr

The matrix S = M^-½P is called the matrix of modal shapes, since it encodes characteristic shape of each mode. Writing x = Sr, and denoting the column vectors of S as S_i we get that

x = A₁S₁cos(ω₁t + φ₁) + A₂S₂cos(ω₂t + φ₂) + A₃S₃cos(ω₃t + φ₃) + ...

Thus, the net displacement of each mass, the components of x, are linear superpositions of the displacements brought about by each mode of oscillation. For instance, in our example of three masses connected to eachother with springs, we would find that one mode for which S₁ = [ -1,0,1 ] ^T. This mode is antisymmetric with the outer masses oscillating 90° out of phase with one another, and the inner mass remaining stationary. This mode has its own frequency ω₁ and is independent of and uninfluenced by the other modes.

In order to solve the problem in its entirety, then, we need determine the values of A_i and φ_i. This can be done as long as the initial conditions for the system are given, or any other independent conditions that fix these constants. Normally the initial conditions are given in terms of the system's initial displacements and velocities.

In the above development we used an approach that allowed us to find both the modal shapes and natural frequencies. If we were only after the frequencies ω_i we would have used the equation Mx" - Kx = 0 directly, by assuming that x_i was of the form A_i cos(ω_it + φ_i). In this method we consider each mode separately, by assuming that the amplitudes of all the other modes are zero (by some appropriate initial conditions). We can then substitute into the above equation to obtain

ω_i²Mx + Kx = 0

(ω_i²M + K)x = 0.

Since non-trivial solutions require that det(ω_i²M + K) = 0, we can find { ω_i } as the roots of the characteristic equation.