Dynamic Reliability Analysis of Large-Span Structures under Crowd Bouncing Excitation

Abstract

Bouncing is one of the most common human crowd activities on civil infrastructures such as sports stadiums and concert halls, where the audience tends to make their bodies jump up and down to celebrate or participate in sport and musical events. Dynamic loads are thus generated and exerted on the structures, giving unpleasant structural vibration, which may affect the functionality of the structure or even lead to a panic of the crowd. Although researchers have studied human-induced vibration from many perspectives including load models, calculation methods, criteria for serviceability evaluation, etc., there has been minimal work regarding crowd-induced reliability analysis, mainly because the stochastic feature of the crowd load as well as the mechanism describing the crowd–structure interaction is still not clear. In this paper, a framework to calculate crowd-induced structural vibration that considers the crowd–structure interaction effect is proposed and is validated through an experimental test. The dynamic parameters of the bouncing person in the crowd are adopted from a previous statistical study. The feasibility of a probability density evolution method (PDEM) is proved to be effective to calculate structural stochastic vibration under the bouncing crowd. The dynamic reliability of the structure is thus analyzed based on the stochastic responses. Results show that the consideration of the crowd–structure interaction effect significantly affects the dynamic reliability, which is also dependent on various factors including bouncing frequency, failure criteria, limit threshold, human model parameter distribution, etc. This paper provides a foundation for the performance-based vibration serviceability design of large-span structures.

1. Introduction

In recent decades, large-span structures become more and more popular in structural design due to the requirement from building functionalities as well as from an aesthetic point of view, especially for public facilities including sports stadiums, transport stations, etc. While the developed construction technologies and new materials are capable to ensure safety and durability, large-span structures such as sports stadiums usually suffer from unpleasant vibrations when human crowd activities such as crowd bouncing take place, leading to serviceability problems of the structure [1]. The well-known incident of the Millennium Bridge triggered research on the human-induced structural serviceability problem [2]. Researchers have studied this problem from various perspectives, including the human-induced load models [3,4], structural calculation methods [5,6], comfort criteria [7], and vibration control technologies [8]. The human-induced vibration serviceability problem has nowadays become an important and sometimes dominant issue that must be considered in the design stage of the large-span structures.

Although the research regarding the structural vibration serviceability problem has gained tremendous popularity among many researchers in the past several decades, the dynamic reliability problem of such vibration has rarely been investigated. In most design codes around the world, the serviceability of large-span structures is satisfied according to the allowable response method [9,10,11,12,13,14], indicating that when the human-induced structural response exceeds a predetermined value, the structure is determined as a failure from a serviceability point of view. However, it is already widely acknowledged that human-induced excitation is highly stochastic, featured by the so-called inter- and intra-subject variability [15]. To fully reflect the stochastic effect, the concept of dynamic reliability needs to be adopted. The failure of the structure should not be based on whether the structural response exceeds the limit value or not. Instead, the dynamic reliability of the structure under the crowd excitation needs to be higher than a predetermined target value.

It is known that the crowd–structure interaction (CSI) has a significant effect on structural vibration, especially when the mass of the crowd is not negligible compared with the mass of the structure, which, unfortunately, is always the case for the light-weight large-span structures [16,17,18,19]. To develop the dynamic reliability analysis method for the human-induced structural vibration, an analytical method that quantitatively considers the interaction effect is necessary. Some researchers use a single degree of freedom (DOF) system to represent the crowd and the structure and thus establish a 2 DOF equation of motion of the crowd–structure interaction system [20], which is unable to consider higher vibration modes that may interact with higher harmonics of the crowd load. Nimmen et al. proposed a detailed crowd model as well as a simplified method that makes it possible to evaluate the interaction effect for walking excitation [21,22]. However, there is still no widely accepted solution to this problem.

One more factor that may affect the CSI result is that the dynamic parameters of the human bodies in the crowd are not easily determined. Researchers have found that the human body could be represented by a spring-mass damper (SMD) model with its own mass, stiffness, and damping [5,23]. The values of these dynamic parameters are usually investigated through modal analysis on the empty and the crowd-occupied structure [16,24]. However, when the structure is occupied by an active crowd, the operational modal analysis becomes doubtful because the structure is under forced vibration. Moreover, the dynamic parameters of the human body are highly stochastic, and a large number of tests are necessary to obtain their probability distribution as well as their dependency on the motion frequency. In 2017 and 2019, the particle filter technique is adopted to estimate the dynamic properties of walking [25] and bouncing people [26], giving a reasonable description of the SMD model parameters, which can be used for further analysis of the CSI effect on the structural response calculation.

Due to the huge computational cost for reliability analysis, an efficient technique to calculate stochastic vibration is usually necessary, especially when a refined structural finite-element model is adopted. For such cases, the widely adopted Monte Carlo simulation (MCS), although versatile for stochastic vibration calculation, has its drawback of high computational cost in practical applications [27]. On the other hand, the probability density evolution method (PDEM) is popular for its capability of reducing the computational cost by the selection of limited representative points in the calculation process [28] and has the potential to be used in the dynamic reliability analysis of crowd-induced vibration.

