
In this lecture we will plot the free vibration response of the system assuming the harmonic force is removed at $t=11$ seconds. To tackle this question we need to determine the constants of integration when the initial conditions are not zero. First we must determine the position and velocity of the system at the instant the harmonic force is removed at $t=11$ s. We can do this by evaluating the equations for $u(t=11)$ and $\dot{u}(t=11)$ derived on the previous lecture.

```python
t = 11 #Time instant

#Position
ut7 = math.e**(-xi*omega_n*t)*(A*np.sin(omega_d*t) + B*np.cos(omega_d*t)) #Transient response
us7 = O*((1-beta**2)*np.sin(omega*t) - 2*xi*beta*np.cos(omega*t)) #Steady-state response
u = ut7+us7 #Combined response
print(f'The position of the system at t=11 s is {round(u,4)} m')

#Velocity
term_1 = math.e**(-xi*omega_n*t)*omega_d*np.cos(omega_d*t)
- np.sin(omega_d*t)*xi*omega_n*math.e**(-xi*omega_n*t)

term_2 = -1*math.e**(-xi*omega_n*t)*omega_d*np.sin(omega_d*t)
- np.cos(omega_d*t)*xi*omega_n*math.e**(-xi*omega_n*t)

v = A*(term_1) + B*(term_2) + C1*(omega*np.cos(omega*t)) - C2*(omega*np.sin(omega*t))
print(f'The velocity of the system at t=11 s is {round(v,4)} m/s')
```

- `The position of the system at t=11 s is 0.0045 m`
- `The velocity of the system at t=11 s is -0.0191 m/s`

These values of $u(t)$ and $\dot{u}(t)$ now serve as the initial conditions when determining the free vibration response at $t=30$ s. Remember we are only concerned with the free vibration,

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

We must now determine the constants $A$ and $B$ with our new initial conditions. We can restate our expression for $\dot{u}(t)$ from the previous lecture.

$$
\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]
\tag{2}
$$

Evaluating both expressions for the initial conditions $u(t=11) =0.0045$ and $\dot{u}(t=11) = -0.0191$ results in two simultaneous equations,

```python
#Numerical evaluating the terms in square brackets...
a = math.e**(-xi*omega_n*t)*np.sin(omega_d*t) #Term in brackets beside A in equation for u(t)
b = math.e**(-xi*omega_n*t)*np.cos(omega_d*t) #Term in brackets beside B in equation for u(t)
print(f'a = {round(a,4)}')
print(f'b = {round(b,4)}')
print(f'term 1 = {round(term_1,3)}') #Term in brackets beside A in equation for v(t)
print(f'term 2 = {round(term_2,3)}') #Term in brackets beside B in equation for v(t)
```

- `a = 0.2167`
- `b = 0.289`
- `term 1 = 1.317`
- `term 2 = -1.029`

Therefore the system of equations, in martix form is,

$$
\begin{bmatrix}
    0.2167 & 0.289 \\
    1.317 & -1.029
\end{bmatrix}\times
\begin{bmatrix}
    A \\
    B
\end{bmatrix} =
\begin{bmatrix}
    0.0045 \\
    -0.0191
\end{bmatrix}
\tag{3}
$$

Solving for A and B...

```python
#Matrix
matrix = np.mat([[a,b],[term_1,term_2]])

#Vector
vector = np.array([[u],[v]])
constants = matrix.I*vector
A = constants[0].item()
B = constants[1].item()
print(f'A = {round(A,3)}')
print(f'B = {round(B,3)}')
```

- `A = -0.001`
- `B = 0.017`

Now we can evaluate the equation for free vibration at $t=30$ s.

```python
t = 30 #Time instant
a = math.e**(-xi*omega_n*t)*np.sin(omega_d*t) #Term in brackets beside A in equation for u(t)
b = math.e**(-xi*omega_n*t)*np.cos(omega_d*t) #Term in brackets beside B in equation for u(t)

u30 = (A*a + B*b).item()
print(f'The free vibration response at t=30s is {round(u30,4)} m')
```

- `The free vibration response at t=30 s is 0.0008 m`

As a final step, let's plot the free vibration response up to $t=30$ seconds. First define a time vector from $t=11$ seconds up to $t=30$ seconds.

```python
tmax = 30 #(sec) The max time
delta_t = 0.001 #(sec) The time-step
nPoints = tmax/delta_t #Number of equally spaced data points
t = np.linspace(11,tmax, int(nPoints)) # Time vector
```

And finally, calculate and plot the response with a horizontal red line indicating the magnitude of the response at $t=30$ seconds.

```python
ut = math.e**(-xi*omega_n*t)*(A*np.sin(omega_d*t) + B*np.cos(omega_d*t)) #Transient response

#Plotting (11-30s)
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,2])
axes.plot(t,ut,label='Transient response')
axes.plot([11, tmax],[u30, u30], 'r--', lw=0.75)
axes.set_xlim([11,30])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Displacement (m)')
axes.set_title('Displacement-time history (11-30 seconds)')
axes.grid()
axes.legend(loc='lower right')
plt.show()
```

<Image
	src="/images/courses/fundamentals-of-structural-dynamics-with-python/worked-example-3-part-c/img1.jpg"
	alt="Displacement-time history ($11-30$ seconds). | EngineeringSkills.com"
	width="1000"
	height="673"
/>
<Caption>Fig 1. Displacement-time history ($11-30$ seconds).</Caption>

That concludes this exercise. You should now have a good understanding of harmonic excitation and how to practically apply the equations derived in this section. In the next section we'll move away from the world of closed-form analytical solutions and start to investigate numerical techniques for determining dynamic response.


Please log in or enroll to continue