
In this lecture we'll simulate the dynamic response of our 3-storey shear building using modal superposition, Fig. 1. Recall that we said this structure was subject to a harmonic load applied to the top floor,

$$
p_3(t) = F\sin(\omega t)
$$

where $F = 5000\:N$ and $\omega=2\pi$. Note that when we simulated the motion of this building previously, we applied arbitrary levels of damping. For the purposes of this modal analysis we won't include damping as we have yet to discuss this in detail (next section).

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/modal-superposition-forced-vibration/img1.jpg"
	alt="Typical 2D plane frame (left), shear building dynamic model that approximates the dynamic behaviour of a typical frame within the structure. | EngineeringSkills.com"
	width="500"
	height="369"
	caption="Fig 1. Typical 2D plane frame (left), shear building dynamic model that approximates the dynamic behaviour of a typical frame within the structure."
/>

### Define parameters and modal matrix

We can start by defining the stiffness matrix, damping matrix and mass normalised modal matrix, all previously calculated.

```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
F3 = 5000 #(N) Force amplitude
w = 2*math.pi #(rads/s) Force frequency
```

### Decoupling the equations of motion

Recall that we decouple the equations using the modal matrix by evaluating the modal mass, stiffness and force from the equation,

$$
\underbrace{
\begin{bmatrix}
\Phi
\end{bmatrix}^T
\begin{bmatrix}
M
\end{bmatrix}
\begin{bmatrix}
\Phi
\end{bmatrix}}_{\text{diagonal modal mass}}
\begin{Bmatrix}
\ddot{Y}
\end{Bmatrix} +
\underbrace{
\begin{bmatrix}
\Phi
\end{bmatrix}^T
\begin{bmatrix}
K
\end{bmatrix}
\begin{bmatrix}
\Phi
\end{bmatrix}}_{\text{diagonal modal stiffness}}
\begin{Bmatrix}
Y
\end{Bmatrix} =
\underbrace{
\begin{bmatrix}
\Phi
\end{bmatrix}^T
\begin{Bmatrix}
P
\end{Bmatrix}}_{\text{modal force vector}}
\tag{1}
$$

```python
modalM = np.matrix(Phi.transpose()*M*Phi) # Modal mass matrix
modalK = np.matrix(Phi.transpose()*K*Phi) # Modal stiffness matrix
print('The modal mass matrix is:')
print(np.matrix.round(modalM,3))
print('The modal stiffness matrix is:')
print(np.matrix.round(modalK,3))
```

The modal mass matrix is:

$$
\begin{bmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1\\
\end{bmatrix}
$$

The modal stiffness matrix is:

$$
\begin{bmatrix}
27.28 & 0 & 0\\
0 & 256.341 & 0\\
0 & 0 & 587.694\\
\end{bmatrix}
$$

In order to calculate a modal force, we first construct the vector of force amplitudes and then transform this using the modal matrix.

```python
P = np.mat([0, 0, 5000]).reshape([3,1])
modalP = np.matrix(Phi.transpose()*P) # Modal force vector

print('The modal force vector is:')
print(np.matrix.round(modalP,3))
```

The modal force vector is:

$$
\begin{Bmatrix}
24.133\\
-18.809\\
16.874\\
\end{Bmatrix}
$$

Therefore, our coupled dynamic system can be considered as the summation of the following 3 independent SDoF systems...

SDoF system 1, independent differential equation describing the response (in modal coordinates) of mode 1:

$$
\ddot{Y}_1(t) + 27.28 Y_1(t) = 24.133 \sin(2\pi t)
\tag{2}
$$

SDoF system 2, independent differential equation describing the response (in modal coordinates) of mode 2:

$$
\ddot{Y}_2(t) + 256.341 Y_2(t) = -18.809 \sin(2\pi t)
\tag{3}
$$

SDoF system 3, independent differential equation describing the response (in modal coordinates) of mode 3:

$$
\ddot{Y}_3(t) + 587.694 Y_3(t) = 16.874 \sin(2\pi t)
\tag{4}
$$

Now we can employ any method we want to solve each independent SDoF system. In this case the excitation is harmonic and so a closed-form analytical solution can be used. However, if the excitation force was an effective ground motion force derived from a ground acceleration, we could just as easily employ a numerical solution technique like the piecewise exact method (we'll see this later).

The [steady-state response to harmonic excitation](/course/fundamentals-of-structural-dynamics-with-python/characterising-the-steady-state-response) was derived in the previous Fundamentals of Structural Dynamics course and is,

$$
u(t) = \frac{P_0}{k}\left[\frac{1}{(1-\beta^2)^2+(2\xi\beta)^2}\right][(1-\beta^2)\sin(\omega t) - 2\xi\beta\cos(\omega t)]
\tag{5}
$$

where $\beta = \omega/\omega_n$ and $\xi$ is the damping ratio. Using this equation would be fine if damping was present in the system, because the transient component of the response would quickly decay to zero leaving only the stead-state response.

However, since in this case there is no damping, the transient behaviour never decays to zero and we never reach a constant amplitude steady-state response i.e. the resulting motion is not harmonic. So, we need to employ the closed-form solution for the case with zero damping. It can be readily confirmed that the solution is given by,

$$
u(t) = \frac{P_0}{k}\left[\frac{\sin(\omega t) - \beta \sin(\omega_n t)}{1-\beta^2}\right]
\tag{6}
$$

### Calculating the modal responses

First we 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 = modalM[0,0] # Modal mass
mk1 = modalK[0,0] # Modal stiffness
mf1 = math.sqrt(mk1/mm1) # Modal frequency (also obtained previosuly as eigenvalue)
beta1 = w/mf1 #Frequency ratio for mode 1
mp1 = modalP[0,0] # Modal Force

# Mode 2 SDOF parameters
mm2 = modalM[1,1] # Modal mass
mk2 = modalK[1,1] # Modal stiffness
mf2 = math.sqrt(mk2/mm2) # Modal frequency (also obtained previosuly as eigenvalue)
beta2 = w/mf2 #Frequency ratio for mode 2
mp2 = modalP[1,0] # Modal Force

# Mode 3 SDOF parameters
mm3 = modalM[2,2] # Modal mass
mk3 = modalK[2,2] # Modal stiffness
mf3 = math.sqrt(mk3/mm3) # Modal frequency (also obtained previosuly as eigenvalue)
beta3 = w/mf3 #Frequency ratio for mode 3
mp3 = modalP[2,0] # Modal Force
```

Calculate the modal response using the equation above.

```python
Y1 = ((mp1/mk1)/(1-beta1**2))*(np.sin(w*t) - beta1*np.sin(mf1*t))
Y2 = ((mp2/mk2)/(1-beta2**2))*(np.sin(w*t) - beta2*np.sin(mf2*t))
Y3 = ((mp3/mk3)/(1-beta3**2))*(np.sin(w*t) - beta3*np.sin(mf3*t))
```

Plot the modal responses.

```python
fig, axes = plt.subplots(figsize=(12,6),nrows=1,ncols=1)
axes.plot(t,Y1,'-',label='mode 1 response, $\omega = 5.22$ rads/sec')
axes.plot(t,Y2,'-',label='mode 2 response, $\omega = 16.01$ rads/s')
axes.plot(t,Y3,'-',label='mode 3 response, $\omega = 24.24$ rads/s')
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Modal displacement')
axes.set_title('Modal 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-forced-vibration/img2.jpg"
	alt="Response in modes 1, 2 and 3. | EngineeringSkills.com"
	width="1000"
	height="531"
	caption="Fig 2. Response in modes 1, 2 and 3."
/>

Note the relative amplitudes; the response in mode 1 with a frequency ratio of $\beta = 1.2$ is significantly larger than the response in modes 2 and 3 with frequency ratios of $0.39$ and $0.26$ respectively. From what we know of the dynamic magnification factor for SDoF systems, this is to be expected.

Now we can combine these modal responses into a single record of structural response using the modal superposition process,

$$
\begin{Bmatrix}
u_1\\
u_2\\
u_3
\end{Bmatrix}=
\begin{bmatrix}
\phi_{1,1} & \phi_{1,2} & \phi_{1,3}\\
\phi_{2,1} & \phi_{2,2} & \phi_{2,3}\\
\phi_{3,1} & \phi_{3,2} & \phi_{3,3}\\
\end{bmatrix}
\begin{Bmatrix}
Y_1\\
Y_2\\
Y_3
\end{Bmatrix}
\tag{7}
$$

Notice how the structural response of each mass contains contributions from each mode of vibration. Based on the previous plot, we can expect mode 1 to make the largest contribution.

```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
```

Finally we can plot the structural response for each mass.

```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-forced-vibration/img3.jpg"
	alt="Structural response of mass (floor) 1, 2 and 3. | EngineeringSkills.com"
	width="1000"
	height="519"
	caption="Fig 3. Structural response of mass (floor) 1, 2 and 3."
/>

Again note that the structural response does not reach a steady-state due to the lack of damping in the system. This is an unrealistic scenario as all systems will exhibit some form of damping. We'll address this in the next section of the course. In the next lecture we'll show that this response is equivalent to the response obtained by directly integrating the coupled equations of motion with zero damping.


Please log in or enroll to continue

Please log in or enroll to continue