Notebook

Universidade Federal do Rio Grande do Sul (UFRGS)¶

Programa de Pós-Graduação em Engenharia Civil (PPGEC)

PEC00025: Introduction to Vibration Theory¶

Class 14 - Vibration of beams¶

1. The vibrating beam equation
2. Free vibration solution
3. Vibration modes and frequencies
4. Solution by approximation
5. Assignment

Prof. Marcelo M. Rocha, Dr.techn. (ORCID)
Porto Alegre, RS, Brazil

In [1]:

# Importing Python modules required for this notebook
# (this cell must be executed with "shift+enter" before any other Python cell)

import numpy as np
import matplotlib.pyplot as plt

1. The vibrating beam equation ¶

The static analysis of a beam under Bernoulli's hypothesis (plane sections remain plane after beam deformation) leads to the well known differential equations:

$\begin{align*} \frac{dQ}{dx} &= -q(x) \\ \frac{dM}{dx} &= Q(x) \\ EI \psi^\prime &= M(x) \end{align*}$

where $q(x)$ is the distributed transversal loading, $Q(x)$ is the shear, $M(x)$ is the bending moment, $\psi(x)$ is the section rotation, and $EI$ is the flexural stiffness (regarded hereinafter as constant along the beam length).

Disregarding shear strains, $\gamma(x) = 0$ , implies that the section rotation is approximated as:

$\psi \approx -w^\prime(x)$

what implies that:

$\begin{align*} EI w^{\prime\prime} &\approx -M(x) \\ EI w^{\prime\prime\prime\prime} &\approx q(x) \end{align*}$

where $w(x)$ is the beam transversal displacement, also called elastic line. The solution of this last differential equation is straighforward once the load $q(x)$ and the boundary conditions (two for each beam extremity) are specified.

We shall now include the inertial forces in these equations, as well as regard now the section mean shear strain,

$\gamma(x) = \psi(x) + w^\prime(x) = \frac{Q(x)}{GA_{\rm s}}$

as relevant, where $GA_{\rm s}$ is the shear stiffness (also regarded hereinafter as constant along the beam length). Although $\gamma(x)$ is indeed negligible for actual slender beams, the following analysis may also be applied to other slender structures, like tall buildings, trusses, etc., for which we may define equivalent stiffnesses, $\left(EI\right)_{\rm eq}$ and $\left(GA_{\rm s}\right)_{\rm eq}$ . The dynamic equilibrium equations now become:

$\begin{align*} Q^\prime &= -q + \mu \ddot{ w} \\ M^\prime &= Q - i_\mu \ddot{\psi} \end{align*}$

where $\mu$ is the beam mass per unit length and $i_\mu$ is the cross section rotational inertia per unit length. Derivating the equation for $\gamma$ and replacing the solution $EI \psi^{\prime} = M(x)$ gives the elastic line equation accounting for shear strains:

$w^{\prime\prime} = -\frac{M}{EI} + \frac{Q^\prime}{GA_{\rm s}}$

Now we replace the shear (with inertial loads):

$w^{\prime\prime} = -\frac{M}{EI} + \frac{-q + \mu \ddot{ w}}{GA_{\rm s}}$

and derivate the whole equation replacing for $M^\prime$ :

$w^{\prime\prime\prime} = \frac{i_\mu \ddot{\psi} - Q}{EI} + \frac{\mu \ddot{w}^\prime - q^\prime}{GA_{\rm s}}$

The angular acceleration, $\ddot{\psi}$ , may be safely disregarded, for the rotations are usually very small. Derivating the equation a last time and replacing for $Q^\prime$ finally gives:

$EI w^{\prime\prime\prime\prime} = q - \mu \ddot{ w} + \frac{EI}{GA_{\rm s}} \, \left(\mu \ddot{ w}^{\prime\prime} - q^{\prime\prime} \right)$

which is the dynamic elastic line equation for a constant section beam under forced vibration due to dynamic load $q(x,t)$ , with shear deformation accounted for (although plane section hypothesis still kept along). The last term may be disregarded whenever the shear stiffness is much larger that the bending stiffness.

2. Free vibration solution ¶

In this section, we take the vibrating beam equation derived above, disregard the shear deformation and look for free vibration solution, which implies that $q(x, t) = 0$ . The beam equation becomes simply:

$EI w^{\prime\prime\prime\prime} = - \mu \ddot{w}$

To solve this equation we separate the time and space independent variables through the hypothesis:

