
We can obtain the natural frequencies and mode shapes for a system by examining its undamped free vibration. Consider the 2DoF shear building, Fig 1. We'll scale back up to a 3DoF structure soon, but a trimmed down model will do for this demonstration.

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/finding-natural-frequencies-and-mode-shapes/img1.jpg"
	alt="2DoF shear building dynamic model without damping. | EngineeringSkills.com"
	width="500"
	height="282"
	caption="Fig 1. 2DoF shear building dynamic model without damping."
/>

Our objective is to determine the natural frequencies and mode shapes for this shear building. We start with the equations of motion,

$$
m_1\ddot{u}_1 + (k_1+k_2)u_1-k_2u_2 = 0
\tag{1}
$$

$$
m_2\ddot{u}_2 + k_2u_2 - k_2u_1 = 0
\tag{2}
$$

Next we make the assumption that the free vibration response of each mass is harmonic,

$$
u_1(t) = U_1\sin(\omega t)
\tag{3}
$$

$$
u_2(t) = U_2\sin(\omega t)
\tag{4}
$$

We differentiate these assumed solutions in order to substitute them back into the original differential equations,

$$
\ddot{u}_1(t) = - \omega^2u_1(t)
\tag{5}
$$

$$
\ddot{u}_2(t) = - \omega^2u_2(t)
\tag{6}
$$

Substituting the solutions back into the equations of motion yields,

$$
(k_1 + k_2 - m_1\omega^2)U_1 - k_2U_2=0
\tag{7}
$$

$$
(k_2-m_2\omega^2)U_2 - k_2U_1 = 0
\tag{8}
$$

We now remove the amplitude terms $U_1$ and $U_2$ by first dividing through by $U_2$,

$$
(k_1 + k_2 - m_1\omega^2)\frac{U_1}{U_2} - k_2=0
\tag{9}
$$

$$
(k_2-m_2\omega^2) - k_2\frac{U_1}{U_2} = 0
\tag{10}
$$

By substituting the preceding equation into its predecessor we obtain,

$$
(k_1 + k_2 - m_1\omega^2)(k_2-m_2\omega^2) - k_2^2=0
\tag{11}
$$

Which expands to,

$$
m_1m_2\omega^4 - (k_2m_1 + k_2m_2 + k_1m_2)\omega^2 + k_1k_2 = 0
\tag{12}
$$

This is known as the **characteristic equation** for the system and is a quadratic in $\omega^2$. When solved, this quadratic results in two values of $\omega^2$ and therefore **two natural frequencies**, $\omega_1$ and $\omega_2$. For each value of $\omega$ we can solve for the corresponding ratio of amplitudes, $U_1/U_2$ above. These relative amplitudes represent the **mode shape** or displacement patterns of the structure in each mode. Each frequency, $\omega$ has a corresponding mode shape.

### Calculating natural frequencies

Let's use the technique described above to calculate the natural frequencies for the shear building from the previous section, with the top floor removed and no damping.

```python
h1 = 6.5 #(m) Story height 0-1
h2 = 4.5 #(m) Story height 1-2
b = 5 #(m) Frame width
w = 6 #(m) Out of plane frame spacing
```

```python
gamma_concrete = 2450 #(kg) Mass of concrete per m^3
gamma_UC254 = 89 #(kg) Mass of 254UC89 per m length
gamma_UC203 = 60 #(kg) Mass of 203UC60 per m length
miscAllowance = 10 #(%) Additional mass allowance per floor for beams and misc structure
slabDepth = 0.25 #(m) Depth of concrete slab
```

```python
m1_slab = b*w*slabDepth*gamma_concrete #(kg) Mass of level 1 concrete slab
m1_cols = h1*gamma_UC254 + h2*gamma_UC203 #(kg) Mass of half of columns above and below level 1
m1_misc = (miscAllowance/100)*(m1_slab+m1_cols) #(kg) Mass of floor beams and misc structure
m1 = m1_slab + m1_cols + m1_misc #(kg) Total lumped mass for level 1

m2_slab = m1_slab #(kg) Mass of level 2 concrete slab
m2_cols = h2*gamma_UC203 # Mass of half of columns below level 2
m2_misc = (miscAllowance/100)*(m2_slab+m2_cols) #(kg) Mass of floor beams and misc structure
m2 = m2_slab + m2_cols + m2_misc #(kg) Total lumped mass for level 2
```

```python
Ix_254 = 14270*10**-8 #(m^4) Second moment of area for 254UC89
Ix_203 = 6125*10**-8 #(m^4) Second moment of area for 203UC60
E = 200*10**9 #(N/m^2) Modulus of Elasticity

k1 = 2 * (12*E*Ix_254/h1**3) #(N/m) Combined column stiffness level 0-1
k2 = 2 * (12*E*Ix_203/h2**3) #(N/m) Combined column stiffness level 1-2
```

Now we can construct the coefficients in the characteristic equation.

