
In this lecture we will consider how to handle the free vibration of our NDoF system using a modal approach. The only additional element we need to consider in this case is how to transfer the system initial conditions into modal coordinates. We start with the statement of modal superposition we derived previously which relates geometric displacements with modal displacements.

$$
\begin{Bmatrix}
u
\end{Bmatrix}  = \sum_{k=1}^n
\begin{Bmatrix}
\phi
\end{Bmatrix}_k \:Y_k
\tag{1}
$$

where $Y_k$ is the response in mode $k$ and $\begin{Bmatrix}\phi_k\end{Bmatrix}$ is the mass normalised eigenvector for mode $k$.

Expanding the right hand side of this equation yields,

$$
\begin{Bmatrix}
u
\end{Bmatrix}=
Y_1
\begin{Bmatrix}
\phi
\end{Bmatrix}_1 +
Y_2
\begin{Bmatrix}
\phi
\end{Bmatrix}_2 +...+
Y_n
\begin{Bmatrix}
\phi
\end{Bmatrix}_n
\tag{2}
$$

In order to evaluate a modal coordinate $Y_N$, given an initial set of geometric displacements $\begin{Bmatrix}u\end{Bmatrix}$, we first pre-multiply both sides by $\begin{Bmatrix}\phi\end{Bmatrix}_n^T\begin{bmatrix}M\end{bmatrix}$, giving,

$$
\begin{Bmatrix}
\phi
\end{Bmatrix}_n^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
u
\end{Bmatrix}=
\begin{Bmatrix}
\phi
\end{Bmatrix}_n^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
\phi
\end{Bmatrix}_1 Y_1+
\begin{Bmatrix}
\phi
\end{Bmatrix}_n^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
\phi
\end{Bmatrix}_2 Y_2+...+
\begin{Bmatrix}
\phi
\end{Bmatrix}_n^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
\phi
\end{Bmatrix}_n Y_n
\tag{3}
$$

Because the mode shape vectors are orthogonal with respect to mass as indicated previously, all terms on the right hand side evaluate to zero except for the term containing $\begin{Bmatrix}\phi\end{Bmatrix}_n^T\begin{bmatrix}M\end{bmatrix}\begin{Bmatrix}\phi\end{Bmatrix}_n$. Therefore simplifying and rearranging we obtain our expression for the modal coordinate $Y_n$ as a function of the geometric nodal displacement vector $\begin{Bmatrix}u\end{Bmatrix}$,

$$
Y_n = \frac{\begin{Bmatrix}
\phi
\end{Bmatrix}_n^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
u
\end{Bmatrix}}{\begin{Bmatrix}
\phi
\end{Bmatrix}_n^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
\phi
\end{Bmatrix}_n }
\tag{4}
$$

Note that for mass normalised mode shapes, the denominator evaluates to 1. The same expression holds for the derivative of $\dot{Y}_n$ , i.e. a velocity initial condition can be transferred to modal coordinates using the expression,

$$
\dot{Y}_n = \frac{\begin{Bmatrix}
\phi
\end{Bmatrix}_n^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
\dot{u}
\end{Bmatrix}}{\begin{Bmatrix}
\phi
\end{Bmatrix}_n^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{Bmatrix}
\phi
\end{Bmatrix}_n }
\tag{5}
$$

### Calculating initial conditions in modal coordinates

Let's assume that the structure is released from a stationary displacement defined by the following initial conditions,

$$
u(0) =
\begin{Bmatrix}
0.012\\
0.019\\
0.021
\end{Bmatrix}\text{m}
\tag{6}
$$

$$
\dot{u}(0) =
\begin{Bmatrix}
0\\
0\\
0
\end{Bmatrix}\text{m/s}
\tag{7}
$$

We need to transform these initial displacements into modal coordinates using the equations defined above. But first we need to redefine some of the structural properties calculated previously.

```python
#Define masses and the mass matrix
m1 = 21145.85
m2 = 20773.5
m3 = 20476.5
M = np.matrix([[m1,0,0],
               [0, m2,0],
               [0,0,m3]])
```

```python
#Define stiffnesses and the stiffness matrix
k1 = 2494164.77
k2 = 3226337.45
k3 = 4593750.0
K = np.matrix([[k1+k2, -k2, 0],
               [-k2, k2+k3, -k3],
               [0, -k3, k3]])
```

```python
#Define the mass normalised modal matrix
Phi = np.array(
[[ 0.00265937,  0.00577103,  0.00262936],
 [ 0.00423975,  0.00053654, -0.00546579],
 [ 0.00482667, -0.00376177,  0.00337473]])
```