$w(x,t) = w_0 \sin \left(\omega t + \theta \right)\, \varphi(x)\,$

which is pretty much alike we have previously done for multiple degrees of freedom systems. The free vibration equilibrium equation then become:

$\varphi^{\prime\prime\prime\prime} - p^4 \varphi = 0$

where we have defined:

$p^4 = \left(\frac{\mu}{EI}\right) \omega^2$

It can be shown that, in the general case, the space dependent function $\varphi(x)$ has the form:

$\varphi(x) = C_1 \left(\cos px + \cosh px \right) + C_2 \left(\cos px - \cosh px \right) + C_3 \left(\sin px + \sinh px \right) + C_4 \left(\sin px - \sinh px \right)$

The corresponding space derivatives will be required to apply the boundary conditions:

$\begin{align*} \varphi^\prime(x) = p^1&\left[C_1 \left(-\sin px + \sinh px \right) + C_2 \left(-\sin px - \sinh px \right) + C_3 \left( \cos px + \cosh px \right) + C_4 \left( \cos px - \cosh px \right)\right] \\ \varphi^{\prime\prime}(x) = p^2&\left[C_1 \left(-\cos px + \cosh px \right) + C_2 \left(-\cos px - \cosh px \right) + C_3 \left(-\sin px + \sinh px \right) + C_4 \left(-\sin px - \sinh px \right)\right] \\ \varphi^{\prime\prime\prime}(x) = p^3&\left[C_1 \left( \sin px + \sinh px \right) + C_2 \left( \sin px - \sinh px \right) + C_3 \left(-\cos px + \cosh px \right) + C_4 \left(-\cos px - \cosh px \right)\right] \\ \varphi^{\prime\prime\prime\prime}(x) = p^4&\left[C_1 \left( \cos px + \cosh px \right) + C_2 \left( \cos px - \cosh px \right) + C_3 \left( \sin px + \sinh px \right) + C_4 \left( \sin px - \sinh px \right)\right] \end{align*}$

The last equation above proves that the assumed general solution is correct.

Now, to have a particular solution for the vibrating beam the kinematic boundary conditions must be applied. Let us assume a cantilever beam, fixed at the left end ( $x = 0$ ) and free at the right end ( $x = L$ ).

The corresponding boundary conditions are:

$\begin{align*} \varphi(0) &= 0 \\ \varphi^ \prime(0) &= 0 \\ \varphi^{\prime\prime}(L) &= 0 \\ \varphi^{\prime\prime\prime}(L) &= 0 \end{align*}$

The two last conditions implies that bending moment and shear force are zero at the right end, respectively. Applying these conditions at the corresponding derivatives:

$\begin{align*} \varphi(0) &= C_1 \left( 1 + 1 \right) + C_2 \left( 1 - 1 \right) + C_3 \left( 0 + 0 \right) + C_4 \left( 0 - 0 \right) = 0\\ \varphi^\prime(0) &= p \left[C_1 \left(-0 + 0 \right) + C_2 \left(-0 - 0 \right) + C_3 \left( 1 + 1 \right) + C_4 \left( 1 - 1 \right)\right] = 0 \end{align*}$

which implies that $C_1 = 0$ and $C_3 = 0$ . The other two conditions become:

$\begin{align*} \varphi^{\prime\prime}(L) &= p^2 \left[C_2 \left(-\cos pL - \cosh pL \right) + C_4 \left(-\sin pL - \sinh pL \right)\right] = 0\\ \varphi^{\prime\prime\prime}(L) &= p^3 \left[C_2 \left( \sin pL - \sinh pL \right) + C_4 \left(-\cos pL - \cosh pL \right)\right] = 0 \end{align*}$

These two equations can be put into matrix form as:

$\left[ \begin{array}{cc} \left(-\cos pL - \cosh pL \right) & \left(-\sin pL - \sinh pL \right) \\ \left( \sin pL - \sinh pL \right) & \left(-\cos pL - \cosh pL \right) \end{array} \right] \left[ \begin{array}{c} C_2 \\ C_4 \end{array} \right] = \left[ \begin{array}{c} 0 \\ 0 \end{array} \right]$

In order to obtain a non trivial (non zero) solution for the unknown coefficients $C_2$ and $C_4$ , the determinant of the coefficients matrix must be zero. This condition yields a nonlinear equation to be solved for $pL$ . We can use the HP Prime for this purpose, as shown in the figure below.