```python
#Coefficients in characteristic equation
coeff1 = m1*m2 #Coefficient of omega^4
coeff2 = -(k2*m1 + k2*m2 + k1*m2) #Coefficient of omega^2
coeff3 = k1*k2 #Coefficient of constant

print(f'Coefficient of omega^4 = {round(coeff1)}')
print(f'Coefficient of omega^2 = {round(coeff2)}')
print(f'Constant = {round(coeff3)}')
```

Therefore the characteristic equation works out to be,

$$
(439273\times 10^3)\omega^4 - (-187058\times 10^6)\omega^2 + 8047017\times 10^6
\tag{13}
$$

Next we solve for the roots of the quadratic.

```python
#Solving the quatratic
a = coeff1
b = coeff2
c = coeff3

omega1_squared = (-b - math.sqrt(b**2 - 4*a*c)) / (2 * a) #Root 1
omega2_squared = (-b + math.sqrt(b**2 - 4*a*c)) / (2 * a) #Root 2

omega1 = math.sqrt(omega1_squared) #(rads/sec)First angular natural frequency
omega2 = math.sqrt(omega2_squared) #(rads/sec) Second angular natural frequency
f1 = omega1/(2*pi) #(Hz)
f2 = omega2/(2*pi) #(Hz)

print(f'The first natural frequency is {round(omega1,1)} rads/s or {round(f1,1)} Hz')
print(f'The second natural frequency is {round(omega2,1)} rads/s or {round(f2,1)} Hz')

```

- `The first natural frequency is $7.0 \:rads/s$ or $1.1 \:Hz$`
- `The second natural frequency is $19.4 \:rads/s$ or $3.1 \:Hz$`

### Calculating the mode shapes

Now that we have the natural frequencies, we can determine their corresponding mode shapes by substituting $\omega_1$ and $\omega_2$ into the equations in our derivation above containing the ratio $U1/U2$. In this case, we use equation 10,

```python
U1U2_mode1 = (1/k2)*(k2-m2*omega1_squared)
U1U2_mode2 = (1/k2)*(k2-m2*omega2_squared)

print(f'The ratio U1/U2 for omega of {round(omega1,1)} rads/s is {round(U1U2_mode1,3)}')
print(f'The ratio U1/U2 for omega of {round(omega2,1)} rads/s is {round(U1U2_mode2,3)}')
```

- `The ratio $U1/U2$ for omega of $7.0 \:rads/s$ is $0.687$`
- `The ratio $U1/U2$ for omega of $19.4 \:rads/s$ is $-1.429$`

It is typical to assign a value of one to one of the amplitudes. In this case we will assign a value of $1$ to $U_1$ in each case and scale the other value appropriately.

```python
U1_omega1 = 1
U2_omega1 = 1/U1U2_mode1

U1_omega2 = 1
U2_omega2 = 1/U1U2_mode2

print(f'For {round(omega1,1)} rads/s, the mode shape is by U1/U2 of {round(U1_omega1,2)}/{round(U2_omega1,2)}')
print(f'For {round(omega2,1)} rads/s, the mode shape is by U1/U2 of {round(U1_omega2,2)}/{round(U2_omega2,2)}')
```

- For omega of $7.0 \:rads/s$, the mode shape is described by $U1/U2$ of $1/1.45$
- For omega of $19.4 \:rads/s$, the mode shape is described by $U1/U2$ of $1/-0.7$

### Plotting the mode shapes

It's helpful to visualise the mode shapes with a simple plot.

```python
#Plotting
fig, axes = plt.subplots(figsize=(12,8),nrows=1,ncols=2)

axes[0].plot([0,0,0],[0, h1, h1+h2],'-k')
axes[0].plot([0,U1_omega1, U2_omega1],[0, h1, h1+h2])
axes[0].plot([U1_omega1],[h1],'-o',markerfacecolor='r', markeredgecolor='k', markersize=13)
axes[0].plot([U2_omega1],[h1+h2],'-o',markerfacecolor='r', markeredgecolor='k', markersize=13)
axes[0].set_xlabel('Relative amplitude')
axes[0].set_ylabel('Height (m)')
axes[0].set_title('Mode shape for $\omega_1$ =  {one} rads/sec'.format(one=round(omega1,1)))
axes[0].set_xlim([-1.6,1.6])
axes[0].grid()

axes[1].plot([0,0,0],[0, h1, h1+h2],'-k')
axes[1].plot([0,U1_omega2, U2_omega2],[0, h1, h1+h2])
axes[1].plot([U1_omega2],[h1],'-o',markerfacecolor='r', markeredgecolor='k', markersize=13)
axes[1].plot([U2_omega2],[h1+h2],'-o',markerfacecolor='r', markeredgecolor='k', markersize=13)
axes[1].set_xlabel('Relative amplitude')
axes[1].set_ylabel('Height (m)')
axes[1].set_title('Mode shape for $\omega_2$ =  {one} rads/sec'.format(one=round(omega2,1)))
axes[1].set_xlim([-1.6,1.6])
axes[1].grid()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/finding-natural-frequencies-and-mode-shapes/img2.jpg"
	alt="Mode shapes for mode 1 (left) and mode 2 (right). | EngineeringSkills.com"
	width="500"
	height="282"
	caption="Fig 2. Mode shapes for mode 1 (left) and mode 2 (right)."
/>

### Simulating modal oscillations

At this point it might be helpful to animate the modal oscillations. To do this we first need to generate an oscillation history for each mass in each mode. This is the modal oscillation of each mass we'll play back in our animation.

```python
# Generate a time vector
tmax = 20 #(sec) The max time
delta_t = 0.01 #(sec) The timestep
nPoints = tmax/delta_t #Number of equally spaced data points
t = np.linspace(0,tmax, int(nPoints)) # Time vector
```

```python
#Generate an oscillation for each mass in each mode

