
Two lectures ago we introduced the three-storey shear building in Fig. 1 and said that it was subject to a harmonic force $p_3(t)=F\sin(\omega t)$ applied to the third floor. Since then we've derived the coupled equations of motion and worked out how to solve them using the central difference method of numerical integration. In this lecture we'll finally solve the system of equations.

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/solving-the-coupled-eom-by-direct-integration/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."
/>

We can start by specifying the basic geometry of the structure. In addition to the dimensions identified in Fig. 1 we will also specify an out of plane frame spacing, $w$.

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

We'll keep our code as flexible as possible and calculate the masses based on material unit weights and the dimensions specified above. This way we can edit the input parameters as we like and easily observe the change in dynamic behaviour.

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

We'll combine the mass of the concrete floor slab, the supporting steel beams and half of the column mass above and below each floor, into a single lumped mass at each level. The calculation is a little tedious but it's worth it for the flexibility it provides to make updates later.

```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 + h3*gamma_UC203 #(kg) Mass of half of columns above and 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

m3_slab = m1_slab #(kg) Mass of level 3 concrete slab
m3_cols = h3*gamma_UC203 #(kg) Mass of half of columns below level 3
m3_misc = (miscAllowance/100)*(m3_slab+m3_cols) #(kg) Mass of floor beams and misc  structure
m3 = m3_slab + m3_cols + m3_misc #(kg) Total lumped mass for level 3

# Mass matrix
M = np.matrix([[m1,0,0],[0, m2,0],[0,0,m3]])
```

Next we can calculate the column stiffness. All columns are modelled as being fixed against rotation at each end therefore the stiffness of one column is,

$$
k_c = \frac{12EI}{L^3}
\tag{1}
$$

Since two columns contribute to the stiffness of each story, the story stiffness $i$ is $k_i = 2k_c$. We will also assume the columns are bending about their major principal axes.

```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
k3 = 2 * (12*E*Ix_203/h3**3) #(N/m) Combined column stiffness level 2-3

#Stiffness matrix
K = np.matrix([[k1+k2, -k2, 0],[-k2, k2+k3, -k3],[0, -k3, k3]])
```

I've eluded several times now to the fact that we're going to examine damping in more detail later. I'm afraid the suspense continues. For the purpose of this analysis, we'll specify arbitrary damping coefficients. Later we'll see how to specify specific damping coefficients that correlate with more intuitive levels of damping; but there is more ground to cover before we get there.

```python
c1 = 50000 #(Ns/m) Damping coefficient
c2 = 50000 #(Ns/m) Damping coefficient
c3 = 50000 #(Ns/m) Damping coefficient

#Damping matrix
C = np.matrix([[c1+c2, -c2, 0],[-c2, c2+c3, -c3],[0, -c3, c3]])
```

Finally we can specify the external force magnitude and frequency.

```python
F3 = 5000 #(N) External force magnitude
w = 2*math.pi #(rads/s) External force angular frequency
```

### Numerical integration

Start by establishing a time vector from $t_0$ to $t_{\text{max}}$ with time-step $h$.

```python
t0 = 0 #(s)
tmax = 20 #(s)
h=0.005 #(s)
t = np.linspace(t0,tmax, int(tmax/h)) #Time vector
t = np.round(t,2) #Tidy up trailing decimal places
```

```python
#Initialise containers to hold the displacement values for each mass
U1 = np.zeros(len(t))
U2 = np.zeros(len(t))
U3 = np.zeros(len(t))
```

Since the system starts from rest, the initial position, velocity and acceleration for each mass is zero.

```python
#Initial positions
u_n =  np.mat([0, 0, 0]).reshape([3,1]) #Position vector at time point n
u_nm1 = np.mat([0, 0, 0]).reshape([3,1])#Position vector at time point n-1
```

Now we can implement the central difference scheme explained in the previous lecture in a time-stepping loop. Note that for clarity within the code, terms in the central difference equation have been relabelled as shown below,

$$
\underbrace{\left(\frac{1}{h^2}\mathbf{M} + \frac{1}{2h}\mathbf{C} \right)}_{\mathbf{x}} \mathbf{u}_{n+1} = \mathbf{p}_n + \underbrace{\left(-\mathbf{K}+\frac{2}{h^2}\mathbf{M} \right)}_{\mathbf{y}}\mathbf{u}_n + \underbrace{\left(\frac{1}{2h}\mathbf{C} - \frac{1}{h^2}\mathbf{M} \right)}_{\mathbf{z}}\mathbf{u}_{n-1}
\tag{2}
$$

So we have,

$$
\mathbf{x}\:\mathbf{u}_{n+1} = \mathbf{p}_n +\mathbf{y}\:\mathbf{u}_n + \mathbf{z}\:\mathbf{u}_{n-1}
\tag{3}
$$

allowing us to solve for $\mathbf{u}_{n+1}$ using,

$$
\mathbf{u}_{n+1} = \mathbf{x}^{-1}(\mathbf{p}_n +\mathbf{y}\:\mathbf{u}_n + \mathbf{z}\:\mathbf{u}_{n-1})
\tag{4}
$$

Now we can step through the time window using a `for` loop.

```python
for n, ti in enumerate(t):

  #Construct the force vector
  pMag = F3*np.sin(w*ti) #Force magnitude applied to mass 3 at time point n
  p_n = np.mat([0, 0, pMag]).reshape([3,1]) #Force vector for time point n

  #Central difference formula for position at time point n+1
  x = (1/h**2)*M + (1/(2*h))*C
  y = -K + (2/h**2)*M
  z = (1/(2*h))*C -(1/h**2)*M
  u_np1 = x.I*(p_n + y*u_n + z*u_nm1)

  #Save the position data
  if(n+1)<len(t):
      U1[n+1] = u_np1[0]
      U2[n+1] = u_np1[1]
      U3[n+1] = u_np1[2]

  #Update u_n and u_nm1 for next iteration
  u_nm1 = u_n
  u_n = u_np1

print('Finished numerical integration')
```

Now let's plot the response of each mass.

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1.5])
axes.plot(t,U1,label='Mass/floor 1')
axes.plot(t,U2,label='Mass/floor 2')
axes.plot(t,U3,label='Mass/floor 3')
axes.set_xlim([0,tmax])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Structural Displacement (m)')
axes.set_title('Forced vibration response')
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/solving-the-coupled-eom-by-direct-integration/img2.jpg"
	alt="Forced vibration response of each mass/floor level. | EngineeringSkills.com"
	width="1000"
	height="521"
	caption="Fig 2. Forced vibration response of each mass/floor level."
/>

Everything is as we might expect. We observe the influence of a transient component of response that decays away due to the arbitrary damping specified above. This leaves us with a steady-state response at the forcing frequency; all things we should expect from a harmonic excitation, given what we know about SDoF systems already.

So, let's review; we have solved the equations of motion which is what we set out to do. But, the underlying dynamic characteristics of the system remain unclear. For example without simulating the response we can't tell what forcing frequencies our structure is particularly susceptible to.

We commented previously that our three-DoF structure has three natural frequencies. But what are they? Once we know these, we unlock some very valuable information about our structure and its response to dynamic loading. This leads us to the large topic of modal analysis we'll cover in the next section. Modal analysis will provide us with an alternative means of simulating the response of our structure, but along the way it also reveals information about the underlying dynamic characteristics of the system.

Before we wrap up this section, we'll briefly review some tools for solving differential equations. These will be useful to know and also provide us with a validation comparison for the response we've just calculated.


Please log in or enroll to continue

Please log in or enroll to continue