From the above statement, three difficulties that hinder the development of the dynamic reliability for crowd-induced serviceability problems can be roughly seen, including (a) practical mechanism reflecting CSI effect, (b) stochastic feature of human body parameters, and (c) efficient stochastic vibration calculation technique. Therefore, the aim of this research is to develop an efficient framework for dynamic reliability analysis of structures suffering crowd bouncing load from the perspective of vibration serviceability and provide analytical results for the reference of engineering practice. In this paper, the governing dynamic equations with and without CSI are first derived and proposed. Each individual is regarded as a SMD model, and thus the effect of physical properties of each bouncing person is included in the analysis. An experiment quantitatively showing the CSI effect by including physical properties of the human body is carried out for validation and the effect of CSI is further investigated by numerical simulation. PDEM is then used in this study to calculate the stochastic response under crowd excitation efficiently, and the dynamic reliability is thus analyzed. The randomness originating from each bouncing person is included using the probability distribution extracted from a previous study [26]. Finally, the dependency of the dynamic reliability on various factors including the CSI effect, reliability criteria, bouncing frequency, limit threshold, and parameter distribution is investigated for engineering practice.

2. Analytical Model of Crowd–Structure Interaction

2.1. Human–Structure Coupled System

In this section, an analytical model of the crowd–structure interaction is given. The description of this analytical model starts from the SMD model representing single person bouncing on an arbitrary structure, which is illustrated by Figure 1. The equation of motion for single-person bouncing is expressed as Equation (1):

$$m_0{\ddot{u}}_0+c_0\left({\dot{u}}_0-{\dot{v}}_0\right)+k_0\left(u_0-v_0\right)=p_0\left(t\right)$$

(1)

where m₀, c₀, and k₀ are the mass, damping, and stiffness of the human body, u₀ is the displacement of the human system, v₀ is the structural displacement response of the structure at the human excitation point, and p₀(t) is the time history of the biomechanical, which could be defined by Fourier series as:

$$p_0\left(t\right)=a_0+{\displaystyle \sum _i\left(a_i\cos \left(2\pi f_{b}t\right)+b_i\sin \left(2\pi f_{b}t\right)\right)}$$

(2)

in which a₀, a_n, b_n are Fourier coefficients, f_b is the motion frequency, and n is the order of the biomechanical force. The nth biomechanical load factor (BLF) is defined by:

$${\mathrm{BLF}}_{0n}=\frac{\sqrt{a_n^2+b_n^2}}{m_{0}g}$$

(3)

From the force equilibrium at the surface of the structure, the contact force F₀(t) between the human body and the structure is equal to the human’s inertia force, as shown by Equation (4):

$$F_0\left(t\right)=-m_0{\ddot{u}}_0$$

(4)

The equation of motion of the structure is written for each mode according to the modal decomposition principle. For the jth mode of the structure, its equation of motion is:

$$M_{\mathrm{s},}{}_j{\ddot{q}}_j+C_{\mathrm{s},}{}_j{\dot{q}}_j+K_{\mathrm{s},}{}_j{q}_j=F_0\left(t\right){\varphi }_{j0}$$

(5)

where M_s,j, C_s,j, and K_s,j are the mass, damping, and stiffness of the jth decomposed mode, q_j is the jth modal coordinate, and ϕ_j0 is the jth mode shape value at the point of the bouncing excitation.

Through modal superposition, the physical displacement response at the excitation point is calculated as:

$$v_0={\displaystyle \sum _j{q}_j{\varphi }_{j0}}$$

(6)

It is clearly observed that Equations (1), (4) and (5) are coupled together, showing that the human mechanical system parameters have an influence on the structural response, which in return affects the human body response as well as the contact force.

2.2. Crowd–Structure Coupled System

The above derivation is extended to the case where multiple bouncing people excite the structure, which happens more frequently in engineering practice (e.g., audience celebrating in a grandstand of a sports stadium) and is more likely to give vibration serviceability problems compared to the single-person case. The crowd bouncing case is illustrated in Figure 2.

For the crowd bouncing case, each person in the crowd is represented by a SMD model with a pair of biomechanical forces. For the ith person in the bouncing crowd, the equation of motion is indicated by:

$$m_i{\ddot{u}}_i+c_i\left({\dot{u}}_i-{\dot{v}}_i\right)+k_i\left(u_i-v_i\right)=p_i\left(t\right)$$

(7)

in which the symbols share the same definition as in Equation (1).

Considering the multiple excitation points from the crowd, the equation of motion of the jth mode of the structure becomes:

$$M_{\mathrm{s},}{}_j{\ddot{q}}_j+C_{\mathrm{s},}{}_j{\dot{q}}_j+K_{\mathrm{s},}{}_j{q}_j=-{\displaystyle \sum _i{m}_i{\ddot{u}}_i{\varphi }_{ji}}$$

(8)