```python
#Define the modal stiffness matrix
modalK = np.array([[27.28, 0, 0],
                   [0, 256.341, 0],
                   [0, 0, 587.694]])
```

Now we can transfer our initial conditions to modal coordinates.

```python
#Construct a vector of initial conditions
u0 = np.matrix([0.012, 0.019, 0.021]).reshape([3,1])

# Transfer initial conditions to modal coordinates
Y0_mode_1 = Phi[:,0].transpose()*M*u0
Y0_mode_2 = Phi[:,1].transpose()*M*u0
Y0_mode_3 = Phi[:,2].transpose()*M*u0

print(f'The initial modal displacement for mode 1 is {round(Y0_mode_1[0,0],3)}')
print(f'The initial modal displacement for mode 2 is {round(Y0_mode_2[0,0],3)}')
print(f'The initial modal displacement for mode 3 is {round(Y0_mode_3[0,0],3)}')
```

- `The initial modal displacement for mode 1 is 4.424`

- `The initial modal displacement for mode 2 is 0.059`

- `The initial modal displacement for mode 3 is -0.039`

### Calculating the free vibration modal response

The free vibration modal response can now be determined using standard closed-form solution. Although our structure doesn't yet include damping, for generality we'll base the following on the damped free vibration response.

$$
u(t) = e^{-\xi\omega_nt}[A\sin(\omega_dt) + B\cos(\omega_dt)]
\tag{8}
$$

where, $\omega_d$ is defined as the damped angular natural frequency and given by,

$$
\omega_d = \sqrt{\frac{k}{m} - \left(\frac{c}{2m}\right)^2}
\tag{9}
$$

To determine the unknown constants $A$ and $B$, we must use the system initial conditions. Let's assume that the position and velocity of the mass at time $t=0$ can be denoted by $u_0$ and $v_0$. This way we retain the variables and can substitute actual values in later for $u_0$ and $v_0$.

#### Initial condition #1

Replacing the left side of the response equation above with $u_0$ and substituting $t=0$ into the right side yields,

$$
B = u_0
\tag{10}
$$

#### Initial condition #2

To make use of the second initial condition, we need to first differentiate our expression for the displacement to get the corresponding expression for velocity. First we can rewrite the response equation,

$$
u(t) = A\underbrace{\left[e^{-\xi\omega_nt}\:\sin(\omega_d t) \right]}_{f_1(t)} + B\underbrace{\left[e^{-\xi\omega_nt}\:\cos(\omega_d t) \right]}_{f_2(t)}
\tag{11}
$$

We've divided the equation into two functions so we can process the differentiation in two parts,

$$
\frac{\mathrm{d}}{\mathrm{d}t}f_1(t) = e^{-\xi\omega_nt}\omega_d\cos(\omega_dt) - \sin(\omega_dt)\: \xi\omega_n \: e^{-\xi\omega_nt}
\tag{12}
$$

$$
\frac{\mathrm{d}}{\mathrm{d}t}f_2(t) = -e^{-\xi\omega_nt}\omega_d\sin(\omega_dt) - \cos(\omega_dt)\:\xi\omega_n\:e^{-\xi\omega_nt}
\tag{13}
$$

Recombining these two expressions we obtain our expression for the velocity, $\dot{u}(t)$,

$$
\begin{split}
\dot{u}(t) =  & \:A\left[e^{-\xi\omega_nt}\omega_d\cos(\omega_dt) - \sin(\omega_dt)\: \xi\omega_n \: e^{-\xi\omega_nt} \right] + \\  & \:B\left[ -e^{-\xi\omega_nt}\omega_d\sin(\omega_dt) - \cos(\omega_dt)\:\xi\omega_n\:e^{-\xi\omega_nt}\right] + \\
\end{split}
\tag{14}
$$

Now we can replace the left side of the equation with $v_0$ and substitute $t=0$ into the right side. When we also substitute in the expression we've already found for $B$ and rearrange for $A$ we finally get,

$$
A = \frac{v_0 + u_0\xi\omega_n}{\omega_d}
\tag{15}
$$

Now we can substitute these expressions for $A$ and $B$ back into the original solution to reveal the final expression for the damped free vibration as a function of the initial conditions, $u_0$ and $v_0$ oscillation,

$$
\boxed{u(t) = e^{-\xi\omega_n t}\left[\frac{v_0 + u_0\xi \omega_n}{\omega_d}\sin(\omega_d t) + u_0\cos(\omega_d t) \right]}
\tag{16}
$$

