
As we said at the outset, there are two approaches to solving the Duhamel integral (or any integration really), analytical solution and numerical solution. By far the most practically useful method is numerical solution. Most of the interesting problems, from an engineering perspective, require numerical solution. Although a numerical solution will always be an approximation, any error introduced usually pales into insignificance when compared against other uncertainties bakes into our modelling.

Numerical solution comes with the huge advantages of simplicity and versatility. I've always been more of a programmer than mathematician so I may be biased. Having said all of this, we will demonstrate an analytical solution first, for completeness; even then we'll use Python to do the hard work.

## Analytical solution to ramp loading

Consider a load that increases linearly from $0$to a magnitude of $P_0$ between $t=0$ and $t=t_1$. We can represent the loading with a simple mathematical expression,

$$
p(t) = \frac{P_0 t}{t_1}
\tag{1}
$$

The fact that the applied loading can be captured with such a simple expression is our first clue that we might be able to solve this analytically. Assigning some arbitrary values to $P_0$ and $t$ we can generate a plot of the loading.

```python
# DEPENDENCIES and DEFAULTS
import math #.....................................Math functionality
import numpy as np #..............................Numpy for working with arrays
import matplotlib.pyplot as plt #.................Plotting functionality

P0 = 1000 #(N) Max load (arbitrary value for visualisation only)
t1 = 10 #(sec) Rise time (arbitrary value for visualisation only)
delT = 0.1 #(sec) Time-step
t = np.arange(0, t1+delT, delT) #Time vector
p = P0*t/t1 #Force vector

#Set up figure and axes
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1])
axes.plot(t,p)

#Housekeeping
axes.set_xlabel('time (sec)')
axes.set_ylabel('Force (N)')
axes.set_title('Applied ramp load')
axes.set_xlim([0,t1])
axes.set_ylim([0,P0])
plt.grid()
plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/analytical-solution/img1.jpg"
	alt="Ramp load. | EngineeringSkills.com"
	width="1000"
	height="375"
	caption="Fig 1. Ramp load."
/>

Next, we can insert our expression for the applied loading directly into the Duhamel integral expression.

$$
u(t) = \frac{1}{m\omega_d}\int_0^t \frac{P_0 \tau}{t_1} e^{-\xi\omega_n(t-\tau)}\sin \omega_d (t-\tau)\mathrm{d}\tau
\tag{2}
$$

In this case, we are only considering the response of the system during the application of the ramp load and not the free vibration that would follow. As such, damping will not have an appreciable influence on the response. We can therefore simplify the task by considering an undamped system. As a result, our expression reduces to,

$$
u(t) = \frac{P_0}{t_1 m\omega_n}\int_0^t \tau \sin \omega_n (t-\tau)\mathrm{d}\tau
\tag{3}
$$

Now we need to integrate this expression; for this we'll use SymPy.

```python
import sympy as sym #Import SymPy
sym.init_printing() #So we print nicely formatted equations

#Define some symbols
m = sym.Symbol('m')
w = sym.Symbol('w')
P0 = sym.Symbol('P0')
t1 = sym.Symbol('t1')
tau = sym.Symbol('tau')
t = sym.Symbol('t')

#Construct the function to integrate
f = tau * sym.sin(w*t-w*tau)

# Use SymPy to get the definite itegral with respect to tau between 0 and t
defInt = sym.integrate(f,(tau,0,t))

print('SymPy generated the following expression for the definite integral:')
sym.simplify(defInt) #Attempt to siplify the definite integral
```

SymPy generated the following expression for the definite integral:

$$
\begin{align*}
0 & \:\:\: \text{(for w=0)}\\
\frac{1}{w^2}(tw-\sin(tw)) & \:\:\: \text{(otherwise)}
\end{align*}
$$

Clearly we're not interested in the case where $w_n=0$, so rearranging the other expression and combining it with the constants from our response expression (equation 3) yields,

$$
u(t) = \frac{P_0}{t_1 m \omega_n}\left(\frac{t}{\omega_n}-\frac{\sin(\omega_n t)}{\omega_n^2}\right)
\tag{4}
$$

We can further rearrange this,

$$
u(t) = \frac{P_0}{k}\left[\frac{t}{t_1} - \frac{\sin(\omega_n t)}{\omega_n t_1}\right]
\tag{5}
$$

To complete the picture, we can plot the response. When plotting we will normalise the time axis by the rise time $t_1$ and normalise the displacement axis by the maximum static displacement, $P_0/k$.

```python
#Redefine constants following symbolic math operations
P0 = 1000 #(N) Max load
t1 = 10 #(sec) Rise time
delT = 0.1 #(sec) Time-step
t = np.arange(0, t1+delT, delT) #Time vector
```

We'll perform the response calculation for a range of different SDoF systems with the same mass but different stiffness, yielding a range of natural frequencies. For each system we'll specify the period of the system as a proportion of the rise time of the load.

```python
m=20 #Mass of the system
periodRange = [0.3,0.4,0.5] #Range of system periods as a propertion of rise time

#Initialise a figrue to plot onto
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1])

for pr in periodRange:
    T = pr*t1
    wn = 2*math.pi/T
    k=m*wn**2
    u = (P0/k)*((t/t1)-((np.sin(wn*t))/(wn*t1)))
    axes.plot(t/t1,u/(P0/k), label=f"T = {pr}t1")

#Housekeeping
axes.set_xlabel('t/t1')
axes.set_ylabel('Displacement ratio')
axes.set_title('SDoF system response to ramp loading')
axes.legend(loc='upper left')
axes.set_xlim([0,1])
plt.grid()
plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/analytical-solution/img2.jpg"
	alt="Response to ramp loading | EngineeringSkills.com"
	width="1000"
	height="378"
	caption="Fig 2. Response to ramp loading."
/>

This example is a relatively trivial case but it demonstrates the analytical integration route. As I mentioned, you're not likely to deploy analytical integration too often, so we'll leave this method behind and focus on numerical integration from here on.


Please log in to continue