Vibration of Multi Degree of Freedom Systems

Data: $L = 1.5$ m, $d = 0.05$ m, $E = 200\times10^9$ Pa, masses 20, 30, 15 kg at $x = 0.4, 0.8, 1.2$ m.

$$I = \frac{\pi d^4}{64} = \frac{\pi (0.05)^4}{64} = 3.0680e-07\ \text{m}^4, \qquad EI = 61359\ \text{N·m}^2$$

Influence coefficients

For a unit load at distance $a$ from the left end (with $b = L - a$), the deflection at $x \le a$ is

$$\alpha(x,a) = \frac{b\,x\,(L^2 - b^2 - x^2)}{6EIL}, \qquad \alpha_{ij} = \alpha_{ji}\ (\text{Maxwell})$$

The influence-coefficient matrix (m/N, $\times10^{-6}$):

$$[\alpha] = 10^{-6}\begin{bmatrix} 0.701 & 0.811 & 0.435 \\ 0.811 & 1.136 & 0.661 \\ 0.435 & 0.661 & 0.469 \end{bmatrix}$$

(i) Rayleigh's method

Static deflections under the weights $W_i = m_ig$: $y_i = \sum_j\alpha_{ij}W_j$.

$$\omega_n^2 = \frac{g\sum m_iy_i}{\sum m_iy_i^2} = \frac{9.81 \times 31.7554\times10^{-3}}{16.16626\times10^{-6}} = 19270\ \text{rad}^2/\text{s}^2$$ $$\omega_n = 138.82\ \text{rad/s}, \qquad f_n = 22.09\ \text{Hz}$$

(ii) Dunkerley's method

The frequency with each mass acting alone is $\omega_i = \dfrac{1}{\sqrt{\alpha_{ii}m_i}}$:

$$\frac{1}{\omega_n^2} = \sum\frac{1}{\omega_i^2} = 5.5136e-05 \Rightarrow \omega_n = 134.67\ \text{rad/s}, \qquad f_n = 21.43\ \text{Hz}$$

Comparison

Rayleigh's method gives an upper bound, Dunkerley's method a lower bound on the true fundamental frequency. A check by solving the eigenvalue problem exactly gives $f = 22.09$ Hz, which lies between them.

Answer: Rayleigh: $f_n \approx 22.09$ Hz; Dunkerley: $f_n \approx 21.43$ Hz.

Data: $m_1 = 4$, $m_2 = 2$, $m_3 = 1$ kg; $k = 2000$ N/m for each spring.

Flexibility and dynamic matrices

For this chain, the flexibility influence coefficient is $a_{ij} = \dfrac{\min(i,j)}{k}$:

$$[a] = \frac{1}{2000}\begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 2 \\ 1 & 2 & 3 \end{bmatrix}\ \text{m/N}$$

The free vibration equation $\{X\} = \omega^2[a][M]\{X\}$ is rewritten with the dynamic matrix $[D] = [a][M]$:

$$[D]\{X\} = \lambda\{X\}, \qquad \lambda = \frac{1}{\omega^2}$$ $$[D] = \frac{1}{2000}\begin{bmatrix} 4 & 2 & 1 \\ 4 & 4 & 2 \\ 4 & 4 & 3 \end{bmatrix}$$

Iteration

Start with the trial vector $\{1, 1, 1\}^T$, multiply by $[D]$ and normalise by dividing by the first element. The first element of $[D]\{X\}$ divided by the first element of $\{X\}$ is the estimate of $\lambda$.

The mode shape and $\lambda$ have converged. From the last cycle,

$$\lambda_1 = 4.4199\times10^{-3}\ \text{s}^2, \qquad \omega_1 = \frac{1}{\sqrt{\lambda_1}} = 15.04\ \text{rad/s}, \qquad f_1 = \frac{\omega_1}{2\pi} = 2.394\ \text{Hz}$$

The first mode shape (normalised to $X_1 = 1$) is $\{1,\ 1.548,\ 1.745\}^T$. All three masses move in phase, with the largest motion at the free end.

(An exact eigenvalue solution gives $\omega_1 = 15.04$ rad/s, which confirms the result.)

Answer: $\omega_1 \approx 15.04$ rad/s ($f_1 \approx 2.394$ Hz); mode shape $\{1,\ 1.548,\ 1.745\}$.

Rayleigh-Ritz method

It improves Rayleigh's quotient by using a trial function with several adjustable constants.

1. Assume the mode shape as a series of admissible functions (satisfying the geometric boundary conditions):

$$y(x) = c_1\phi_1(x) + c_2\phi_2(x) + \cdots + c_n\phi_n(x)$$

2. Write the Rayleigh quotient $\omega^2 = \dfrac{U_{max}}{T^*_{max}}$, where $U_{max}$ is the maximum strain energy and $T^*_{max}$ is the kinetic energy expression with $\omega^2$ taken out.
3. Make $\omega^2$ stationary with respect to every constant:

$$\frac{\partial}{\partial c_i}\left(U_{max} - \omega^2T^*_{max}\right) = 0, \qquad i = 1,2,\ldots,n$$

This gives $\left([K] - \omega^2[M]\right)\{c\} = 0$, an eigenvalue problem of order $n$. 4. The $n$ roots approximate the lowest $n$ natural frequencies, always as upper bounds. More terms give better accuracy, and the higher modes are less accurate.

Finite difference method

It replaces the derivatives in the differential equation by differences between displacements at equally spaced stations.

For a beam divided into stations a distance $h$ apart:

$$\left.\frac{d^2y}{dx^2}\right|_i \approx \frac{y_{i+1} - 2y_i + y_{i-1}}{h^2}, \qquad \left.\frac{d^4y}{dx^4}\right|_i \approx \frac{y_{i+2} - 4y_{i+1} + 6y_i - 4y_{i-1} + y_{i-2}}{h^4}$$

Procedure:

1. Divide the system into $n$ segments and write the governing equation, e.g. $EI\dfrac{d^4y}{dx^4} = \omega^2\rho Ay$, at every interior station.
2. Apply the boundary conditions (e.g. $y = 0$ and $y'' = 0$ at a simple support) using extra fictitious stations if necessary.
3. The result is $[A]\{y\} = \omega^2[B]\{y\}$, solved as an eigenvalue problem.

Accuracy improves as $h$ decreases, but the order of the matrix grows. For a uniform beam, 4 to 6 segments already give the first frequency within a few percent.