There will be infinite solutions $\alpha_k = \left( pL \right)_k$ , $k = 1, 2, \dots, \infty$ , each one associated to a vibration frequency and a modal shape, $\left[ \omega_k, \varphi_k(x) \right]$ . The natural vibration frequencies are obtained by recalling the definition of $p$ , what finally gives:

$\omega_k = \left( \frac{\alpha_k}{L} \right) ^2 \sqrt{\frac{EI}{\mu}}$

For instance, the fundamental vibration frequency, $f_{\rm n}$ , is given by:

$f_{\rm n} \approx \frac{1}{2\pi} \left( \frac{1.8751}{L} \right) ^2 \sqrt{\frac{EI}{\mu}}$

which is a very useful formula for estimating the fundamental frequency of slender piles and towers with constant cross section.

In [2]:

x = np.linspace(0, 10, 1000)

3. Vibration modes and frequencies ¶

The table below was taken from a Sonderausdruck (special edition) of the german Betonkalender (concrete almanac), 1988. It summarizes the solutions for some other support conditions of slender beams. If more accuracy is desiredfor the $\alpha_k$ constants, one can solve the so-called characteristic equation with the help of a calculator like the HP Prime, ou HP 50G.

The characteristic equations are the determinants of respective coefficients matrix, which can also be solved with the fsolve() method from scipy, as shown below.

In [12]:

def char_eq(x):
    
    x = x[0]
    A = np.array([[-np.cos(x)-np.cosh(x), -np.sin(x)-np.sinh(x)], 
                  [ np.sin(x)-np.sinh(x), -np.cos(x)-np.cosh(x)]])
 
    return np.linalg.det(A)              # from coefficientes matrix

#   return np.cos(x)*np.cosh(x) + 1      # from characteristic equation

#-----------------------------------------------------------------------

from scipy.optimize import fsolve

ak = fsolve(char_eq, 1)

print('Cantilever beam frequency parameter: {0}'.format(ak[0]))

Cantilever beam frequency parameter: 1.8751040687119394

Observe that the result is exactly the same (within the required precision) as previously obtained with the HP Prime. One can use directly the characteristic equation, or one can program the determinant calculation using np.linalg.det() in the user function char_eq above.

4. Solution by approximation ¶

The solutions for the beam vibration frequencies presented above may also be calculated by means of Rayleigh quotient, as long as an educated guess for the $\varphi(x)$ function is assumed. As an example, let us take a simply supporte beam, for which we assume:

$\varphi(x) = 4x \left( \frac{L - x}{L^2} \right)$

with $\varphi(L/2) = 1$ , which obviously is not the correct modal shape.

We also assume that the beam is subjected to a constant distributed load, $q$ corresponding to its self weight. The maximum displacement at the beam center is known to be:

$w_{\rm max} = \frac{5 q L^4}{384EI}$

In the following we shall estimate both the fundamental vibration frequency and the maximum displacement with the assumed modal shape. The reference kinetic energy is given by:

$T_{\rm ref} = \frac{1}{2} \int_0^L {\mu \varphi ^2(x) \, dx}$

while for the elastic potential energy we need the curvature function, $\varphi^{\prime\prime}(x)$ :

$V = \frac{1}{2} \int_0^L {EI \left[ \varphi^{\prime\prime}(x) \right] ^2 \, dx} = \frac{1}{2} \int_0^L {q w(x) \, dx}$

On the other hand, the modal properties are evaluated with a continuous version of the same formula presented on Class 11. In particular, the modal mass and the modal load are:

$\begin{align*} \vec\phi_k^{\intercal}{\mathbf M} \vec\phi_k &\implies M_k = \int_0^L{\mu \, \varphi^2(x) \, dx} \\ \vec\phi_k^{\intercal} \vec{F} &\implies F_k = \int_0^L{ q \, \varphi(x) \, dx} \end{align*}$

The static modal response can be calculated as:

$u_k = F_k/K_k$

where $K_k = \omega_k^2 M_k$ .

Let us apply this approach to the beam example above. For the assume modal shape we have the curvature:

$\varphi^{\prime\prime}(x) = -\frac{8}{L^2}$

Hence:

$\begin{align*} T_{\rm ref} &= \frac{1}{2} \int_0^L {\mu \left[ \frac{4x(L - x)}{L^2} \right]^2 \, dx} = \frac{ 4}{15 }\mu L \\ V &= \frac{1}{2} \int_0^L { EI \left[ -\frac{8} {L^2} \right]^2 \, dx} = \frac{32}{L^3} EI \end{align*}$

