
In this lecture our objective is to demonstrate equivalence between direct integration of the equations of motion and modal superposition. Now, we've already done this for free and forced vibration of an **undamped system**. However, now we'll focus on doing the same thing for a damped system, which is much more useful, since all systems have at least some damping. Since the presence of damping is more obvious in the free vibration case, we'll deal with that here.

We'll assume that we want to model the free vibration of our 3-storey structure with the same levels of damping specified in the last lecture, i.e. we want the modal damping to be $\xi_1 = 4\%$, $\xi_2 = 2\%$ and $\xi_3 = 2\%$. As we saw in the previous lecture, this requires the following damping matrix to be used in the direct integration,

```python
C = np.matrix([
  [13856.35354784, -3181.13438688, 28.00019939],
  [-3181.13438688, 15822.2691131, -4521.03912885],
  [28.00019939, -4521.03912885, 12511.79848761]])
```

### Solution by direct integration

First we will redefine the system parameters.

```python
#Define masses
m1 = 21145.85 #(kg)
m2 = 20773.5 #(kg)
m3 = 20476.5 #(kg)
```

```python
#Define stiffnesses and the stiffness matrix
k1 = 2494164.77 #(N/m)
k2 = 3226337.45 #(N/m)
k3 = 4593750.0 #(N/m)
```

We will use `ODEINT` since we've already done much of the setup work for this in section 3. This time we'll pass in the full damping matrix to the `systemEquations` function and assign its coefficients within the equations of motion.

```python
# Use ODEINT to solve the differential equations
from scipy.integrate import odeint
```

```python
def systemEquations(inits,t,params,):
"""
Defines the differential equations for the coupled system

var: vector of the state variables: var = [u1,y1,u2,y2,u3,y3]
    u1 = mass 1 position
    y1 = mass 1 velociity
    u2 = mass 2 position
    y2 = mass 2 velocity
    u3 = mass 3 position
    y3 = mass 3 velocity
t: time
params:  vector of parameters: params = [m1,m2,m3,k1,k2,k3,c1,c2,c3]
"""

# Unpack position and velocity variables and the system parameters
u1, y1, u2, y2, u3, y3 = inits # Position and velocity variables
m1, m2, m3, k1, k2, k3, C, P, w = params # System parameters

#Unpack the damping coefficients
c11 = C[0,0]
c12 = C[0,1]
c13 = C[0,2]
c21 = C[1,0]
c22 = C[1,1]
c23 = C[1,2]
c31 = C[2,0]
c32 = C[2,1]
c33 = C[2,2]

#Define the force vector (legacy)
F = F3*np.sin(w*t)

# Define a system of first order differential equations
#(introduce dummy variables y1 and y2 for this)
f = [y1,
     ( -y1*c11 - y2*c12 -(k1+k2)*u1 + k2*u2)/m1,
     y2,
     (-y1*c21 - y2*c22 - y3*c23 + u1*k2 - u2*(k2+k3) + u3*k3)/m2,
    y3,
    (F - y2*c32 - y3*c33 + u2*k3 - u3*k3)/m3]

return f
```

Note the change in some damping force signs within the equations of motion compared to when we used this function previously. Remember that the calculated damping coefficients already have the signs _baked in_.

Now we can establish some initial conditions, construct a time vector and call the ODE solver to calculate our solution.

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

```python
# Initial conditions
x1 = 0.012 #initial displacement of mass 1
y1 = 0.0 #initial velocity of mass 1
x2 = 0.019 #initial displacement of mass 2
y2 = 0.0 #initial velocity of mass 2
x3 = 0.021 #initial displacement of mass 2
y3 = 0.0 #initial velocity of mass 2

#Set the forcing term magnitude and frequency to zero
F3 = 0
w=0

# Pack up the parameters and initial conditions:
params = [m1, m2, m3, k1, k2, k3, C, F3, w]
inits = [x1, y1, x2, y2, x3, y3]

# Call the ODE solver.
sol = odeint(systemEquations, inits, t, args=(params,))

U1_directFree = sol[:,0]
U2_directFree = sol[:,2]
U3_directFree = sol[:,4]
```

Next we can visualise the response.