Using this expression, we can model the damped free vibration for each mode.

### Plotting the modal responses

Establish a time vector.

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

Declare the modal properties for each SDoF system.

```python
# Mode 1 SDOF parameters
mm1 = 1 # Modal mass
mk1 = modalK[0,0] # Modal stiffness
mf1 = math.sqrt(mk1/mm1) # Modal frequency (also obtained previosuly as eigenvalue)
wd1 = mf1 #Damped natural frequency (equal to natural frequency)
xi1 = 0 #Damping ratio

# Mode 2 SDOF parameters
mm2 = 1 # Modal mass
mk2 = modalK[1,1] # Modal stiffness
mf2 = math.sqrt(mk2/mm2) # Modal frequency (also obtained previosuly as eigenvalue)
wd2 = mf2 #Damped natural frequency (equal to natural frequency)
xi2 = 0 #Damping ratio

# Mode 3 SDOF parameters
mm3 = 1 # Modal mass
mk3 = modalK[2,2] # Modal stiffness
mf3 = math.sqrt(mk3/mm3) # Modal frequency (also obtained previosuly as eigenvalue)
wd3 = mf3 #Damped natural frequency (equal to natural frequency)
xi3 = 0 #Damping ratio
```

Calculate the modal response using the equation above and note that no damping is present.

```python
Y1 = Y0_mode_1.item()*np.cos(wd1*t)
Y2 = Y0_mode_2.item()*np.cos(wd2*t)
Y3 = Y0_mode_3.item()*np.cos(wd3*t)
```

Plot the modal responses.

```python
fig, axes = plt.subplots(figsize=(12,10),nrows=3,ncols=1)
axes[0].plot(t,Y1,label='Mode 1 free vibration response')
axes[0].set_xlim([0,20])
axes[0].set_xlabel('Time (sec)')
axes[0].set_ylabel('Modal displacement')
axes[0].set_title('Mode 1 free vibration response')
axes[0].grid()
axes[0].legend(loc='lower right')

axes[1].plot(t,Y2,label='Mode 2 free vibration response')
axes[1].set_xlim([0,20])
axes[1].set_xlabel('Time (sec)')
axes[1].set_ylabel('Modal displacement')
axes[1].set_title('Mode 2 free vibration response')
axes[1].grid()
axes[1].legend(loc='lower right')

axes[2].plot(t,Y3,label='Mode 3 free vibration response')
axes[2].set_xlim([0,20])
axes[2].set_xlabel('Time (sec)')
axes[2].set_ylabel('Modal displacement')
axes[2].set_title('Mode 3 free vibration response')
axes[2].grid()
axes[2].legend(loc='lower right')
plt.tight_layout()
plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/modal-superposition-free-vibration/img1.jpg"
	alt="Free vibration modal response in modes 1, 2 and 3 (top to bottom). | EngineeringSkills.com"
	width="1000"
	height="832"
	caption="Fig 1. Free vibration modal response in modes 1, 2 and 3 (top to bottom)."
/>

Note that since each mode is described by a freely oscillating system, the modal responses are at the natural frequency of the mode. Now we can recombine the modal responses to obtain the structural response for each mass.

```python
# Recombining modal responses to obtain the structural response for each mass
u1 = Phi[0,0]*Y1 + Phi[0,1]*Y2 + Phi[0,2]*Y3
u2 = Phi[1,0]*Y1 + Phi[1,1]*Y2 + Phi[1,2]*Y3
u3 = Phi[2,0]*Y1 + Phi[2,1]*Y2 + Phi[2,2]*Y3
```

Next, plot the structural response.

```python
fig, axes = plt.subplots(figsize=(12,6),nrows=1,ncols=1)
axes.plot(t,u1,'-',label='mass 1 (first floor)')
axes.plot(t,u2,'-',label='mass 2 (second floor)')
axes.plot(t,u3,'-',label='mass 3 (third floor)')
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Displacement (m)')
axes.set_title('Structural response')
axes.set_xlim([0,20])
axes.grid()
axes.legend(loc='lower right')
plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/modal-superposition-free-vibration/img2.jpg"
	alt="Free vibration structural response of masses 1, 2 and 3 (top to bottom). | EngineeringSkills.com"
	width="1000"
	height="519"
	caption="Fig 2. Free vibration structural response of masses 1, 2 and 3 (top to bottom)."
/>

Mode 1 clearly dominates the free vibration response. Note again that the lack of damping makes this a very unrealistic response. We would expect this response to decay exponentially to zero.


Please log in or enroll to continue

Please log in or enroll to continue