
One of the key assumptions in this project is that there is no interaction between pedestrians (i.e. no traffic effects) and that here is no human-structure interaction. As we've said already, this greatly simplifies the simulation task because now we can rely on superposition. So, to simulate the dynamic influence of a pedestrian crowd all we need to do is superimpose the responses generated by $N$ pedestrians walking individually across the bridge at different points in time.

For our simulation, we will generate normally distributed random walking velocities. We could also generate random body-weights, in fact we could randomise any of the pedestrian parameters. We'll stick with walking speed for now as this demonstrates the concept and has the biggest impact.

Our approach will be to loop through each pedestrian, generate their GRF and simulate the response of the structure. After all responses have been calculated, we'll superimpose them to determine the overall structural response during the simulation window.

Let's start by assuming that $N$ pedestrians cross the bridge in a 30 minute window. For simplicity we'll assume that the arrival times (onto the bridge) are uniformly distributed in this window. An alternative would be to assume a Poisson distribution of arrival times. This is something you can look into if developing a more rigorous simulation.

```python
N = 100 #Number of pedestrians that cross the bridge in the time window
window = 30*60 #(s) Simulation window
buffer = 200 #(s) Additional seconds to allow simulation of response beyond window length (late finishers)
mp = 80 #(kg) Pedestrian mass
G = 9.81*mp #(N) Static weight of pedestrian
```

Next we can use Numpy to generate $N$ random start times and walking velocities.

```python
#Random variables
tStart = np.random.uniform(low=0.0, high=window, size=N) #Uniformly distributed start times
Vp = np.random.normal(loc=1.3, scale=0.125, size=N) #Normally distributed walking velocities
```

Now we can loop through the crowd and calculate the modal force and response generated by each pedestrian. We'll save the force and response for each person as a separate row in a 2D matrix for forces and responses. Note that this is not the most efficient way of running this simulation as we need to call our Duhamel function $N$ times. It would be more efficient to calculate the combined crowd force and call the Duhamel function once. But coding the less efficient way preserves the individual pedestrian responses. This is a matter of personal preference. If your simulation is taking too long to run, just reduce the duration or number of pedestrians (helpful while you're developing the code).

```python
#Set up the simulation time vector
tMax = window + buffer #(s) Max time
time = np.arange(0, tMax+delT, delT)

#Initialise containers to hold the individual forces and responses calculated for each pedestrian
crowdForce = np.zeros([N,len(time)])
crowdResponse = np.zeros([N,len(time)])

#For each pedestrian...
for i, n in enumerate(np.arange(N)):
    vp = Vp[i] #(m/s) Walking velocity
    startTime = tStart[i] #(s) Start time
    tCross = L/vp #(s) Crossing time
    tEnd = startTime + tCross #(s) Finish time

    fv = 0.35*vp**3 - 1.59*vp**2 + 2.93*vp #(Hz) Pacing frequency
    DLF = 0.41*(fv-0.95) #Dynamic load factor

    timeVector = np.arange(0, tCross+delT, delT) #Time vector for this pedestrian
    Fv = G + abs(G*DLF*np.sin(2*math.pi*(fv/2)*timeVector)) #Static + Dynamic GRF (ignore static component)

    xp = vp*timeVector #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

    response = Duhamel(timeVector, Fn) #Response calculated using the Duhamel integral function

    #Save the GRF and response for this pedestrian at the correct position in the overal simulation records
    iStart = round(startTime/delT) #Index for start time
    crowdForce[i, iStart:iStart+len(Fn)] = Fn
    crowdResponse[i,iStart:iStart+len(Fn)] = response
```

Now we have the modal force and response for each individual pedestrian shifted along the time axis to match their random start time. Let's plot these to get a picture of what's happening.

```python
fig, axes = plt.subplots(figsize=(14,10),nrows=2,ncols=1)

for i in np.arange(len(crowdForce)):
    axes[0].plot(time,crowdForce[i,:],'-')
    axes[1].plot(time,-crowdResponse[i,:],'-')

#Housekeeping
axes[0].plot([window, window],[0,np.max(crowdForce)],'r--')
axes[0].plot([window+buffer, window+buffer],[0,np.max(crowdForce)],'r--')
axes[0].set_xlabel('time (s)')
axes[0].set_ylabel('Force (N)')
axes[0].set_title('Individual Modal forces')
axes[0].set_xlim([0,tMax])
# axes[0].set_xlim([startTime,startTime+tCross])
axes[0].grid()

axes[1].plot([window, window],[0,-np.max(crowdResponse)],'r--')
axes[1].plot([window+buffer, window+buffer],[0,-np.max(crowdResponse)],'r--')
axes[1].set_xlabel('time (s)')
axes[1].set_ylabel('Disp (m)')
axes[1].set_title('Individual Modal responses')
axes[1].set_xlim([0,tMax])
# axes[1].set_xlim([startTime,startTime+tCross])
axes[1].grid()

plt.show()
```

<Image
	src="/images/courses/simulating-crowd-induced-vibrations-using-the-duhamel-integral/dynamic-analysis-bridge-n-pedestrian-crowd/img1.jpg"
	alt="Modal forces (top) and modal responses (bottom), generated by N pedestrians. | EngineeringSkills.com"
	width="1000"
	height="689"
	caption="Fig 1. Modal forces (top) and modal responses (bottom), generated by N pedestrians."
/>

Now we can calculate the combined crowd-induced loading and response simply by superimposing the data for each pedestrian individually.

```python
#Sum across rows of crowdForce and crowdResponse
F_Crowd = sum(crowdForce)
Res_crowd = sum(crowdResponse)
```

We can also make use of our Peaks function to extract the response envelope.

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

Finally, we can plot the crowd loading and total modal response. Remember, this is the modal response so it corresponds to the response at the mid-span location of the beam. The response at any other point along the beam can be obtained directly from this by scaling using the mode shape $\phi$.

```python
fig, axes = plt.subplots(figsize=(14,10),nrows=2,ncols=1)

axes[0].plot(time,F_Crowd,'-')
axes[1].plot(time,-Res_crowd,'-')
axes[1].plot(times,-peaks,'r-')

axes[0].plot([window, window],[0,max(F_Crowd)],'r--')
axes[0].plot([window+buffer, window+buffer],[0,max(F_Crowd)],'r--')
axes[0].set_xlabel('time (s)')
axes[0].set_ylabel('Disp (m)')
axes[0].set_title('Modal responses (varying bridge mass)')
axes[0].set_xlim([0,tMax])
axes[0].grid()

axes[1].set_xlabel('time (s)')
axes[1].set_ylabel('Disp (m)')
axes[1].set_title('Modal responses (varying bridge mass)')
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-n-pedestrian-crowd/img2.jpg"
	alt="Combined modal force (top) and modal response (bottom), generated by N pedestrians. | EngineeringSkills.com"
	width="1000"
	height="694"
	caption="Fig 2. Combined modal force (top) and modal response (bottom), generated by N pedestrians."
/>

Now that we have the crowd-induced dynamic displacement, we could numerically differentiate this vector of numbers twice to obtain the vertical acceleration of the bridge/beam. This would be a good next step as vibration serviceability limits are usually reported in terms of acceleration. We'll skip this step and move on to building an animation of our simulation.


Please log in to continue