```python
fig, axes = plt.subplots(figsize=(12,6),nrows=1,ncols=1)
axes.plot(t,U1_directFree,'-',label='mass 1 (first floor)')
axes.plot(t,U2_directFree,'-',label='mass 2 (second floor)')
axes.plot(t,U3_directFree,'-',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/direct-integration-versus-modal-superposition-for-a-damped-system/img1.jpg"
	alt="Damped free vibration structural response of masses 1, 2 and 3 obtained from direct integration. | EngineeringSkills.com"
	width="1000"
	height="515"
	caption="Fig 1. Damped free vibration structural response of masses 1, 2 and 3 obtained from direct integration."
/>

### Solution by modal superposition

Again we'll rely heavily on what we've already accomplished in previous lectures in section four and in particular [lecture 28]("/course/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/compare-modal-superposition-with-direct-integration-free-vibration).

We can start by redefining some parameters.

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

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

Next we convert our initial displacements into 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).item()
Y0_mode_2 = (Phi[:,1].transpose()*M*u0).item()
Y0_mode_3 = (Phi[:,2].transpose()*M*u0).item()
```

We'll use the equation derived previously to calculate the free vibration response of each mode with initial position and velocity, $u_0$ and $v_0$ respectively,

$$
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{1}
$$

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)
xi1 = 0.04 #Damping ratio
wd1 = mf1*math.sqrt(1-xi1**2) #Damped natural frequency

# 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)
xi2 = 0.02 #Damping ratio
wd2 = mf2*math.sqrt(1-xi2**2) #Damped natural frequency

# 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)
xi3 = 0.02 #Damping ratio
wd3 = mf3*math.sqrt(1-xi3**2) #Damped natural frequency
```

Calculate the modal response using the equation above.

```python
Y1 = np.exp(-xi1*mf1*t)*((((Y0_mode_1*xi1*mf1)/wd1)*np.sin(wd1*t)) + (Y0_mode_1*np.cos(wd1*t)))
Y2 = np.exp(-xi2*mf2*t)*((((Y0_mode_2*xi2*mf2)/wd2)*np.sin(wd2*t)) + (Y0_mode_2*np.cos(wd2*t)))
Y3 = np.exp(-xi3*mf3*t)*((((Y0_mode_3*xi3*mf3)/wd3)*np.sin(wd3*t)) + (Y0_mode_3*np.cos(wd3*t)))
```

Next we can recombine the modal responses.

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

And now plot the structural response of 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/direct-integration-versus-modal-superposition-for-a-damped-system/img2.jpg"
	alt="Damped free vibration structural response of masses 1, 2 and 3 (top to bottom) obtained using modal superposition. | EngineeringSkills.com"
	width="1000"
	height="515"
	caption="Fig 2. Damped free vibration structural response of masses 1, 2 and 3 (top to bottom) obtained using modal superposition."
/>

At first glance all seems to be in order and we see close agreement with the direct integration solution. We can confirm this by comparing both solutions for each mass, as we've done before.

### Compare the modal response with the direct integration response

```python
fig, axes = plt.subplots(figsize=(12,14),nrows=3,ncols=1)
axes[0].plot(t,u1,'-',label='Modal superposition')
axes[0].plot(t,U1_directFree,'--',label='Direct integration')
axes[0].set_xlim([0,tmax])
axes[0].set_ylabel('Structural Displacement (m)')
axes[0].set_title('Mass 1')
axes[0].grid()
axes[0].legend(loc='lower right')

axes[1].plot(t,u2,'-',label='Modal superposition')
axes[1].plot(t,U2_directFree,'--',label='Direct integration')
axes[1].set_xlim([0,tmax])
axes[1].set_ylabel('Structural Displacement (m)')
axes[1].set_title('Mass 2')
axes[1].grid()
axes[1].legend(loc='lower right')

axes[2].plot(t,u3,'-',label='Modal superposition')
axes[2].plot(t,U3_directFree,'--',label='Direct integration')
axes[2].set_xlim([0,tmax])
axes[2].set_xlabel('Time (sec)')
axes[2].set_ylabel('Structural Displacement (m)')
axes[2].set_title('Mass 3')
axes[2].grid()
axes[2].legend(loc='lower right')

plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/direct-integration-versus-modal-superposition-for-a-damped-system/img3.jpg"
	alt="Comparison of damped free vibration structural responses for each mass obtained using direct integration and modal superposition. | EngineeringSkills.com"
	width="1000"
	height="1093"
	caption="Fig 3. Comparison of damped free vibration structural responses for each mass obtained using direct integration and modal superposition."
/>

Again as expected, both response calculations agree. At this point we have demonstrated complete equivalence of the solutions obtained using direct integration and modal superposition, for damped and undamped systems. You should now be comfortable taking either approach in your MDoF dynamic analysis. We will now move forward exclusively with a modal superposition approach for the remaining lectures in the course. Remember, this implies we are dealing with linear structural behaviour; if you encounter non-linear behaviour, you will need to deploy a direct solution approach.


Please log in or enroll to continue

Please log in or enroll to continue