#Mode 1 oscillation
um1m1 = 1*U1_omega1*np.sin(omega1*t) #Motion of mass 1, mode 1
um2m1 = 1*U2_omega1*np.sin(omega1*t) #Motion of mass 2, mode 1

#Mode 2 oscillation
um1m2 = 1*U1_omega2*np.sin(omega2*t) #Motion of mass 1, mode 2
um2m2 = 1*U2_omega2*np.sin(omega2*t) #Motion of mass 2, mode 2
```

Next we import the relevant animation functionality and change the plotting mode for the notebook (we need to do this for animations).

```python
from matplotlib.animation import FuncAnimation
%matplotlib notebook
```

```python
#Animation parameters
delta_t = 0.01
animLength = 10 #(sec)
frameRate = 20 #(5,10,20) frames per second (too high and animation slows down)
plotInterval = 1/frameRate #(sec) time between frame plots
dataInterval = plotInterval/delta_t #Plot moving elements every 'dataInterval-th' point

#Set up plots
fig, (ax1, ax2) = plt.subplots(nrows=1,ncols=2, figsize=(10, 5))
ax1.set_xlabel('Relative amplitude')
ax1.set_ylabel('Height (m)')
ax1.set_title('Mode shape for $\omega_1$ =  {one} rads/sec'.format(one=round(omega1,1)))
ax1.set_xlim([-1.6,1.6])
ax2.set_xlabel('Relative amplitude')
ax2.set_ylabel('Height (m)')
ax2.set_title('Mode shape for $\omega_2$ =  {one} rads/sec'.format(one=round(omega2,1)))
ax2.set_xlim([-1.6,1.6])
ax1.grid()
ax2.grid()

#Plot initial state
line1, = ax1.plot([0,0,0],[0, h1, h1+h2],'-',lw=2)
line2, = ax2.plot([0,0,0],[0, h1, h1+h2],'-',lw=2)
dotMass1Mode1, = ax1.plot([0], h1, 'ro', markeredgecolor='k', markersize=15)
dotMass2Mode1, = ax1.plot([0], h1+h2, 'ro', markeredgecolor='k', markersize=15)
dotMass1Mode2, = ax2.plot([0], h1, 'ro', markeredgecolor='k', markersize=15)
dotMass2Mode2, = ax2.plot([0], h1+h2, 'ro', markeredgecolor='k', markersize=15)

#Function to animate plot objects
def animate(i):
    frm = int(i*dataInterval)

    line1.set_data([0,um1m1[frm],um2m1[frm]], [0, h1, h1+h2])
    dotMass1Mode1.set_data(um1m1[frm],h1)
    dotMass2Mode1.set_data(um2m1[frm],h1+h2)

    line2.set_data([0,um1m2[frm],um2m2[frm]], [0, h1, h1+h2])
    dotMass1Mode2.set_data(um1m2[frm],h1)
    dotMass2Mode2.set_data(um2m2[frm],h1+h2)

    return line1, line2, dotMass1Mode1, dotMass2Mode1, dotMass1Mode2, dotMass2Mode2

#Function to generate the animation
myAnimation = FuncAnimation(fig,
                            animate,
                            frames=int(1 + (animLength/plotInterval) ),
                            interval=plotInterval*1000,  #milliseconds
                            blit=True,
                            repeat=True)

plt.show()
myAnimation.save('2Modes.gif')
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/finding-natural-frequencies-and-mode-shapes/img3.gif"
	alt="Animated mode shapes for mode 1 (left) and mode 2 (right). | EngineeringSkills.com"
	width="1000"
	height="500"
	caption="Fig 3. Animated mode shapes for mode 1 (left) and mode 2 (right)."
/>

When we animate the mode shapes like this, we can more easily intuit their relative contribution to the overall dynamic behaviour of the structure. Having now observed the animation, recall the discussion from the previous lecture in which we considered the modal response of a simple beam.

What we said in that discussion about:

- the proximity of the forcing frequency to the modal frequency and,
- the physical position of the dynamic load on the structure,

govern the response in a mode. Depending on the location and frequency of dynamic loading, some modes will be excited more than others. That discussion maps directly to the vibration modes in this structure. Therefore dynamic force with a frequency of $7\: rads/s$ applied to the top floor will have the maximum dynamic influence on mode $1$. Whereas to have the maximum influence on mode $2$, we should apply a dynamic force with a frequency of $19.5\: rads/s$ to the first floor. Remember that it is the combined response of both modes that gives us the overall dynamic response of the structure.


Please log in or enroll to continue

Please log in or enroll to continue