The Rayleigh quotient results:

$\omega_k^2 = \frac{V}{T_{\rm ref}} = \frac{32EI}{L^3} \frac{15}{4\mu L} = \frac{ 120}{L^4} \left( \frac{EI}{\mu} \right)$

which gives a fundamental frequency comparing to the exact solution as:

$\omega_k \approx \left( \frac{3.31}{L} \right)^2 \sqrt{\frac{EI}{\mu}} \approx \left( \frac{ \pi}{L} \right)^2 \sqrt{\frac{EI}{\mu}}$

with an error of approximatelly 11%. The modal shape approximation may also be used to estimate the displacement at beam center, for which we calculate the modal mass and the modal load as:

$\begin{align*} M_k &= \int_0^L{\mu \, \left[ \frac{4x(L - x)}{L^2} \right]^2 \, dx} = \frac{8}{15}\mu L\\ F_k &= \int_0^L{ q \, \left[ \frac{4x(L - x)}{L^2} \right] \, dx} = \frac{2}{ 3} q L \end{align*}$

The modal stiffness is then:

$K_k = \frac{ 120}{L^4} \left( \frac{EI}{\mu} \right) \cdot \frac{8}{15}\mu L = \frac{64EI}{L^3}$

and the modal displacement is:

$u_k = \frac{2}{3} qL \cdot \frac{L^3}{64EI} = \frac{4 q L^4}{384EI} \approx \frac{5 q L^4}{384EI}$

The modal displacement is already the displacement at the beam center, for $\varphi(L/2) = 1$ . The implied error is hence 20%, not bad for an arbitrary elastic line shape. The following scripts show how to numerically accomplish these calculations.

In [8]:

L  =  1         # bar length (m)
EI =  2.6       # bending stiffness (Nm2)
mu =  0.260     # mass per unity length (kg/m)
q  =  mu*9.81   # distributed load is self weight (N/m)

# Proposed modal shape for first mode (second order polynomial)
x  =  np.linspace(0, L, 200)
qk =  4*x*(L - x)/L/L          # guessed modal shape
q0 =  np.sin(np.pi*x/L)        # exact modal shape!!!

plt.figure(1, figsize=(10,4))
plt.plot(x, qk, 'b', x, q0, 'r')
plt.legend(('Trial', 'Exact'))

plt.xlim( 0.0, L  );  plt.xlabel('x');
plt.ylim(-0.5, 1.5); plt.ylabel('phi(x)');
plt.title('Proposed modal shape')
plt.grid(True)

The same calculation could be carried out with the correct modal frequency, with much more accurate results:

In [9]:

wk =  ((np.pi/L)**2)*np.sqrt(EI/mu)     # exact fundamental frequency
fk =    wk/(2*np.pi)
Mk =    np.sum(mu*qk*qk)                # modal mass from guessed modal shape
Kk =    wk*wk*Mk                        # improved modal stiffness

print('Available fundamental vibration frequency: {0:5.2f} Hz'.format(fk))
print('Modal mass (integrated over bar length):   {0:5.1f} kg'.format(Mk))
print('Modal stiffness (from mass and frequency): {0:5.0f} N/m'.format(Kk))

Available fundamental vibration frequency:  4.97 Hz
Modal mass (integrated over bar length):    27.6 kg
Modal stiffness (from mass and frequency): 26880 N/m

In [10]:

Fk =  np.sum(q*qk)         # modal force
uk =  Fk/Kk                # modal displacement

wp =  np.max(uk*qk)        # approximated elastic line
w0 = (5*q*L**4)/(384*EI)   # exact elastic line

print('Maximum displacement approximation: {0:6.2f}mm'.format(1000*wp))
print('Theoretical maximum displacement:   {0:6.2f}mm'.format(1000*w0))

Maximum displacement approximation:  12.59mm
Theoretical maximum displacement:    12.77mm

In [11]:

w    =  uk*qk
V    =  np.sum( q*w  )/2    # potential energy calculated with external work
Tref =  np.sum(mu*w*w)/2

wk   =  np.sqrt(V/Tref)
fk   =  wk/(2*np.pi)

print('Fundamental frequency from Rayleigh quotient: {0:5.2f} Hz'.format(fk))

Fundamental frequency from Rayleigh quotient:  4.97 Hz

In [ ]: