
Since we're assuming that the bridge has only one modal frequency in the range of interest, we only need to simulate a single vibration mode. Conveniently for us and through the magic of modal analysis, this means that the dynamic analysis of the bridge reduces to simulating the response of a SDoF system. If we needed to simulate the combined response of multiple modes, this is easily done through modal superposition.

Our first task is to turn our GRF, $F_v$, into a modal force, $F_n$. This simply means multiplying it by the unity-normalised mode shape value at the pedestrian's location, $x_p$,

$$
F_n = F_v\phi (x_p)
\tag{1}
$$

Since our pedestrian is walking with a constant velocity, we can easily determine their position, $x_p$, as a function of time. From here we can determine the value of the mode shape as a function of time and directly multiply this by the GRF calculated earlier to obtain the modal load applied to the SDoF modal model as a function of time.

```python
xp = vp*time #Pedestrian position as a function of time
phi = np.sin(math.pi*xp/L) #Mode shape at pedestrian's location
Fn = Fv*phi #Modal force experienced by SDoF system
```

At this point it's helpful to visualise the modal force. Remember, this is the force applied to the first bridge/beam mode with a half-sine mode shape. We can think of this force as the influence of the GRF on mode one which is at a maximum at the centre-span position.

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1])
axes.plot(time,Fn, '-', label='Modal load')

#Housekeeping
axes.set_xlabel('time (s)')
axes.set_ylabel('Force (N)')
axes.set_title('Modal Force')
axes.set_xlim([0,tMax])
plt.grid()
plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/dynamic-analysis-bridge-1-pedestrian/img1.jpg"
	alt="Modal force. | EngineeringSkills.com"
	width="1000"
	height="378"
	caption="Fig 1. Modal force."
/>

We can clearly see that the dynamic influence of the pedestrian's GRF is increasing as they move towards the center of the bridge and diminishing thereafter. Now it's just a matter of providing this force along with the time vector directly into our Duhamel function. Since our function relies on globally defined variables, we need to be careful to redefine these for our bridge before calling the function.

```python
M = 2000 #(kg/m) Mass per unit length]
m = 0.5*M*L #(kg) Modal mass of mode 1
xi = 0.025 #Damping ratio
fn = 2.5 #(Hz) Bridge modal frequency
wn = 2*math.pi*fn #(rads/s) Angular modal frequency
wd = wn*math.sqrt(1-xi**2) #(rads/s) Damped angular modal frequency
```

```python
response = Duhamel(time, Fn) #Response calculated using the Duhamel integral function
```

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1])
axes.plot(time,-response, '-')

#Housekeeping
axes.set_xlabel('time (s)')
axes.set_ylabel('Disp (m)')
axes.set_title('Modal response (static+dynamic)')
axes.set_xlim([0,tMax])
plt.grid()
plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/dynamic-analysis-bridge-1-pedestrian/img2.jpg"
	alt=" Modal response (static and dynamic components). | EngineeringSkills.com"
	width="1000"
	height="382"
	caption="Fig 2. Modal response (static and dynamic components)."
/>

After plotting the modal response we can see the oscillation behaviour superimposed on top of the static deflection. We can visualise this more clearly if we separate out the static and dynamic components of the GRF and simulate their influence separately.

```python
Fn_static = G*phi #Static component of GRF
Fn_dynamic = abs(G*DLF*np.sin(2*math.pi*(fv/2)*time))*phi #Dynamic component of GRF

#Response due to constant magnitude moving load
response_static = Duhamel(time, Fn_static)
#Response due to footsteps (with static load component removed)
response_dynanmic = Duhamel(time, Fn_dynamic)
```

```python
fig, axes = plt.subplots(figsize=(14,10),nrows=2,ncols=1)
axes[0].plot(time,-response_dynanmic, '-', label='Dynamic')
axes[0].plot(time,-response_static, 'r-', label='Static')
axes[0].set_xlabel('time (s)')
axes[0].set_ylabel('Disp (m)')
axes[0].set_title('Modal response (separate components)')
axes[0].legend(loc='lower right')
axes[0].set_xlim([0,tMax])
axes[0].grid()

axes[1].plot(time,-response, '-', label='Combined')
axes[1].set_xlabel('time (s)')
axes[1].set_ylabel('Disp (m)')
axes[1].set_title('Modal response (combined)')
axes[1].legend(loc='lower right')
axes[1].set_xlim([0,tMax])
axes[1].grid()

plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/dynamic-analysis-bridge-1-pedestrian/img3.jpg"
	alt="Static and dynamic response separated (top), combined modal response (bottom). | EngineeringSkills.com"
	width="1000"
	height="728"
	caption="Fig 3. Static and dynamic response separated (top), combined modal response (bottom)."
/>

It will help us to more clearly visualise the response if we write a quick function to plot the displacement envelope. This will be particularly helpful later on when our response is derived from the influence of many pedestrians walking simultaneously. We'll do this by cycling through the displacement record and extracting the oscillation peaks. A peak will be identified when the slope of the displacement changes sign from positive to negative.

This won't be a perfect strategy as it will also identify 'minor' peaks due to frequencies other than the dominant oscillation frequency but it will be close enough for visualisation purposes. Again, we'll package the code into a function that can be easily called multiple times.

```python
def Peaks(disp, time):
  #Initialise containers to hold peaks and their times
  peaks = np.empty([1,0])
  times = np.empty([1,0])

  #Calculate slopes for each data point
  slopes = np.zeros(len(disp))
  for i, u in enumerate(disp):
      if(i<len(disp)-1):
          slopes[i] = disp[i+1]-disp[i]

  #Cycle through all slopes and pick out peaks
  for i, s in enumerate(slopes):
      if (i<len(slopes)-1):
          if(slopes[i+1]<0 and slopes[i]>0):
              peaks = np.append(peaks,disp[i])
              times = np.append(times,time[i])

  return [peaks, times]
```

```python
peaks, times = Peaks(response,time)
```

Now we can plot the response again but also plot an envelope line that spans between oscillation peaks.

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1])
axes.plot(time,-response, '-', label='Response')
axes.plot(times,-peaks,'r-', label='Response envelope')

#Housekeeping
axes.set_xlabel('time (s)')
axes.set_ylabel('Disp (m)')
axes.set_title('Modal response and envelope')
axes.legend(loc='lower right')
axes.set_xlim([0,tMax])
plt.grid()
plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/dynamic-analysis-bridge-1-pedestrian/img4.jpg"
	alt=" Modal response (static and dynamic components) with displacement envelope. | EngineeringSkills.com"
	width="1000"
	height="382"
	caption="Fig 4. Modal response (static and dynamic components) with displacement envelope."
/>

Now that we have a neater way of visualising the response, i.e. the response envelope, we can run a number of simulations to visualise the influence of any one of our simulation parameters. Let's consider the influence of bridge mass per unit length for example. We'll start by calculating the bridge stiffness based on our original mass and frequency parameters and then just loop through the different masses.

```python
k = m*wn**2 #(N/m) Original system stiffness
Masses = [1750, 2000, 2250] #(kg/m) Masses per unit length to test

fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1])

for M in Masses:
    m = 0.5*M*L #(kg) Modal mass of mode 1
    wn = math.sqrt(k/m)
    wd = wn*math.sqrt(1-xi**2) #(rads/s) Damped angular modal frequency

    response = Duhamel(time, Fn) #Response calculated using the Duhamel integral function
    peaks, times = Peaks(response,time)

    axes.plot(times,-peaks,'-', label=f'Response envelope M={M}')

#Housekeeping
axes.set_xlabel('time (s)')
axes.set_ylabel('Disp (m)')
axes.set_title('Modal responses (varying bridge mass)')
axes.legend(loc='lower right')
axes.set_xlim([0,tMax])
plt.grid()
plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/dynamic-analysis-bridge-1-pedestrian/img5.jpg"
	alt="Enveloped of response for three simulations with different bridge mass. | EngineeringSkills.com"
	width="1000"
	height="382"
	caption="Fig 5. Enveloped of response for three simulations with different bridge mass."
/>

As we might expect, as the mass increases the modal frequency reduces bringing it closer to the pacing frequency. This obviously leads to larger amplitude oscillations as the frequency ratio approaches one. We could repeat this process and test the sensitivity of our response to any of the input parameters - the process is more or less the same. Instead we'll do something a little more interesting and simulate the influence of multiple pedestrians crossing the bridge at the same time.


Please log in to continue