where ϕ_ji indicates the jth mode shape value at the location of the ith bouncing person.

From the modal superposition principle, the structural displacement at the location of the ith bouncing person is expressed:

$$v_i={\displaystyle \sum _j{q}_j{\varphi }_{ji}}$$

(9)

The combination of Equations (7)–(9) governs the crowd–structure coupling system. It is shown that the structural response is affected by all bouncing people in the crowd. In return, the human model parameters of each bouncing person in the crowd have an influence on the structural response.

2.3. Governing Dynamic Equations with and without Interaction Effect

The governing equations expressed by Equations (7)–(9) need to be solved to obtain the structural responses. For the case of N people bouncing on a structure decomposed by Y vibration modes, the total number of degree-of-freedom is N + Y. The equations rewritten in the matrix form is expressed by:

$$\left[\begin{array}{cc}m& O\\ Z_{\mathrm{m}}& M_{\mathrm{s}}\end{array}\right]\left[\begin{array}{c}\ddot{U}\\ \ddot{q}\end{array}\right]+\left[\begin{array}{cc}c& Z_{\mathrm{c}}\\ O& C_{\mathrm{s}}\end{array}\right]\left[\begin{array}{c}\dot{U}\\ \dot{q}\end{array}\right]+\left[\begin{array}{cc}k& Z_{\mathrm{k}}\\ O& K_{\mathrm{s}}\end{array}\right]\left[\begin{array}{c}U\\ q\end{array}\right]=\left[\begin{array}{c}P\\ O\end{array}\right]$$

(10)

in which

$$m={\left[\begin{array}{ccc}{m}_1& 0& 0\\ 0& \ddots & 0\\ 0& 0& m_N\end{array}\right]}_{N\times N},c={\left[\begin{array}{ccc}{c}_1& 0& 0\\ 0& \ddots & 0\\ 0& 0& c_N\end{array}\right]}_{N\times N},k={\left[\begin{array}{ccc}{k}_1& 0& 0\\ 0& \ddots & 0\\ 0& 0& k_N\end{array}\right]}_{N\times N}$$

(11)

$$M_{\mathrm{s}}={\left[\begin{array}{ccc}{M}_{\mathrm{s},1}& 0& 0\\ 0& \ddots & 0\\ 0& 0& M_{\mathrm{s},Y}\end{array}\right]}_{Y\times Y},C_{\mathrm{s}}={\left[\begin{array}{ccc}{C}_{\mathrm{s},1}& 0& 0\\ 0& \ddots & 0\\ 0& 0& C_{\mathrm{s},Y}\end{array}\right]}_{Y\times Y},K_{\mathrm{s}}={\left[\begin{array}{ccc}{K}_{\mathrm{s},1}& 0& 0\\ 0& \ddots & 0\\ 0& 0& K_{\mathrm{s},Y}\end{array}\right]}_{Y\times Y}$$

(12)

$$Z_{\mathrm{m}}=\left[\begin{array}{ccc}{m}_1{\varphi }_{11}& \cdots & m_N{\varphi }_{1N}\\ ⋮& \ddots & ⋮\\ m_1{\varphi }_{Y1}& \cdots & m_N{\varphi }_{YN}\end{array}\right],Z_{\mathrm{c}}=-\left[\begin{array}{ccc}{c}_1{\varphi }_{11}& \cdots & c_1{\varphi }_{Y1}\\ ⋮& \ddots & ⋮\\ c_N{\varphi }_{1N}& \cdots & c_N{\varphi }_{YN}\end{array}\right],Z_{\mathrm{k}}=-\left[\begin{array}{ccc}{k}_1{\varphi }_{11}& \cdots & c_k{\varphi }_{Y1}\\ ⋮& \ddots & ⋮\\ k{\varphi }_{1N}& \cdots & k_N{\varphi }_{YN}\end{array}\right]$$

(13)

$$U={\left[\begin{array}{ccc}{u}_1& \cdots & u_N\end{array}\right]}^{\mathrm{T}},P={\left[\begin{array}{ccc}{p}_1\left(t\right)& \cdots & p_N\left(t\right)\end{array}\right]}^{\mathrm{T}},q={\left[\begin{array}{ccc}{q}_1\left(t\right)& \cdots & q_Y\left(t\right)\end{array}\right]}^{\mathrm{T}}$$

(14)

and O indicate a zero matrix. In this work, a three-order Fourier series model is adopted to simulate p_i(t).

From Equations (10)–(14), it is observed that the crowd and the structure are coupled together through the matrices Z_c and Z_k. For the purpose of comparing the structural response with and without the CSI effect, the equation of motion without the coupling term is obtained by replacing them with zero matrices, as shown by Equation (15):

$$\left[\begin{array}{cc}m& O\\ Z_{\mathrm{m}}& M_{\mathrm{s}}\end{array}\right]\left[\begin{array}{c}\ddot{U}\\ \ddot{q}\end{array}\right]+\left[\begin{array}{cc}c& O\\ O& C_{\mathrm{s}}\end{array}\right]\left[\begin{array}{c}\dot{U}\\ \dot{q}\end{array}\right]+\left[\begin{array}{cc}k& O\\ O& K_{\mathrm{s}}\end{array}\right]\left[\begin{array}{c}U\\ q\end{array}\right]=\left[\begin{array}{c}P\\ O\end{array}\right]$$

(15)

Equations (10) and (15) can be easily solved by the typical numerical algorithm such as the Runge–Kutta algorithm or the Newmark algorithm once the parameters are determined. It is noteworthy that the above analysis applies to all types of structures since only the modal parameters are needed.

3. Numerical and Experimental Test for Structural Response Calculation

In this section, a numerical example is given to calculate the structural acceleration response following the procedure proposed above. The crowd parameter settings and the structural analytical model are briefly introduced and the structural responses with and without the consideration of the CSI effect are compared.

3.1. Analytical Model of a Large-Span Structure

A large-span structure with a size of 10 m × 6 m is adopted as the structural model for the numerical example in this section. The floor is with line supports two long and two short sides. The dynamic properties are acquired through modal tests conducted in advance and are listed in

Table 1. Dynamic properties of the structural model.

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6
Modal mass (kg)	8583	2587	9625	2423	2898	4900
Frequency (Hz)	3.500	6.150	6.750	14.120	15.190	18.100
Damping ratio (%)	0.374	0.514	0.614	0.913	0.666	1.497

. The mode shapes of the structure are depicted in Figure 3. Considering crowd bouncing usually results in linear vibration, the modal decomposition method is adopted using the modal parameters in Table 1 for efficiency instead of a finite element model.

3.2. Numerical Example Showing CSI Effect

As stated in the introduction section, the randomness of the structural responses originates from the stochastic physical parameters of the crowd. In a previous study, the physical parameters of a bouncing person as well as their probability distribution have been investigated. It is reported that the physical parameters of a bouncing person, including natural frequency, damping ratio, and BLFs, follow a skew-normal distribution [26] defined by its location parameter μ, scale parameter σ, and shape parameter α, as shown by Equation (16):

$$h\left(x\right)=\frac{2}{\sigma \sqrt{2\pi }}{e}^{-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}}{\displaystyle {\int }_{-\infty }^{\alpha \left(\frac{x-\mu }{\sigma }\right)}\frac{1}{\sqrt{2\pi }}{e}^{-\frac{t^2}{2}}\mathrm{d}t}$$

(16)

where h(x) describes the probability density function (PDF) of parameter x that represents either of the physical parameters of a bouncing person.

Moreover, the mass of crowd follows a normal distribution defined by parameter μ and σ according to Ref [29]. The coefficients of the above distribution when the bouncing frequency f_b equals 1.75 Hz are listed in

Table 2. Distribution coefficients of parameters for the analytical model, adapted from Refs. [26,29].

	Natural Frequency (Hz)	Damping Ratio (%)	BLF1	BLF2	BLF3	Mass (kg)
μ	1.73	0.13	0.20	0.15	0.01	62.80
σ	0.37	0.16	0.10	0.13	0.06	10.90
α	−3.52	4.00	2.61	0.75	8.01	——

.

In this paper, random variables representing the physical parameters of each bouncing person in the crowd are generated following the predetermined probability distribution. In this manner, the crowd is coupled with the structure in the way described in Section 2. Structural responses under the bouncing crowd with six different crowd sizes were calculated. Considering the size of the structure, the number of people in each crowd is decided to be 1, 3, 9, 15, 21, and 27, and their geometrical distribution of this crowd (i.e., the excitation points of the crowd) is illustrated in Figure 4. Because the structure has a fundamental frequency of around 3.500 Hz (See Table 1), a bouncing frequency of 1.750 Hz, which can be easily achieved for the audience activities in a sports stadium or in a concert hall, is analyzed in this numerical example.

As a typical illustration of the CSI effect, the structural mid-span acceleration with nine bouncing people is calculated according to Equations (10) and (15), for structural response with and without considering CSI, respectively. The calculated responses are given in Figure 5. Note that the human model parameters of the bouncing crowd are generated following the distribution given in Table 2. It is shown that the response with CSI is smaller than that without CSI, which is explained by the fact that the people in the bouncing crowd act as additional dampers that absorb energy from the structural vibration. A CSI index ε is defined following Equation (17) to show the error if the CSI effect is considered.

$$\epsilon =\frac{{\mathrm{RMS}}_{\mathrm{withoutCSI}}-{\mathrm{RMS}}_{\mathrm{withCSI}}}{{\mathrm{RMS}}_{\mathrm{withCSI}}}$$

(17)

In this calculation, the root mean square (RMS) value of the time history without CSI is calculated to be 4.78 m/s², while the response considering CSI is 2.95 m/s². Following the above definition, the CSI index is calculated to be 62.03%, clearly showing that large errors may occur if the CSI effect is not properly considered in the structural response prediction at the design stage.

The CSI index for different numbers of bouncing people is calculated for each crowd size. Because the physical parameters of the people in the bouncing crowd are random variables (See Table 2 for their distribution), the structural responses with and without the CSI effect are also random. The above process is repeated three times to show the increasing trend of the CSI index against the increase in number of people, as shown by Figure 6. With more bouncing people, the CSI effect tends to become larger, mainly because the mass ratio between the crowd and the structure becomes higher. The results show that the error from the negligence of CSI can reach up to 150% for such a large-span structure.

3.3. Formatting of Mathematical Components

To validate the procedure of the structural response calculation under crowd bouncing excitation, an experimental test was conducted on a prestressed concrete plate with a weight of 16.5 t. The dynamic properties of this plate have been given in Section 3.1. Accelerometers were attached at the bottom of the plate to capture the structural acceleration. A test participant with a weight of 59.7 kg was asked to bounce at the mid-point of the plate with a bouncing frequency of 1.75 Hz under the guidance of a metronome. The physical parameters of this participant had been identified in advance using an inverse analysis technique described in Ref. [26] and are listed here in

Table 3. Physical parameters of the bouncing person.

Parameter	Natural Frequency (Hz)	Damping Ratio (%)	BLF1	BLF2	BLF3
Value	1.49	41	0.345	0.267	0.046

. The overview of the experimental setup is shown in Figure 7.

The test participant was asked to bounce on the structure three times, each of which lasted around 30 s. The structural mid-span acceleration of Test I was plotted in Figure 8 together with the prediction given by Equations (10) and (15). It was expected that the measured structural acceleration was lower than the prediction because in the real case the test participant could not keep his bouncing frequency as a constant, while in the calculation procedure the bouncing frequency was fixed at 1.75 Hz to assure resonance. However, it is clear that the prediction with the consideration of CSI effect, i.e., responses given by Equation (10), is much closer and sometimes equals to the measured acceleration than the one without CSI effect.

The maximum value of the mid-span acceleration of each bouncing test was listed in

Table 4. Comparison of predicted and measured mid-span RMS acceleration.

Case	No CSI	With CSI	Test 1	Test 2	Test 3
Value (m/s2)	1.48	1.14	0.93	0.80	0.87

and compared with the prediction with and without the CSI effect. It is shown that although both predicted values are higher than the measured ones, the prediction with the CSI effect is much closer to the measured values, indicating that the proposed method for structural response calculation is reasonable.

4. PDEM-Based Stochastic Vibration Analysis and Its Verification

4.1. Formatting of Mathematical Components

Because of the variability of the human-induced load, the structural responses under crowd excitation should be considered from the perspective of random vibration. In view of this, the governing equation expressed by Equation (10) is rewritten in the form of Equation (18):

$$M\left(\Theta \right)\ddot{X}\left(\Theta ,t\right)+C\left(\Theta \right)\dot{X}\left(\Theta ,t\right)+K\left(\Theta \right)X\left(\Theta ,t\right)=F\left(\Theta ,t\right)$$

(18)

where

$$X={\left[\begin{array}{cc}U& q\end{array}\right]}^{\mathrm{T}},F={\left[\begin{array}{cc}P& O\end{array}\right]}^{\mathrm{T}}$$

(19)

and K, C, and M represent the matrices in Equation (10) before X and its derivatives.

The vector Θ in Equation (18) characterizes all random variables involved in the crowd–structure coupling system. For example, if the dynamic parameters of the structural model are considered deterministic, this vector will only contain human model parameters and has the form of Equation (20):

$$\Theta =[\begin{array}{c}\begin{array}{ccccccccc}{m}_1& \cdots & m_N& c_1& \cdots & c_N& k_1& \cdots & k_N\end{array}\end{array}{\begin{array}{ccccccccc}{ \mathrm{BLF}}_1{}_1& \cdots & {\mathrm{BLF}}_1{}_N& {\mathrm{BLF}}_2{}_1& \cdots & {\mathrm{BLF}}_2{}_N& {\mathrm{BLF}}_3{}_1& \cdots & {\mathrm{BLF}}_3{}_N\end{array}]}^{\mathrm{T}}$$

(20)

where BLF_ij denotes the ith order coefficient of the bouncing load for the jth person, and the others share the same definition with Equation (11). Therefore, the solution of Equation (18) is the function of Θ and t, and is expressed as:

$$X=H\left(\Theta ,t\right)$$

(21)

It is noteworthy that any physical quantities, i.e., acceleration, bending moment, etc., of the structural system could be expressed as a function of the solution X and/or its derivatives, as shown by Equation (22):

$$Z=F_Z\left(X\right)=H_Z\left(\Theta ,t\right)$$

(22)

where Z is the physical quantity of interest of the system and H_z is a deterministic vector operator that describes the physical mechanism of the system expressed by Equation (18).

According to the density evolution theory [28], the evolutionary joint PDF of (Z^T, Θ^T) ^T, denoted as p_ZΘ(z, θ, t), is governed by Equation (23).

$$\frac{\partial p_{Z\Theta }\left(z,\mathsf{\theta },t\right)}{\partial t}+{\displaystyle \sum _j{\dot{Z}}_j\left(\mathsf{\theta },t\right)}\frac{\partial p_{Z\Theta }\left(z,\mathsf{\theta },t\right)}{\partial z_j}=0$$

(23)

If only one physical quantity Z is of interest, Equation (23) will reduce to a one-dimensional form with the initial condition as shown by Equations (24) and (25), which is used more often in most situations.

$$\frac{\partial p_{Z\Theta }\left(z,\mathsf{\theta },t\right)}{\partial t}+\dot{Z}\left(\mathsf{\theta },t\right)\frac{\partial p_{Z\Theta }\left(z,\mathsf{\theta },t\right)}{\partial z}=0$$

(24)

$$p_{Z\Theta }\left(z,\mathsf{\theta },t_0\right)=\delta \left(z-z_0\right)p_{\Theta }\left(\mathsf{\theta }\right)$$

(25)

Equation (24) governs the PDF evolution of the PDF of the structural response which evolves with time and is thus denoted as the generalized density evolution equation (GDEE). The purpose to solve the stochastic vibration problem is to obtain the PDF of the response from that of the system parameters. For the crowd-induced structural vibration, the acceleration response is usually much related to the serviceability of the structure in the view of human comfort. Therefore, the physical quantity Z in this study is represented by the structural acceleration ${\ddot{v}}_i\left(t\right)$.

4.2. Procedure to Numerically Solve GDEE

For most engineering problems, the governing equation shown by Equation (24) needs to be solved numerically. In this process, the coefficient of the partial differential equation, $\dot{Z}\left(\mathsf{\theta },t\right)$, which can be calculated from Equation (22), should be obtained in the first place. This equation is solved by a point evolution method whose details are found in [27]. The steps using this method are briefly reviewed.

A number of n_pt representative points in the distribution domain Ω_Θ for vector Θ are firstly selected and are denoted as:

$${\mathsf{\theta }}_q=\left(\begin{array}{cccc}{\mathsf{\theta }}_{q,1}& {\mathsf{\theta }}_{q,2}& \cdots & {\mathsf{\theta }}_{q,N_{\mathsf{\theta }}}\end{array}\right),q=1,2,\cdots ,n_{\mathrm{pt}}$$

(26)

where N_θ denotes the dimension of vector Θ. The assigned probability corresponding to the qth representative point is:

$$P_q={\displaystyle {\mathsf{\int }}_{{\mathsf{\Omega }}_q}{p}_{\Theta }\left(\mathsf{\theta }\right)\mathrm{d}\mathsf{\theta }}$$

(27)

where Ω_q is the subdomain for Ω_Θ that satisfies Equation (28):

$$\{\begin{array}{c}{\mathsf{\Omega }}_q\cup {\mathsf{\Omega }}_p=\oslash ,\forall p\ne q\\ {\cup }_{q=1}^{n_{\mathrm{pt}}}{\mathsf{\Omega }}_q={\mathsf{\Omega }}_{\Theta }\end{array}$$

(28)

In each subdomain Ω_q, it is assumed that the coefficient $\dot{Z}\left(\mathsf{\theta },t\right)$ is invariant with respect to θ, as expressed as follows:

$$\dot{Z}\left(\mathsf{\theta },t\right)=\dot{Z}\left({\mathsf{\theta }}_q,t\right)$$

(29)

The governing equation is transformed to Equation (30) by integrating Equation (24) over Ω_q with respect to θ:

$$\frac{\partial p_q\left(z,t\right)}{\partial t}+\dot{Z}\left({\mathsf{\theta }}_q,t\right)\frac{\partial p_q\left(z,t\right)}{\partial z}=0,q=1,2,\cdots ,n_{\mathrm{pt}}$$

(30)

where

$$p_q\left(z,t\right)={\displaystyle {\int }_{{\mathsf{\Omega }}_q}{p}_{Z\Theta }\left(z,\mathsf{\theta },t\right)\mathrm{d}\mathsf{\theta }}$$

(31)

The corresponding initial condition is also transformed to:

$$p_q\left(z,t_0\right)=\delta \left(z-z_0\right)P_q,q=1,2,\cdots ,n_{\mathrm{pt}}$$

(32)

which could be solved through the finite difference method, giving the numerical solution of p_q(z, t). The PDF of Z at each time instant can be calculated by summation of the PDF of each representative point, as shown by Equation (33):

$$p_Z\left(z,t\right)={\displaystyle {\int }_{{\mathsf{\Omega }}_{\Theta }}{p}_{Z\Theta }\left(z,\mathsf{\theta },t\right)\mathrm{d}\mathsf{\theta }}={\displaystyle \sum _{q=1}^{n_{\mathrm{pt}}}{\displaystyle {\int }_{{\mathsf{\Omega }}_q}{p}_{Z\Theta }\left(z,\mathsf{\theta },t\right)\mathrm{d}\mathsf{\theta }}}={\displaystyle \sum _{q=1}^{n_{\mathrm{pt}}}{p}_q\left(z,t\right)}$$

(33)

The PDF expressed in the above equation serves as the basis to evaluate the structural dynamic reliability in further analysis.

4.3. Response Calculation through PDEM and Its Verification

Because the human physical parameters of the bouncing person in a crowd are random parameters, the structural acceleration under such a crowd excitation needs to be considered as a random process. In calculation, the GF-discrepancy method is used for point selection [30]. A difference scheme known as total variation diminishing (TVD) is adopted to solve the GDEE expressed by Equation (24), whose details are found in Ref [31], to obtain the PDF of the structural responses. The mean value and the standard deviation time histories are plotted in Figure 9 for the case described in Section 3.2, to characterize the randomness. As in Figure 9a, the mean time history shares a similar shape with its representative time history shown in Figure 5, while the standard deviation can reach as large as 1.2 m/s², indicating that the randomness of the physical parameters in the crowd could not be neglected.

Furthermore, the PDF evolution surface of the mid-span acceleration responses is obtained and illustrated in Figure 10. The extreme value of the PDF decreases over time within the first ten seconds, indicating that the structural responses are gradually stimulated, and the variation keeps increasing. After around 15 s, the extreme value becomes nearly constant, mainly because the structural responses have reached the stable stage. The PDF thus starts to evolve with a regulated pattern, as already illustrated in Figure 5 and Figure 8.

The above results are obtained from the PDEM method on the basis of 1000 deterministic analyses, i.e., 1000 points are selected using the GF-discrepancy method. To verify its accuracy, the distribution at some time instants is calculated through Monte-Carlo simulation (MCS) using 100,000 examples. The random distribution is compared in the form of cumulative density function (CDF) as shown by Figure 11. Good consistency is observed from the comparison of CDF from PDEM and MCS, and the results obtained by PDEM are thus demonstrated to be correct.

Based on the results of the random vibration analysis above, the dynamic reliability of the structure is evaluated, which will be discussed in the next section.

5. Dynamic Reliability Analysis for Crowd-Induced Structural Vibration

5.1. Failure Criteria

For the calculation of structural dynamic reliability, reasonable failure criteria need to be defined. Because the human comfort on the structure is mainly related to structural acceleration, most design codes require that the maximum acceleration response should not exceed a predefined limit to maintain the vibration serviceability of the structure. Once the acceleration exceeds such a limit, the structure is considered a failure from the perspective of vibration serviceability. If the structural acceleration at the mid-point of the structure is interested, the reliability could be expressed by the most frequently adopted first-passage criterion shown in Equation (34):

$$R\left(t\right)=P\left\{\left|{\ddot{v}}_{\mathrm{mid}}\left(\tau \right)\right|⩽\left[\ddot{v}\right],\tau \in \left[0,t\right]\right\}$$

(34)

in which R(t) is the time-dependent dynamic reliability, ${\ddot{v}}_{\mathrm{mid}}$ is the acceleration at the mid-point of the structure, and $\left[\ddot{v}\right]$ is the predefined limit threshold.

However, the mid-point of the structure may not always be the location where the largest structural acceleration occurs, especially when higher vibration modes are excited. Furthermore, the locations on which people are bouncing or standing are more interested because the structural acceleration at these locations is directly perceived by people. In this manner, it is more appropriate to define the dynamic reliability as the probability that the acceleration responses at a series of selected locations do not exceed the predefined limit threshold, as expressed by Equation (35):

$$R\left(t\right)=P\left\{{\displaystyle \underset{l=1}{\overset{L}{\cap }}\left|{\ddot{v}}_l\left(\tau \right)\right|⩽\left[\ddot{v}\right]},\tau \in \left[0,t\right]\right\}$$

(35)

where L is the number of locations of interest, and ${\ddot{v}}_l$ indicates the structural acceleration at the lth location.

5.2. Calculation of Dynamic Reliability

The reliability of a structure is usually defined as the probability that the structure finished its expected function within a certain time period. Li proposed a method to calculate the dynamic reliability by considering the structure as a probability dissipative system in 2020 [32], whose GDEE is in the form of:

$$\frac{\partial p_{U\Theta }\left(u,\mathsf{\theta },t\right)}{\partial t}+\dot{U}\left(\mathsf{\theta },t\right)\frac{\partial p_{U\Theta }\left(u,\mathsf{\theta },t\right)}{\partial u}=-\mathscr{H}\left[f\left(U\left(\mathsf{\theta },t\right)\right)\right]\cdot p_{U\Theta }\left(u,\mathsf{\theta },t\right)$$

(36)

where

$$\mathscr{H}\left[f\left(U\left(\mathsf{\theta },t\right)\right)\right]=\{\begin{array}{c}0,f\left(U\left(\mathsf{\theta },t\right)\right)\in {\mathsf{\Omega }}_{\mathrm{S}}\\ 1,f\left(U\left(\mathsf{\theta },t\right)\right)\in {\mathsf{\Omega }}_{\mathrm{D}}\end{array}$$

(37)

and u can be the time history of any physical quantity of the system, f (·) is a general function that links u to the structural responses of interest, and Ω_S and Ω_D indicate the safety domain and the failure domain of the system, respectively.

$$R\left(t\right)={\displaystyle {\int }_{-\infty }^{\infty }{p}_U\left(u,t\right)\mathrm{d}u}$$

(38)

5.3. Reliability Analysis Results

In this sub-section, the time-dependent dynamic reliability of the structure in Section 3.3 is analyzed following the procedure described in the previous section. Factors that may affect the reliability are analyzed, including CSI effect, reliability criteria, bouncing frequency, limit threshold, and parameter distribution. Note that the results in this section correspond to the same excitation described in Section 3.2, i.e., a crowd of nine people bouncing at a frequency f_b of 1.75 Hz, unless otherwise specified. Figure 12 exhibits a series of figures that show the relation between dynamic reliability and time. It is observed that all curves decrease with time, in accordance with the absorbing feature of the first-passage criteria.

Figure 12a shows the dependency of the dynamic reliability on the adopted failure criteria, in which the limit threshold equals 5.0 m/s² for both situations. If only the mid-point structural acceleration is interested, the reliability becomes stable around 0.21. However, if the structural accelerations at all bouncing locations are considered, the reliability decreases to 0.18. This phenomenon indicates that the mid-point may not always be the location with the largest acceleration response. It is better to consider the acceleration at more locations during vibration serviceability design.

The bouncing frequency is another important factor that highly affects the dynamic reliability of the structure. If the bouncing frequency or its harmonics are close to the fundamental frequencies of the structure, the structural responses tend to become much larger, and the dynamic reliability is thus decreased. In Figure 12b, it is clearly observed that the dynamic reliability significantly fluctuates with the bouncing frequency. When f_b equals 1.75 Hz, whose harmonics give rise to the resonant response of the structure, the dynamic reliability drops to its minimum value, while for other bouncing frequencies, the dynamic reliability becomes much higher.

It is also observed from Figure 12c that if the CSI effect is considered, the dynamic reliability is largely increased, especially for the case when f_b = 1.75 Hz, which is explained by the fact that the CSI mostly affects the resonant vibration. On the other hand, the serviceability problem usually occurs when the structure is under resonant excitation. Therefore, it is important to consider the CSI effect in the vibration serviceability design to avoid under-estimation of the dynamic reliability.

The influence of human parameter distribution is also checked in this study, as shown by Figure 12d. The line styles of the first column of the legend correspond to the results obtained by assuming the human body model parameters follow the distribution in Table 2, i.e., people in the bouncing crowd have independent parameters, while that of the second column assume that all people in the crowd share the same human model parameters, which equal to the mean value of their distributions, respectively. Results show that the ignorance of human model parameter distribution may lead to either an over-estimate or under-estimate of the structural reliability.

It is apparent that the dynamic reliability also changes with the predefined allowable limit. Figure 12e clearly shows the increasing trend of the dynamic reliability against the predefined threshold value.

6. Concluding Remarks

The dynamic reliability of a large-span structure is evaluated through PDEM based on the crowd–structure coupled model. The uncertainties of the coupling system come from the random distribution of the model parameters of the human body and the biomechanical forces in the bouncing crowd. The governing equations of motion of the coupled system with and without the CSI effect are given. The features of the random vibration of the coupled system are studied through PDEM and the results lead to the evaluation of dynamic reliability from the perspective of vibration serviceability. Two failure criteria in terms of one-point and multiple-point acceleration passage are adopted for the reliability evaluation, and the affecting factors on the reliability are discussed. The following conclusions are drawn from this study:

(1)

The proposed calculation procedure to consider the interaction effect in the crowd–structure coupled system can well predict structural responses and is validated to be reasonable through an experimental test.

(2)

The CSI highly affects the structural responses, especially when the crowd size is large.

(3)

Through comparison with traditional MCS, the PDEM is tested as capable to conduct human-induced random vibration analysis, which is the foundation of dynamic reliability calculation considering the CSI effect. The PDEM has a great potential when a refined model needs to be adopted or a more complex situation is considered, where MCS becomes unavailable due to its high computational cost.

(4)

The dynamic reliability of the large-span structure in terms of vibration serviceability is affected by many factors, including failure criteria, excitation frequency, limit threshold, distribution of human model parameters, and CSI effect.

The current design codes usually require that the structural vibration under human activities does not exceed a predefined limit to maintain the vibration serviceability of the structure. However, because of the uncertainties in the crowd–structure coupling system, this requirement is actually satisfied with a certain degree of probability rather than in a deterministic manner, thus leading to the concept of reliability. This paper presents the procedure of dynamic reliability analysis for large-span structures in terms of vibration serviceability. In future studies, the threshold of vibration perception of each person in the crowd could serve as another random variable in the system. This paper provides a foundation for achieving performance-based vibration serviceability design in the future.