# Nonlinear Stochastic Analysis of Footbridge Lateral Vibration Based on Probability Density Evolution Method

Authors: Zheng Yang; Buyu Jia; Quansheng Yan; Xiaolin Yu; Yinghao Zhao
Journal: Shock and Vibration
 (2019)

---

## Abstract

Footbridge lateral vibration remains an unsolved problem and is characterized by the following: (1) pedestrians are sensitive to bridge vibration, which causes the pedestrian’s excitation being dependent on the bridge vibration; (2) pedestrian lateral excitation is a stochastic process rather than a perfect periodic load. Therefore, footbridge lateral vibration is essentially a complex nonlinear stochastic vibration system. Thus far, an effective method of dealing with such nonlinear stochastic vibration of footbridges remains lacking. A framework based on the probability density evolution (PDE) method is presented. For the mathematical model, the parameter resonance model is used to describe the pedestrian-bridge interaction while treating the pedestrian lateral excitation as a narrow-band process. For the analysis method, PDE is used to solve the nonlinear stochastic equations in combination with the number theoretical and finite difference methods. The proposed method establishes a new approach in studying footbridge lateral vibration. First, PDE based on the small sample strategy avoids the large amount of computation. Second, the randomness of both structural parameters and pedestrian lateral excitation could be taken into consideration by the proposed method. Third, based on the probability results with rich information, the serviceability, dynamic reliability, and random stability analyses are realized in a convenient manner.


---

## Body

## 1. Introduction

In 2000, the famous London Millennium Bridge was urgently closed down because of the large lateral vibration on its opening day [1]. This event became an influential symbol of the pedestrian-induced vibration of bridges and caused an in-depth discussion about the pedestrian-induced vibration of footbridges. Scholars began to study the underlying mechanisms, such as nonlinear vibration, parametric resonance, and pedestrian-bridge interaction, resulting in a variety of models. The first type is the linear model, such as the linear direct resonance model, which holds that pedestrian excitation is independent of bridge vibration. In this model, pedestrian lateral excitation is expressed by a simple harmonic function, and the large lateral vibration of footbridges is assumed to be caused by the direct resonance between pedestrian lateral walking frequency and bridge frequency. Fujino et al. [2] conducted a systematic study on the T-bridge in Japan by using this type of model. The second type is the vibration-dependent model, wherein pedestrian excitation depends on bridge vibration, and this dependency can be described by known models. A few dozen models are based on this type of method, such as the Dallard et al. [1], Nakamura [3], Piccardo and Tubino [4], and Ingólfsson et al. [5, 6] models. The Dallard, Nakamura, and Ingólfsson models belong to nonlinear velocity-dependent models. In this type of model, the velocity-dependent term is regarded as the application of additional negative damping to the structure. When the additional negative damping exceeds the real structural damping, vibration will diverge. The Piccardo model is a parametric resonance model based on the classical Matthew equation, by which the critical condition triggering vibration divergence can be obtained through stability analysis. The third type is a coupling model mixed with pedestrian motion and bridge vibration. This model also considers the influence of bridge vibration on pedestrian excitation, but this influence is described by the coupling equations of pedestrian motion and bridge vibration rather than a known empirical parameter model. Examples of coupling models are the Roberts [7], Newl [8], hybrid Van der Pol/Rayleigh [9, 10], and Macdonald et al. [11, 12] models. The Macdonald model, which is essentially an inverted pendulum model, considers the pedestrian lateral excitation as the inertia force of the human body mass center controlled by the balance strategy. Unlike other general models, this model insists that pedestrians maintain their comfort and balance by adjusting the position of the pace rather than the step frequency. Recently, on the basis of the inverted pendulum model proposed by Macdonald, Carroll et al. [13, 14] repeated the treadmill experiment of Ricciardelli and Pizzimenti [15] by using a 3D human motion capture technology and analyzed the pedestrian self-excitation characteristics caused by the pedestrian-structure interaction. Similarly, Bocian et al. [16], also based on the inverted pendulum model, carried out an experiment of pedestrian walking on the treadmill. The experiment was characterized by using a virtual reality simulation technology to reproduce the actual surrounding environment of pedestrians walking to reduce the influence of indoor environment on pedestrian.On the contrary, the randomness in pedestrian excitation is a factor that cannot be ignored because of its remarkable effect on bridge vibration, which may lead to results inconsistent with those of deterministic cases. However, research in this field is inadequate. The continuous pedestrian excitation obtained from actual observations is not an ideal periodic force but rather a narrow-band stochastic process. In the frequency domain, the Fourier spectra are not distributed in discrete frequency points obtained from the perfect periodic load but with a certain distributed width around the main harmonics, resulting in a reduced response, as compared with that of a perfect periodic force. At present, the analysis of the randomness of pedestrian excitation is focused on vertical direction [17–20]. By contrast, few investigations are available on the randomness of pedestrian lateral excitation. Ricciardelli and Pizzimenti [15] tested the lateral force caused by pedestrians walking on a treadmill. On the basis of the test data, they converted the Fourier spectrum to power spectrum and fitted the first five-order power spectrums by using the Gaussian shape. On the basis of Ricciardelli’s results, Ingólfsson and Georgakis [5] used a discrete-time Gaussian Markov process to simulate the random coefficients of equivalent damping and inertia forces; the body weight, walking frequency, step length, walking speed, and arrival time were considered random variables to represent the intersubject variability. The narrow-band process of pedestrian lateral excitation was further tested and studied by Racic and Brownjohn [21, 22]. They used multiple Gaussian curves to fit the lateral power spectrum and restored the amplitude and phase process. Bocian et al. [23] used the data of the British population to obtain the statistical values of pedestrian parameters (e.g., gait length, step frequency, and body weight) involved in the inverted pendulum model; they conducted a lateral stability analysis of the bridge in terms of probability. However, in Bocian’s model, the randomness of pedestrian excitation is actually not reflected through the stochastic process.Now, it well known that pedestrians act as autonomous dynamic systems interacting with the footbridge, and that the excitation term in the dynamic equation will contain the bridge response. This means that the pedestrian-induced lateral vibration is a nonlinear dynamic system. Moreover, as mentioned above, the random nature is another apparent characteristic of pedestrian load. The real consecutive pedestrian load is not a perfect periodic load but rather a narrow-band stochastic process. Therefore, the pedestrian-induced lateral vibration of footbridges is essentially a nonlinear stochastic vibration problem, which means that the traditional theory of linear stochastic vibration is no longer applicable. For nonlinear stochastic vibration problems, only a small number of cases have the closed-form solutions with certain specific conditions, but most of these problems in reality can only be solved by discrete numerical methods [5], which may require a large number of Monte Carlo (MC) numerical simulations.In the field of footbridge vibration, in addition to the aforementioned models and methods, postprocessing analysis (e.g., serviceability, dynamic reliability, and stability) is also the focus of researchers (especially designers). Serviceability is a crucial index to measure the normal performance of footbridges. When the number of people on the bridge is small, vibration could keep stable, but such vibration may affect the pedestrian walking and make the pedestrian feel uncomfortable and annoyed and result in the impairment of the bridge’s performance for normal usage. Therefore, it is necessary to investigate the serviceability of pedestrians when they walk on the footbridge. At present, most of the vibration serviceability assessment on the pedestrian structure are concentrated on the subjects of floor vertical vibration [24–26] and the footbridge vertical vibration [27–29], most of which belong to the linear vibration system. However, there are few reports about the serviceability of pedestrian impacting by the lateral vibration of the footbridge. Regarding the serviceability analysis of footbridges, the existing methods constantly use the concept of the peak acceleration or root mean square (RMS) acceleration reaching a certain threshold. In general, RMS acceleration is obtained directly from the acceleration power spectrum (or through the frequency response function (FRF) based on the force power spectrum) [17, 28]. This method is only suitable for linear systems, and the randomness of structural parameters is ignored.Serviceability analysis is actually a question of whether the bridge response (acceleration) exceeds a certain limit. Naturally, a deeper question arises: what is the probability of the bridge lateral vibration crossing a critical value? In other words, the probability of bridge lateral vibration remaining in a certain safe area must be ascertained. This process involves the issue of dynamic reliability. The dynamic reliability of bridges under various dynamic loads (e.g., seismic action, wind load, and train load) has been widely studied [30–35]. However, the dynamic reliability of pedestrian-induced lateral vibration in the footbridge has been rarely reported, especially when nonlinear stochastic vibration is considered. Unlike the conventional dynamic reliability problem, a nonlinear stochastic dynamic reliability problem is explored in this study. The methods for solving the nonlinear stochastic dynamic reliability problem can be divided into two types: numerical simulation method and analytical method. The Monte Carlo simulation method can be regarded as a simple and available numerical simulation method, but the huge amount of computation limits its applicability [5, 23, 36]. In terms of analytical methods, the conventional means for dynamic reliability of the nonlinear stochastic system usually need the using of Kolmogorov backward equation developed from the theory of the first-passage process [37–39]. This type of method requires the joint PDFs of the response process and its derivative process combined with a specific assumption of the crossing event, which makes it difficult to solve the multidimensional system. Moreover, this technique can only be used to analyze the single random (SR) case (i.e., the randomness of structural parameters cannot be considered).In this study, the probability density evolution (PDE) method based on the small sample strategy proposed by Li et al. [40–42] is used to avoid the large amount of computation caused by the MC numerical method. Meanwhile, the randomness in structural parameters and pedestrian lateral excitation are considered synthetically. The proposed method is successfully applied to the numerical examples (the London Millennium Bridge and the Passerelle Simone de Beauvoir bridge), and the probability results with rich information are obtained. Then, several critical analyses related to the footbridge’s normal service performance and ultimate limit performance, including the serviceability, dynamic reliability, and random stability analyses, are realized.The remainder of this paper is organized as follows: Section2 introduces the PDE method. Section 3 presents the nonlinear stochastic model for footbridge lateral vibration and describes the establishment and solution of PDE equations related to the nonlinear stochastic model. Section 4 presents the applications of the proposed method to the real footbridges and the realizations of serviceability, dynamic reliability, and random stability analyses. Finally, the last two sections draw the conclusions and discussions.


## 2. Probability Density Evolution Method

In practical engineering, apart from the randomness in the excitation load, the structural parameters, such as structural mass, modulus, and section properties, may also have strong randomness. Considering the randomness in structural parameters and excitation simultaneously (also known as a composite random issue) is a complex problem that has remained unsolved. The PDE method proposed by Li et al. [40–42] provides an alternative approach for conducting a composite random analysis.LetΘ=Θs,Θp denote all the random variables in the targeted system. Θs represents the random variables of the structural parameters, and Θp refers to the random variables from excitation. In general, a stochastic dynamic system can be expressed as follows:(1)Gy¨,y˙,y,Θp,t=FΘs,t.The structural response of equation (1) can be expressed as(2)y=yΘ,t.Given that the randomness ofΘ, y is a stochastic process. The velocity process of y is given as(3)y˙=∂yΘ,t∂t.The randomness of processy comes completely from Θ, and in the entire evolution process, no randomness disappears and no new randomness is added in the extended system [yt, Θ]. This feature indicates that [yt, Θ] is a probabilistic conservative system. On the basis of the probabilistic conservative principle, we have(4)ddt∫Ωt×ΩΘpyΘy,θ,tdydθ=0,where pyΘy,θ,t is the joint probability density function (PDF) of [yt, Θ] and Ωt×ΩΘ denotes the distribution domain of Θ at time t. After a series of derivations, the generalized density evolution equation can be obtained as(5)∂pyΘy,θ,t∂t+y˙θ,t∂pyΘy,θ,t∂y=0.


## 3. PDE for the Lateral Nonlinear Stochastic Vibration of Footbridges

During the opening day and subsequent on-site tests, a large lateral vibration was observed in the middle span of the London Millennium Bridge. The tests conducted by Arub Corp [1] demonstrated that the lateral fundamental frequency of the middle span is approximately 0.5 Hz, which is considerably lower than 1 Hz required for the pedestrian’s primary resonance; this result indicated that other mechanisms may exist. The parametric resonance model [4] can explain this phenomenon. In the present study, a parametric resonance stochastic model is established, in which the parametric resonance model is considered the basis, whereas the pedestrian lateral load is treated as a narrow-band stochastic excitation process.


### 3.1. Equations of Footbridge Lateral Vibration

Ignoring the moment of inertia, axial force, and shear deformation, the lateral vibration equation of the footbridge based on Euler–Bernoulli beam [43–45] can be expressed as follows:(6)msx∂2yx,t∂t2+csx∂yx,t∂t+∂2∂x2EIx∂2yx,t∂x2=fcx,t,where msx, EIx, csx, and yx,t are bridge mass per unit length, bending stiffness, damping, and lateral displacement along the spatial coordinate x and fcx,t denotes the lateral load per unit length exerted by pedestrians. Considering the proportional damping, equation (6) can be decoupled into the modal equations by Galerkin method. Only the first order yx,t=φxqt is considered here, in which the mode function is φx=sinπx/L with L being the bridge length. The modal equation is written as(7)q¨t+2ωsζq˙t+ωs2qt=1Ms∫0Lfcx,tφxdx=Ft,where qt, φx, ωs, ζ, Ms, and Ft are the modal displacement, modal shape, angular frequency, modal damping ratio, modal mass, and mass normalised modal load.


### 3.2. Pedestrian Lateral Excitation and Its Discretization

The crowd on the bridge can be divided into unsynchronized and synchronous groups. The lateral force produced by the unsynchronized group is assumed to be equivalent to the force produced by the crowd with different phases (uniform distribution) but with the same frequency, which is independent of vibration, and its number can be converted into the equivalent number established by the Fujino method [2]. The synchronous group is related to vibration, which can be quantified by the synchronization coefficient ρ. Hence, the pedestrian lateral excitation can be expressed as(8)fcx,t=1−ρLNmpxgd1+ρmpxgd2yx,tcosωpt,where mpx=Nmp/L is the pedestrian mass per unit length (N, pedestrian number; mp, single pedestrian mass), g is the acceleration of gravity, d1=0.04 and d2=2m−1 are the dynamic loading factors related to the static and self-excited loads, respectively, ρ is the synchronization coefficient, and ωp is the lateral walking frequency of the pedestrian.The randomness of the pedestrian load is now considered. As previously mentioned, the pedestrian load is not an ideal periodic force and is not suitable to be described in the form of a harmonic function. Hence, a stochastic excitation processξt is used to replace the harmonic function in equation (8), which can then be rewritten as(9)fcx,t=1−ρLNmpxgd1+ρmpxgd2yx,tξt.On the basis of equation (7), the modal load Ft is expressed as follows:(10)Ft=1Ms∫0L1−ρLNgd1mpxφxdx+∫0Lρgd2mpxφx2qdxξt.Letμ1=1−ρL/Ngd1/Ms, μ2=ρgd2/Ms, β1=∫0Lmpxφxdx, and β2=∫0Lmpxφx2dx; then, equation (10) can be written as(11)Ft=μ1β1+μ2β2qξt.To perform the numeric analysis in the time domain, the simulation and discretization ofξt are required. Various stochastic process simulation methods based on the power spectrum [46–48] have been widely used in engineering practice. Among them, the spectral expression method has been used frequently. On the basis of the spectral expression method, the pedestrian stochastic excitation process can be modeled as a harmonic function with random phases or frequencies as follows:(12)ξt≈∑i=1NAωicosωit+φi,Aωi=2SFωiΔωπ,where Aωi, ωi, and φi are the amplitude, circular frequency, and phase of the i-th harmonic component, Δω is the increment of ωi, N is the number of harmonic components, and SFω is the PSD of ξt. On the basis of the experimental results of Ricciardelli and Pizzimenti [15], the PSD around the first harmonic of pedestrian lateral excitation ξt can be expressed as(13)SFω=σF22πωasbsexp−2ω/ωp−1bs2,where σF2=W2d2 (W, single pedestrian weight; d, dynamic loading factor) denotes the doubled area of PSD around the first load harmonic and ωp is the reference undisturbed gait frequency. as=0.9 and bs=0.043 are the fitting parameters.In the normal spectral expression method, the harmonic function phase is assumed to be a uniformly distributed random variable to simulate the randomness of excitation, and the harmonic function frequency is a series of divided deterministic variables (ωi=iΔω,i=1,2,…,N). To guarantee the accuracy, the above normal spectral expression method often requires at least 500 harmonic components, resulting in high-dimensional random vectors and large computational workload. Fortunately, an improved approach named the second kind of spectral expression (SSE) could significantly reduce the discrete harmonic components. On the basis of a previous study [49], if Aωi, ωi, and φi in equation (12) satisfy the following conditions:(1)
φi,i=1,2,…,N are mutually independent random variables subject to uniform distribution among 0,2π,(2)
ωi,i=1,2,…,N are mutually independent random variables subject to uniform distribution among ωi−1s,ωis, i=1,2,…,N, respectively,(3)
Aωi=2SFωiΔωi/π, Δωi=ωis−ωi−1s, i=1,2,…,N,then the PSD of harmonic function ξt represented by equation (12) will be equal to SFω. The use of this method can remarkably reduce the number of harmonic components needed for discretization.According to the expression of PSD given in equation (13), the SSE is applied for the simulation and discretization of ξt. The analysis results (Figure 1) show that the accuracy requirements will be met in the SSE if the number of harmonic components is up to N=10 (i.e., the total number of random variables is 2N=20).Figure 1
Simulation ofξt by using the SSE method: (a) one discrete sample generated by SSE; (b) PSDs under different discrete harmonic components and the comparison with real PSD.
(a)
(b)


### 3.3. Establishing the PDE for Lateral Vibration of the Footbridge

On the basis of equations (5) and (7), the PDE for the lateral vibration of the footbridge can be expressed as(14)∂pqΘq,θ,t∂t+q˙θ,t∂pqΘq,θ,t∂q=0.The corresponding initial condition is(15)pqΘq,θ,t0=δq−q0pΘθ,where q0 is the deterministic initial value for qt and δ is the Dirac delta function. By solving equation (14) combined with the initial condition of equation (15), the PDF of structural response q can be obtained as(16)pqq,t=∫ΩΘpqΘq,θ,tdθ.The value ofq˙θ,t should be determined firstly before solving equation (14). It is observed that there is no partial differential in terms of θ in equation (14). This means that for specific θ, equation (14) can be solved by a deterministic numerical method (in this study, the finite difference method will be adopted, see Section 3.5). Therefore, it is necessary to conduct the discretization of θ, whereby a set of discrete representative points θi,i=1,2,…,n in the probability distribution space of θ can be obtained. Once these discrete representative points have been determined, the velocity process of q˙θ,t can be obtained by the deterministic solution of equation (7). Afterwards, piq,θi,t can also be obtained according to equation (14), and equation (16) will be changed to(17)pqq,t=∑i=1npiq,θi,t.The accuracy of equation (17) depends on the discrete representative point θi,i=1,2,…,n. It is not a wise choice to select a sufficient number of θ to ensure the accuracy, as this may result in a huge computation. Therefore, how to select a set of representative points in the probability distribution space of θ to meet the accuracy requirements while reducing the number of points as small as possible will be the next problem to be considered. This will involve the partition of the probability-assigned space and selection of representative points, which are discussed in following Section 3.4.


### 3.4. Partition of the Probability-Assigned Space and Selection of Representative Points

For any subdomainΩh, h=1,2,…,np, which is a partition of domain space ΩΘ, if ∪h=1npΩh=ΩΘ and Ωh∩​Ωwh≠w=∅, then equation (16) can then be rewritten as(18)pqq,t=∑h=1np∫ΩhpqΘq,θ,tdθ=∑h=1npphq,θh,t,where phq,θh,t=∫ΩhpqΘq,θ,tdθ. Let(19)Ph=∫−∞∞phq,θh,tdq=∫−∞∞∫ΩhpqΘq,θ,tdθdq=∫ΩhpΘθdθ.Hence,Ph is a normalised factor and also called the assigned probability because it denotes the probability in subdomain Ωh.By selecting a representative pointθh in subdomain Ωh, equation (14) becomes(20)∂phq,θh,t∂t+q˙θh,t∂phq,θh,t∂q=0,h=1,2,…,np.The corresponding initial condition can be rewritten as(21)phq,θh,t0=δq−q0Ph.By combining equations (20) and (21), a group of equations with number np can be solved to obtain phq,t. Then, the PDF of q can be easily obtained by using equation (18). The representative point θh and assigned probability Ph depend strongly on the partition of the probability-assigned space (i.e., the selection of subdomain Ωh). For a given representative point set, its representative subdomain is preferred to the Voronoi region (Figure 1).As mentioned earlier, for the PDE analysis of structural response, it is necessary to obtain discrete representative points in multidimensional probability spaceΩΘ, and these discrete representative points should be scattered as uniformly as possible. When the dimension of random variables is large, the number theoretical method (NTM) can be used, by which a representative point set uniformly scattered can be generated.In the NTM, the discrepancy is used to measure the uniformity of a set of sample points. For as-dimensional space, let Cs denote a unit cube in s-dimensional space. Furthermore, let n be a positive integer and θk=θ1nk,θ2nk,…,θsnk,1≤k≤n denote a set of sample points of Cs, among which, 0≤θink<1,1≤i≤s. For any η=η1,η2,…,ηs∈Cs, if set Nnη,θ as the number of points which meet the condition 0≤θk,1n<η1,…,0≤θk,sn<ηs, the discrepancy Dn,θ of θk is given as(22)Dn,θ=supς∈CsNnη,θn−η,where η=η1η2…ηs denotes the volume of region 0,η. And when the following condition holds(23)Dn,θ=On−1/2,where O denotes the order of error, then θ can be considered as the NTM point set among Cs. There are several methods that can be used to generate the NTM point set, such as good Latin point (GLP) set, Halton set, Hamersley set, and good point (GP) set. The GP set will be used in this paper, as it has a good applicability to the high-dimensional problem. Let η=η1,η2,…,ηs∈Cs and suppose that there is a set θk expressed as(24)θk=η1k,…,ηsk,1≤k≤n.If the discrepancy ofθk meets(25)Dn,θ≤Oη,εn−1+ε.Thenθ will be called a GP set. The operational character . in equation (24) means the decimal operator. Normally, the good point η can be generated by the cyclotomic field method [50]:(26)η=2cos2πγ,2cos4πγ,…,2cos2πsγ,where γ is a prime number which satisfies γ≥2s+3. By substituting equation (26) into (24), the wanted set of GP points can be obtained. To better understand this method, two simple cases are considered:(1)
Two variables (x1,x2) use the NTM to generate the GP points among the intervals of x1∈−66 and x2∈−412(2)
Three variables (x1,x2,x3) use the NTM to generate the GP points among the intervals of x1∈−66, x2∈−412,andx3∈515.Figure2 presents the GP points for the two cases. It can be observed that the GP points have a good uniformity.Figure 2
GP points (green “o” marker) and Voronoi regions (polyhedral). (a) 2 variables; (b) 3 variables.
(a)
(b)


### 3.5. Finite Difference Method (Total Variation Diminishing (TVD) Scheme Based on the Lax–Wendroff Method)

Equation (20) is a typical convective equation, which has difficulty in obtaining its analytical solution. The finite difference method [51, 52] or the finite element method [53] can be used to solve such partial differential equations. For simplicity, the finite difference method is adopted in this study. To ensure the second-order accuracy and avoid oscillations, the TVD scheme [54] based on the Lax–Wendroff scheme is used in this study. For convenience, phq,θh,t in equation (20) will be abbreviated to p. Initially, the Lax–Wendroff expression of equation (20) is written as(27)pjn+1=pjn−12q˙λ+q˙λΔpj−1/2n−12q˙λ−q˙λΔpj+1/2n−121−q˙λq˙λΔpj+1/2n−Δpj−1/2n,where n denotes the time mesh step index, j represents the space mesh step index, λ=Δt/Δy refers to the mesh ratio between the time step and the space step, q˙ refers to the velocity process term (i.e., q˙θh,t), and Δpj−1/2n=pjn−pj−1n; Δpj+1/2n=pj+1n−pjn. The first and second terms in equation (27) are the numerical fluxes of the upwind and Lax–Wendroff schemes, respectively. To ensure that equation (27) has the TVD scheme (avoid oscillations) and second-order accuracy, equation (27) is rewritten as(28)pjn+1=pjn−12q˙λ+q˙λΔpj−1/2n−12q˙λ−q˙λΔpj+1/2n−121−q˙λq˙λφrj+1/2+,rj+1/2−Δpj+1/2n−φrj−1/2+,rj−1/2−Δpj−1/2n,where rj+1/2+=Δpj+3/2n/Δpj+1/2n, rj+1/2−=Δpj−1/2n/Δpj+1/2n, rj−1/2+=Δpj+1/2n/Δpj−1/2n, and rj−1/2−=Δpj−3/2n/Δpj−1/2n and φ is a function of r+ and r−, which is called the flux limiter. To ensure that equation (28) satisfies the sufficient conditions of TVD, φ must be further restricted. In this study, the following adaptive flux limiter will be used [55]:(29)φr+,r−=ϑ−q˙φsr++ϑq˙φsr−,φsr=max0,min2r,1,minr,2,where ϑx=1,x≥00,x<0. The Courant–Friedrichs–Lewy (C.F.L) criteria for the stability of equation (28) is given by(30)q˙maxλ≤1.The value ofq˙max is estimated by trial calculations, by which the mesh ratio λ can be determined according to equation (30). By doing this, the convergence of the solution of equation (28) will be guaranteed.To illustrate this method, a simple numerical example is considered:(31)∂p∂t+a∂p∂x=0,x∈−∞,∞,t>0,where a=2, and the initial condition px,0 is given as(32)px,0=10x+1,−0.1≤x≤0,−10x+1,0≤x≤0.1,0,others.For comparison, besides the method of TVD based on Lax–Wendroff used in this study, two other methods are also used to solve equation (31), including the upwind method, the Lax-Wendroff method. Figure 3(a) presents the solution at t=0.2, as well as the comparative results among different methods. It can be observed that the upwind method with the inherent accuracy up to first order has the worst performance; the Lax–Wendroff method has an oscillation in the vicinity of the left-hand side due to the effect of dispersion, while the method of TVD based on Lax–Wendroff has nearly the same result as the exact solution. Figure 3(b) presents the entire process of the solution along the timeline obtained by the method of TVD based on Lax–Wendroff, from which the convergence can be confirmed.Figure 3
Solutions of equation (31) at t=0.2. (a) Comparative results among different methods. (b) Entire process of solution by the method of TVD based on Lax–Wendroff.
(a)
(b)


### 3.6. Outline of the Procedure

The proposed method is presented in Figure4, which is based on the following iterative procedure.Step 1. The random variables of structural parameters Θs=θs1,θs2,…,θsn, such as bending stiffness, mass, and damping, are determined. Meanwhile, the SSE method is used to discretize the pedestrian lateral excitation process ξt into harmonic components with N=10, by which the discrete random variables of Θp=θω1,θω2,…,θω10;θφ1,θφ2,…,θφ10 with 2N=20 are defined. θφi,i=1,2,…,N are mutually independent random variables subject to uniform distribution among 0,2π, and θωi,i=1,2,…,N are mutually independent random variables subject to uniform distribution among ωi−1s,ωis, i=1,2,…,N, respectively.Step 2. NTM is used to select the set of discrete representative points θh,h=1,2,…,np in space ΩΘ, where Θ=Θs,Θp. np is generally set to 10 times as the total number of random variables, that is, np=10dimΘs+dimΘs. Meanwhile, the corresponding assigned probability Ph,q=1,2,…,np is determined through the statistical characteristics of each variable and the Voronoi partition region.Step 3. For each θh, the deterministic solution of lateral vibration, equation (7), is performed to yield the velocity process q˙θh,t.Step 4. q˙θh,t obtained in Step 3 is substituted into equation (20) and combined with the initial conditional equation (21). Then, equation (20) is calculated to obtain the numerical solution phy,t,h=1,2,…,np by using the finite difference method (i.e., equation (28)).Step 5. Finally, pqq,t can be obtained by substituting the results from Step 4 into equation (18) and performing a summation.Figure 4
Flowchart of the proposed method.


## 3.1. Equations of Footbridge Lateral Vibration

Ignoring the moment of inertia, axial force, and shear deformation, the lateral vibration equation of the footbridge based on Euler–Bernoulli beam [43–45] can be expressed as follows:(6)msx∂2yx,t∂t2+csx∂yx,t∂t+∂2∂x2EIx∂2yx,t∂x2=fcx,t,where msx, EIx, csx, and yx,t are bridge mass per unit length, bending stiffness, damping, and lateral displacement along the spatial coordinate x and fcx,t denotes the lateral load per unit length exerted by pedestrians. Considering the proportional damping, equation (6) can be decoupled into the modal equations by Galerkin method. Only the first order yx,t=φxqt is considered here, in which the mode function is φx=sinπx/L with L being the bridge length. The modal equation is written as(7)q¨t+2ωsζq˙t+ωs2qt=1Ms∫0Lfcx,tφxdx=Ft,where qt, φx, ωs, ζ, Ms, and Ft are the modal displacement, modal shape, angular frequency, modal damping ratio, modal mass, and mass normalised modal load.


## 3.2. Pedestrian Lateral Excitation and Its Discretization

The crowd on the bridge can be divided into unsynchronized and synchronous groups. The lateral force produced by the unsynchronized group is assumed to be equivalent to the force produced by the crowd with different phases (uniform distribution) but with the same frequency, which is independent of vibration, and its number can be converted into the equivalent number established by the Fujino method [2]. The synchronous group is related to vibration, which can be quantified by the synchronization coefficient ρ. Hence, the pedestrian lateral excitation can be expressed as(8)fcx,t=1−ρLNmpxgd1+ρmpxgd2yx,tcosωpt,where mpx=Nmp/L is the pedestrian mass per unit length (N, pedestrian number; mp, single pedestrian mass), g is the acceleration of gravity, d1=0.04 and d2=2m−1 are the dynamic loading factors related to the static and self-excited loads, respectively, ρ is the synchronization coefficient, and ωp is the lateral walking frequency of the pedestrian.The randomness of the pedestrian load is now considered. As previously mentioned, the pedestrian load is not an ideal periodic force and is not suitable to be described in the form of a harmonic function. Hence, a stochastic excitation processξt is used to replace the harmonic function in equation (8), which can then be rewritten as(9)fcx,t=1−ρLNmpxgd1+ρmpxgd2yx,tξt.On the basis of equation (7), the modal load Ft is expressed as follows:(10)Ft=1Ms∫0L1−ρLNgd1mpxφxdx+∫0Lρgd2mpxφx2qdxξt.Letμ1=1−ρL/Ngd1/Ms, μ2=ρgd2/Ms, β1=∫0Lmpxφxdx, and β2=∫0Lmpxφx2dx; then, equation (10) can be written as(11)Ft=μ1β1+μ2β2qξt.To perform the numeric analysis in the time domain, the simulation and discretization ofξt are required. Various stochastic process simulation methods based on the power spectrum [46–48] have been widely used in engineering practice. Among them, the spectral expression method has been used frequently. On the basis of the spectral expression method, the pedestrian stochastic excitation process can be modeled as a harmonic function with random phases or frequencies as follows:(12)ξt≈∑i=1NAωicosωit+φi,Aωi=2SFωiΔωπ,where Aωi, ωi, and φi are the amplitude, circular frequency, and phase of the i-th harmonic component, Δω is the increment of ωi, N is the number of harmonic components, and SFω is the PSD of ξt. On the basis of the experimental results of Ricciardelli and Pizzimenti [15], the PSD around the first harmonic of pedestrian lateral excitation ξt can be expressed as(13)SFω=σF22πωasbsexp−2ω/ωp−1bs2,where σF2=W2d2 (W, single pedestrian weight; d, dynamic loading factor) denotes the doubled area of PSD around the first load harmonic and ωp is the reference undisturbed gait frequency. as=0.9 and bs=0.043 are the fitting parameters.In the normal spectral expression method, the harmonic function phase is assumed to be a uniformly distributed random variable to simulate the randomness of excitation, and the harmonic function frequency is a series of divided deterministic variables (ωi=iΔω,i=1,2,…,N). To guarantee the accuracy, the above normal spectral expression method often requires at least 500 harmonic components, resulting in high-dimensional random vectors and large computational workload. Fortunately, an improved approach named the second kind of spectral expression (SSE) could significantly reduce the discrete harmonic components. On the basis of a previous study [49], if Aωi, ωi, and φi in equation (12) satisfy the following conditions:(1)
φi,i=1,2,…,N are mutually independent random variables subject to uniform distribution among 0,2π,(2)
ωi,i=1,2,…,N are mutually independent random variables subject to uniform distribution among ωi−1s,ωis, i=1,2,…,N, respectively,(3)
Aωi=2SFωiΔωi/π, Δωi=ωis−ωi−1s, i=1,2,…,N,then the PSD of harmonic function ξt represented by equation (12) will be equal to SFω. The use of this method can remarkably reduce the number of harmonic components needed for discretization.According to the expression of PSD given in equation (13), the SSE is applied for the simulation and discretization of ξt. The analysis results (Figure 1) show that the accuracy requirements will be met in the SSE if the number of harmonic components is up to N=10 (i.e., the total number of random variables is 2N=20).Figure 1
Simulation ofξt by using the SSE method: (a) one discrete sample generated by SSE; (b) PSDs under different discrete harmonic components and the comparison with real PSD.
(a)
(b)


## 3.3. Establishing the PDE for Lateral Vibration of the Footbridge

On the basis of equations (5) and (7), the PDE for the lateral vibration of the footbridge can be expressed as(14)∂pqΘq,θ,t∂t+q˙θ,t∂pqΘq,θ,t∂q=0.The corresponding initial condition is(15)pqΘq,θ,t0=δq−q0pΘθ,where q0 is the deterministic initial value for qt and δ is the Dirac delta function. By solving equation (14) combined with the initial condition of equation (15), the PDF of structural response q can be obtained as(16)pqq,t=∫ΩΘpqΘq,θ,tdθ.The value ofq˙θ,t should be determined firstly before solving equation (14). It is observed that there is no partial differential in terms of θ in equation (14). This means that for specific θ, equation (14) can be solved by a deterministic numerical method (in this study, the finite difference method will be adopted, see Section 3.5). Therefore, it is necessary to conduct the discretization of θ, whereby a set of discrete representative points θi,i=1,2,…,n in the probability distribution space of θ can be obtained. Once these discrete representative points have been determined, the velocity process of q˙θ,t can be obtained by the deterministic solution of equation (7). Afterwards, piq,θi,t can also be obtained according to equation (14), and equation (16) will be changed to(17)pqq,t=∑i=1npiq,θi,t.The accuracy of equation (17) depends on the discrete representative point θi,i=1,2,…,n. It is not a wise choice to select a sufficient number of θ to ensure the accuracy, as this may result in a huge computation. Therefore, how to select a set of representative points in the probability distribution space of θ to meet the accuracy requirements while reducing the number of points as small as possible will be the next problem to be considered. This will involve the partition of the probability-assigned space and selection of representative points, which are discussed in following Section 3.4.


## 3.4. Partition of the Probability-Assigned Space and Selection of Representative Points

For any subdomainΩh, h=1,2,…,np, which is a partition of domain space ΩΘ, if ∪h=1npΩh=ΩΘ and Ωh∩​Ωwh≠w=∅, then equation (16) can then be rewritten as(18)pqq,t=∑h=1np∫ΩhpqΘq,θ,tdθ=∑h=1npphq,θh,t,where phq,θh,t=∫ΩhpqΘq,θ,tdθ. Let(19)Ph=∫−∞∞phq,θh,tdq=∫−∞∞∫ΩhpqΘq,θ,tdθdq=∫ΩhpΘθdθ.Hence,Ph is a normalised factor and also called the assigned probability because it denotes the probability in subdomain Ωh.By selecting a representative pointθh in subdomain Ωh, equation (14) becomes(20)∂phq,θh,t∂t+q˙θh,t∂phq,θh,t∂q=0,h=1,2,…,np.The corresponding initial condition can be rewritten as(21)phq,θh,t0=δq−q0Ph.By combining equations (20) and (21), a group of equations with number np can be solved to obtain phq,t. Then, the PDF of q can be easily obtained by using equation (18). The representative point θh and assigned probability Ph depend strongly on the partition of the probability-assigned space (i.e., the selection of subdomain Ωh). For a given representative point set, its representative subdomain is preferred to the Voronoi region (Figure 1).As mentioned earlier, for the PDE analysis of structural response, it is necessary to obtain discrete representative points in multidimensional probability spaceΩΘ, and these discrete representative points should be scattered as uniformly as possible. When the dimension of random variables is large, the number theoretical method (NTM) can be used, by which a representative point set uniformly scattered can be generated.In the NTM, the discrepancy is used to measure the uniformity of a set of sample points. For as-dimensional space, let Cs denote a unit cube in s-dimensional space. Furthermore, let n be a positive integer and θk=θ1nk,θ2nk,…,θsnk,1≤k≤n denote a set of sample points of Cs, among which, 0≤θink<1,1≤i≤s. For any η=η1,η2,…,ηs∈Cs, if set Nnη,θ as the number of points which meet the condition 0≤θk,1n<η1,…,0≤θk,sn<ηs, the discrepancy Dn,θ of θk is given as(22)Dn,θ=supς∈CsNnη,θn−η,where η=η1η2…ηs denotes the volume of region 0,η. And when the following condition holds(23)Dn,θ=On−1/2,where O denotes the order of error, then θ can be considered as the NTM point set among Cs. There are several methods that can be used to generate the NTM point set, such as good Latin point (GLP) set, Halton set, Hamersley set, and good point (GP) set. The GP set will be used in this paper, as it has a good applicability to the high-dimensional problem. Let η=η1,η2,…,ηs∈Cs and suppose that there is a set θk expressed as(24)θk=η1k,…,ηsk,1≤k≤n.If the discrepancy ofθk meets(25)Dn,θ≤Oη,εn−1+ε.Thenθ will be called a GP set. The operational character . in equation (24) means the decimal operator. Normally, the good point η can be generated by the cyclotomic field method [50]:(26)η=2cos2πγ,2cos4πγ,…,2cos2πsγ,where γ is a prime number which satisfies γ≥2s+3. By substituting equation (26) into (24), the wanted set of GP points can be obtained. To better understand this method, two simple cases are considered:(1)
Two variables (x1,x2) use the NTM to generate the GP points among the intervals of x1∈−66 and x2∈−412(2)
Three variables (x1,x2,x3) use the NTM to generate the GP points among the intervals of x1∈−66, x2∈−412,andx3∈515.Figure2 presents the GP points for the two cases. It can be observed that the GP points have a good uniformity.Figure 2
GP points (green “o” marker) and Voronoi regions (polyhedral). (a) 2 variables; (b) 3 variables.
(a)
(b)


## 3.5. Finite Difference Method (Total Variation Diminishing (TVD) Scheme Based on the Lax–Wendroff Method)

Equation (20) is a typical convective equation, which has difficulty in obtaining its analytical solution. The finite difference method [51, 52] or the finite element method [53] can be used to solve such partial differential equations. For simplicity, the finite difference method is adopted in this study. To ensure the second-order accuracy and avoid oscillations, the TVD scheme [54] based on the Lax–Wendroff scheme is used in this study. For convenience, phq,θh,t in equation (20) will be abbreviated to p. Initially, the Lax–Wendroff expression of equation (20) is written as(27)pjn+1=pjn−12q˙λ+q˙λΔpj−1/2n−12q˙λ−q˙λΔpj+1/2n−121−q˙λq˙λΔpj+1/2n−Δpj−1/2n,where n denotes the time mesh step index, j represents the space mesh step index, λ=Δt/Δy refers to the mesh ratio between the time step and the space step, q˙ refers to the velocity process term (i.e., q˙θh,t), and Δpj−1/2n=pjn−pj−1n; Δpj+1/2n=pj+1n−pjn. The first and second terms in equation (27) are the numerical fluxes of the upwind and Lax–Wendroff schemes, respectively. To ensure that equation (27) has the TVD scheme (avoid oscillations) and second-order accuracy, equation (27) is rewritten as(28)pjn+1=pjn−12q˙λ+q˙λΔpj−1/2n−12q˙λ−q˙λΔpj+1/2n−121−q˙λq˙λφrj+1/2+,rj+1/2−Δpj+1/2n−φrj−1/2+,rj−1/2−Δpj−1/2n,where rj+1/2+=Δpj+3/2n/Δpj+1/2n, rj+1/2−=Δpj−1/2n/Δpj+1/2n, rj−1/2+=Δpj+1/2n/Δpj−1/2n, and rj−1/2−=Δpj−3/2n/Δpj−1/2n and φ is a function of r+ and r−, which is called the flux limiter. To ensure that equation (28) satisfies the sufficient conditions of TVD, φ must be further restricted. In this study, the following adaptive flux limiter will be used [55]:(29)φr+,r−=ϑ−q˙φsr++ϑq˙φsr−,φsr=max0,min2r,1,minr,2,where ϑx=1,x≥00,x<0. The Courant–Friedrichs–Lewy (C.F.L) criteria for the stability of equation (28) is given by(30)q˙maxλ≤1.The value ofq˙max is estimated by trial calculations, by which the mesh ratio λ can be determined according to equation (30). By doing this, the convergence of the solution of equation (28) will be guaranteed.To illustrate this method, a simple numerical example is considered:(31)∂p∂t+a∂p∂x=0,x∈−∞,∞,t>0,where a=2, and the initial condition px,0 is given as(32)px,0=10x+1,−0.1≤x≤0,−10x+1,0≤x≤0.1,0,others.For comparison, besides the method of TVD based on Lax–Wendroff used in this study, two other methods are also used to solve equation (31), including the upwind method, the Lax-Wendroff method. Figure 3(a) presents the solution at t=0.2, as well as the comparative results among different methods. It can be observed that the upwind method with the inherent accuracy up to first order has the worst performance; the Lax–Wendroff method has an oscillation in the vicinity of the left-hand side due to the effect of dispersion, while the method of TVD based on Lax–Wendroff has nearly the same result as the exact solution. Figure 3(b) presents the entire process of the solution along the timeline obtained by the method of TVD based on Lax–Wendroff, from which the convergence can be confirmed.Figure 3
Solutions of equation (31) at t=0.2. (a) Comparative results among different methods. (b) Entire process of solution by the method of TVD based on Lax–Wendroff.
(a)
(b)


## 3.6. Outline of the Procedure

The proposed method is presented in Figure4, which is based on the following iterative procedure.Step 1. The random variables of structural parameters Θs=θs1,θs2,…,θsn, such as bending stiffness, mass, and damping, are determined. Meanwhile, the SSE method is used to discretize the pedestrian lateral excitation process ξt into harmonic components with N=10, by which the discrete random variables of Θp=θω1,θω2,…,θω10;θφ1,θφ2,…,θφ10 with 2N=20 are defined. θφi,i=1,2,…,N are mutually independent random variables subject to uniform distribution among 0,2π, and θωi,i=1,2,…,N are mutually independent random variables subject to uniform distribution among ωi−1s,ωis, i=1,2,…,N, respectively.Step 2. NTM is used to select the set of discrete representative points θh,h=1,2,…,np in space ΩΘ, where Θ=Θs,Θp. np is generally set to 10 times as the total number of random variables, that is, np=10dimΘs+dimΘs. Meanwhile, the corresponding assigned probability Ph,q=1,2,…,np is determined through the statistical characteristics of each variable and the Voronoi partition region.Step 3. For each θh, the deterministic solution of lateral vibration, equation (7), is performed to yield the velocity process q˙θh,t.Step 4. q˙θh,t obtained in Step 3 is substituted into equation (20) and combined with the initial conditional equation (21). Then, equation (20) is calculated to obtain the numerical solution phy,t,h=1,2,…,np by using the finite difference method (i.e., equation (28)).Step 5. Finally, pqq,t can be obtained by substituting the results from Step 4 into equation (18) and performing a summation.Figure 4
Flowchart of the proposed method.


## 4. Case Study

The London Millennium Bridge is the first footbridge built on the Thames River. This structure is a shallow suspension bridge with the spans of 81 m + 144 m + 108 m. During its opening day and subsequent on-site tests, large lateral vibrations in all three spans were observed. The middle span of the Millennium Bridge is selected as the main analysis object, wherein the fundamental frequency and damping ratio are approximately 0.48 Hz and 0.007, respectively. The bending stiffness of the main beamEIb, mass per unit length mb, and damping ratio ζ are considered as structural random variables. Since the report of the statistical characteristics of these structural random variables are very limited, it is assumed that they obey the normal distribution, and their variability are defined according to the general practical experience. The properties of structural random variables are listed in Table 1. The detail of the discrete random variables from pedestrian lateral excitation can be referred to Section 3.2.Table 1
Properties of structural random variables.
Random variable
Distribution type
Mean
Coefficient of variation
Bending stiffness of the main beam—EIb, (N⋅m2)
Normal
8.0383 × 1010
0.05
Mass per unit length—mb, (kg/m)
Normal
2000
0.05
Damping ratio—ζ
Normal
0.007
0.1Figure5 shows a series of probability results of the lateral vibration displacement of the bridge when the number of pedestrians on the bridge is set to Np=120, and the pedestrian central lateral walking frequency is the double of the bridge fundamental frequency (parameter resonance). Figure 5(a) presents the 3D diagram of PDFs during evolution containing the entire probability information of the timeline. To achieve an improved display effect, Figure 5(b) illustrates the PDF contour results (i.e., density flow) in the entire process. The results show that the density flow moves steadily without abrupt changes when the number of pedestrians is relatively small. It can also be found that, with the increase of time, the probability density begins to increase in the large displacement areas but tends to be stable at the later stage. On the basis of the results of Figure 5(a), some certain time points of interest can be obtained to analyze the specified PDF of lateral displacement. The PDF results of three time points (T/3, 2T/3, and T) are shown in Figure 5(c) and compared with the MC results. In the Monte Carlo simulation (MCS), the structural random variables and discrete random variables are sampled according to their distributions, and the dynamic equation of equation (7) is solved under each set of samples, by which the statistic of structural response is implemented to obtain the probability distributions. In Figure 5(c), it can be known that the PDF results at different time instances of the proposed method are consistent with those of MCS. In Figure 5(d) which gives the mean values of displacement along the timeline, the proposed method has a worse performance in the initial stage, but its results closely match with those of MCS in the later stage. The congruous comparative results confirm the effectiveness of the proposed method. Compared with the 10,000 repetitions of dynamic analyses needed in the MCS, the PDE-based approach only needs 230 iterations without compromising the accuracy.Figure 5
Probability results of the lateral displacement (Np=120). (a) PDF evolution process (T = 30 s). (b) Contour of the PDF surface during evolution. (c) PDFs at different time instances. (d) Mean values of displacement.
(a)
(b)
(c)
(d)To further illustrate the availability of the proposed method, the Passerelle Simone de Beauvoir bridge [56], a combined shallow arch bridge located in Paris, is also taken as the analysis object. Its main span length, bending stiffness of the main beam, and mass per unit length are 190 m, 5.665×1011N⋅m2, and 3420 kg/m. Its fundamental frequency and damping ratio are approximately 0.56 Hz and 0.0076, respectively. Like the London Millennium Bridge, the bending stiffness of the main beam, mass per unit length, and damping ratio are considered as structural random variables, and their coefficients of variation are assumed as the same as those in the London Millennium Bridge. Due to space limitations, Figure 6 presents the results of PDFs at different time instances. The results show good agreement with the MCS results, confirming the availability of the proposed method again.Figure 6
PDFs at different time instances (P-bridge).


### 4.1. Influence of Randomness

In comparison with the single random (SR) case that considers the randomness of excitation only, the composite random (CR) case (considers the randomness of the structural parameters and excitation) has always been a difficult issue. To date, no studies are available on the CR analysis of footbridge lateral vibration. The structural parameters of the footbridge may have evident randomness that also considerably influences the results due to the measurement error. Figures7(a) and 7(b) show that the PDF result of the CR case is flatter than that of the SR case, indicating that the former dispersion is higher than the latter. This finding is consistent with the expectation because adding the randomness of structural parameters to the system will likely lead to a high discreteness of the response. Given that only three structural parameters are considered random variables in this study, it is believed that the results will be more different from those of the SR case if more structural parameters are involved. Furthermore, to analyze the influence brought by the random disturbance intensity of excitation, the parameter bs in equation (13) is slightly changed to be 0.35. From Figure 7(c), when bs is reduced (i.e., the random disturbance intensity is increased), the PDF around the zero-value region grows larger, which means that increasing the random disturbance intensity is helpful to keep the small vibration.Figure 7
Influence of randomness underNp=120. (a) Comparative PDF results of the CR and SR cases. (b) Comparative CDF results of the CR and SR cases. (c) Comparative PDF results of different bs (under CR case).
(a)
(b)
(c)


### 4.2. Serviceability and Dynamic Reliability

In this study, the RMS acceleration is obtained in a convenient manner of using all the temporal and spatial probability distributions instead of using the power spectrum method. By replacing the displacementqt with acceleration q¨t in equation (2), the probability distributions of acceleration pq¨q¨,t can be obtained by using the same procedure as that for pqq,t. Then, RMS acceleration can be developed by using the following expressions:(33)Γq¨t=Eq¨t2=∫−∞∞q¨2pq¨q¨,tdq¨,σq¨=1T∫0TΓq¨tdt=1T∫0TEq¨t2dt,where Γq¨t is the original moment of second order.In comparison with the definition of the vertical comfort limit on the footbridge, limited references are available for the lateral direction. Bachmann et al. [57] suggested that the limit of lateral peak acceleration should be no more than 0.2 m/s2, and Dallard et al. [58] thought that the limit set to 0.2–0.4 m/s2 was more appropriate. In terms of code, BS 5400 [59] offered that the maximum lateral acceleration should be less than 0.25 m/s2, and EN 1990 [60] suggested that lateral peak acceleration should be less than 0.14f or 0.15 m/s2. In the present study, a strict limit is adopted, wherein the RMS acceleration may not exceed the value of 0.12fs.The time length is still set to 30 s. Figure8(a) shows the time-history results of Γq¨t1/2 when the number of pedestrians Np are 90, 110, 130, and 150. Figure 8(b) presents the results of σq¨ when the number of pedestrians ranges from 90 to 150. As shown in Figure 8(a), Γq¨t1/2 has a low value, and the corresponding RMS acceleration σq¨ is considerably lower than the limit value of comfort (Figure 8(b)) when the number of pedestrians is small. However, σq¨ rapidly increases when the number of pedestrians increases. When the number of people is up to Np=150, σq¨ exceeds the comfort limit, indicating that the pedestrians will start to feel annoyed due to the bridge vibration. If the number of pedestrians continues to increase to a certain critical value, then a large instability vibration of the bridge may occur. At this time, the pedestrian’s feelings may be beyond the concept of comfort but probably panic. Hence, the analysis will no longer belong to the scope of serviceability but the random stability analysis of bridge vibration.Figure 8
Γ
q
¨
t
1
/
2 and σq¨ under different numbers of pedestrians. (a) Γq¨t1/2vs.T under different Np; (b) σq¨ vs. Np.
(a)
(b)As previously mentioned, the conventional analytical methods are mainly based on the Kolmogorov backward equations, which are difficult to apply due to their strict assumptions and inability to deal with the SR case. By contrast, the PDF results of structural response obtained from the PDE method can easily solve the dynamic reliability while avoiding the use of the backward Kolmogorov equation, thus providing an alternative method for the dynamic reliability analysis of nonlinear stochastic systems. In this study, we will only consider the first-passage failure of the bridge lateral vibration displacement, that is, we will calculate the probability of the bridge lateral vibration displacement not crossing a certain limit value for the first time. The dynamic reliability based on the first passage failure criterion can be expressed as(34)Rt=Pqτ∈Ωs,0≤τ≤t,Ωs denotes the safe domain. If the double symmetry thresholds Ωs⟶qτ≤qlim are considered, equation (34) can be rewritten as(35)Rt=Pqτ≤qlim,0≤τ≤t.Ifqτ∉Ωs, then the structure performance is a failure, and the probability carried by the corresponding event will no longer return to the security domain. This process is equivalent to adding an absorbing boundary condition:(36)pqq,t=0，qτ∉Ωs,0≤τ≤t.On the basis of the boundary conditions of equation (36) and the initial condition equation (21), equation (20) will be resolved by using the finite difference method. Then, a new PDF marked as pq∗q,t with the absorbing boundary condition can be obtained, whereby the dynamic reliability will be subsequently obtained as(37)Rt=∫−∞∞pq∗q,tdq.By using the aforementioned method, the dynamic reliability analysis of footbridge lateral vibration is successfully completed. Figure9 illustrates the results of dynamic reliability under different thresholds (qlim=0.003,0.0035,0.0045) when Np=160. It is found that when qlim=0.003, the value decreases rapidly with time in the initial stage but tends to be stable. In this case, the stable value at the ending time is approximately 0.678, that is, the probability that the lateral vibration does not exceed 0.003 for the first time is 0.678. On the basis of the comparative results of different threshold conditions, the dynamic reliability value will be significantly improved if the passage threshold is increased; such findings are consistent with the expected results. The dynamic reliability of footbridges, especially for those with increasing pedestrian loads and possible damages [61], is a critical index for bridge safety assessment because it can provide a probabilistic sense of evaluation on the normal use or ultimate performance of bridges. Apart from the case of nonlinear stochastic vibration of footbridges, the calculation of dynamic reliability has always been challenging for designers or researchers. However, this study successfully demonstrates that such problems can be solved by using an effective and relatively simple method.Figure 9
Dynamic reliability of bridge lateral vibration.


### 4.3. Random Stability

Random stability is also analyzed on the basis of the obtained probability results. Figure10 shows the probability density flow of lateral displacement under the case of Np=180. In comparison with the case of a smaller number Np=120 (Figure 5(b)), the probability density flow of lateral displacement is no longer stable. In the latter parts of the timeline, the probability density flow begins disorder, and a large amplitude vibration will occur for the first time with a high probability. This phenomenon can also be observed from the joint PDF of displacement and velocity. Figures 11(a) and 11(b) present the joint PDFs of the cases with Np=170 and Np=180, respectively. Figure 11(a) shows that a peak exists around the zero-value region when Np=170, indicating that the probability of small vibration remains high. When the number of people increases to Np=180 (Figure 11(b)), the shape and the number of peaks of joint PDF change. The shape of the single peak in the zero-value region disappears and shifts around, and the overall shape becomes an irregular volcano-shaped peak around the nonzero region, indicating that the possibility of large vibration is relatively high. This result is consistent with the on-site observed result (the number of pedestrians that triggered the divergency is approximately 170) and the results from other researchers [4, 5], which verifies the effectiveness of the proposed method.Figure 10
Contour of the PDF surface during the evolution process (Np=180).Figure 11
Joint PDF of displacement and velocity. (a)Np=170; (b) Np=180.
(a)
(b)


## 4.1. Influence of Randomness

In comparison with the single random (SR) case that considers the randomness of excitation only, the composite random (CR) case (considers the randomness of the structural parameters and excitation) has always been a difficult issue. To date, no studies are available on the CR analysis of footbridge lateral vibration. The structural parameters of the footbridge may have evident randomness that also considerably influences the results due to the measurement error. Figures7(a) and 7(b) show that the PDF result of the CR case is flatter than that of the SR case, indicating that the former dispersion is higher than the latter. This finding is consistent with the expectation because adding the randomness of structural parameters to the system will likely lead to a high discreteness of the response. Given that only three structural parameters are considered random variables in this study, it is believed that the results will be more different from those of the SR case if more structural parameters are involved. Furthermore, to analyze the influence brought by the random disturbance intensity of excitation, the parameter bs in equation (13) is slightly changed to be 0.35. From Figure 7(c), when bs is reduced (i.e., the random disturbance intensity is increased), the PDF around the zero-value region grows larger, which means that increasing the random disturbance intensity is helpful to keep the small vibration.Figure 7
Influence of randomness underNp=120. (a) Comparative PDF results of the CR and SR cases. (b) Comparative CDF results of the CR and SR cases. (c) Comparative PDF results of different bs (under CR case).
(a)
(b)
(c)


## 4.2. Serviceability and Dynamic Reliability

In this study, the RMS acceleration is obtained in a convenient manner of using all the temporal and spatial probability distributions instead of using the power spectrum method. By replacing the displacementqt with acceleration q¨t in equation (2), the probability distributions of acceleration pq¨q¨,t can be obtained by using the same procedure as that for pqq,t. Then, RMS acceleration can be developed by using the following expressions:(33)Γq¨t=Eq¨t2=∫−∞∞q¨2pq¨q¨,tdq¨,σq¨=1T∫0TΓq¨tdt=1T∫0TEq¨t2dt,where Γq¨t is the original moment of second order.In comparison with the definition of the vertical comfort limit on the footbridge, limited references are available for the lateral direction. Bachmann et al. [57] suggested that the limit of lateral peak acceleration should be no more than 0.2 m/s2, and Dallard et al. [58] thought that the limit set to 0.2–0.4 m/s2 was more appropriate. In terms of code, BS 5400 [59] offered that the maximum lateral acceleration should be less than 0.25 m/s2, and EN 1990 [60] suggested that lateral peak acceleration should be less than 0.14f or 0.15 m/s2. In the present study, a strict limit is adopted, wherein the RMS acceleration may not exceed the value of 0.12fs.The time length is still set to 30 s. Figure8(a) shows the time-history results of Γq¨t1/2 when the number of pedestrians Np are 90, 110, 130, and 150. Figure 8(b) presents the results of σq¨ when the number of pedestrians ranges from 90 to 150. As shown in Figure 8(a), Γq¨t1/2 has a low value, and the corresponding RMS acceleration σq¨ is considerably lower than the limit value of comfort (Figure 8(b)) when the number of pedestrians is small. However, σq¨ rapidly increases when the number of pedestrians increases. When the number of people is up to Np=150, σq¨ exceeds the comfort limit, indicating that the pedestrians will start to feel annoyed due to the bridge vibration. If the number of pedestrians continues to increase to a certain critical value, then a large instability vibration of the bridge may occur. At this time, the pedestrian’s feelings may be beyond the concept of comfort but probably panic. Hence, the analysis will no longer belong to the scope of serviceability but the random stability analysis of bridge vibration.Figure 8
Γ
q
¨
t
1
/
2 and σq¨ under different numbers of pedestrians. (a) Γq¨t1/2vs.T under different Np; (b) σq¨ vs. Np.
(a)
(b)As previously mentioned, the conventional analytical methods are mainly based on the Kolmogorov backward equations, which are difficult to apply due to their strict assumptions and inability to deal with the SR case. By contrast, the PDF results of structural response obtained from the PDE method can easily solve the dynamic reliability while avoiding the use of the backward Kolmogorov equation, thus providing an alternative method for the dynamic reliability analysis of nonlinear stochastic systems. In this study, we will only consider the first-passage failure of the bridge lateral vibration displacement, that is, we will calculate the probability of the bridge lateral vibration displacement not crossing a certain limit value for the first time. The dynamic reliability based on the first passage failure criterion can be expressed as(34)Rt=Pqτ∈Ωs,0≤τ≤t,Ωs denotes the safe domain. If the double symmetry thresholds Ωs⟶qτ≤qlim are considered, equation (34) can be rewritten as(35)Rt=Pqτ≤qlim,0≤τ≤t.Ifqτ∉Ωs, then the structure performance is a failure, and the probability carried by the corresponding event will no longer return to the security domain. This process is equivalent to adding an absorbing boundary condition:(36)pqq,t=0，qτ∉Ωs,0≤τ≤t.On the basis of the boundary conditions of equation (36) and the initial condition equation (21), equation (20) will be resolved by using the finite difference method. Then, a new PDF marked as pq∗q,t with the absorbing boundary condition can be obtained, whereby the dynamic reliability will be subsequently obtained as(37)Rt=∫−∞∞pq∗q,tdq.By using the aforementioned method, the dynamic reliability analysis of footbridge lateral vibration is successfully completed. Figure9 illustrates the results of dynamic reliability under different thresholds (qlim=0.003,0.0035,0.0045) when Np=160. It is found that when qlim=0.003, the value decreases rapidly with time in the initial stage but tends to be stable. In this case, the stable value at the ending time is approximately 0.678, that is, the probability that the lateral vibration does not exceed 0.003 for the first time is 0.678. On the basis of the comparative results of different threshold conditions, the dynamic reliability value will be significantly improved if the passage threshold is increased; such findings are consistent with the expected results. The dynamic reliability of footbridges, especially for those with increasing pedestrian loads and possible damages [61], is a critical index for bridge safety assessment because it can provide a probabilistic sense of evaluation on the normal use or ultimate performance of bridges. Apart from the case of nonlinear stochastic vibration of footbridges, the calculation of dynamic reliability has always been challenging for designers or researchers. However, this study successfully demonstrates that such problems can be solved by using an effective and relatively simple method.Figure 9
Dynamic reliability of bridge lateral vibration.


## 4.3. Random Stability

Random stability is also analyzed on the basis of the obtained probability results. Figure10 shows the probability density flow of lateral displacement under the case of Np=180. In comparison with the case of a smaller number Np=120 (Figure 5(b)), the probability density flow of lateral displacement is no longer stable. In the latter parts of the timeline, the probability density flow begins disorder, and a large amplitude vibration will occur for the first time with a high probability. This phenomenon can also be observed from the joint PDF of displacement and velocity. Figures 11(a) and 11(b) present the joint PDFs of the cases with Np=170 and Np=180, respectively. Figure 11(a) shows that a peak exists around the zero-value region when Np=170, indicating that the probability of small vibration remains high. When the number of people increases to Np=180 (Figure 11(b)), the shape and the number of peaks of joint PDF change. The shape of the single peak in the zero-value region disappears and shifts around, and the overall shape becomes an irregular volcano-shaped peak around the nonzero region, indicating that the possibility of large vibration is relatively high. This result is consistent with the on-site observed result (the number of pedestrians that triggered the divergency is approximately 170) and the results from other researchers [4, 5], which verifies the effectiveness of the proposed method.Figure 10
Contour of the PDF surface during the evolution process (Np=180).Figure 11
Joint PDF of displacement and velocity. (a)Np=170; (b) Np=180.
(a)
(b)


## 5. Conclusions

As pedestrian lateral excitation depends on vibration and involves randomness, the pedestrian-induced lateral vibration of footbridges is essentially a complex problem of a nonlinear stochastic vibration system. To solve this problem effectively, a framework based on the PDE method is proposed in this study. First, the parametric resonance stochastic model and the corresponding vibration equations are developed. Second, to avoid the large amount of computation, the PDE, number theoretical, and finite difference methods with the TVD scheme are used to solve the vibration equations at a cost of small number of samples, and the probability results with rich information are obtained. The proposed method is successfully applied to the London Millennium Bridge and the Passerelle Simone de Beauvoir bridge, and its effectiveness is verified by comparing its results with those of the MCS.The serviceability, dynamic reliability, and random stability analyses are conducted on the London Millennium Bridge, based on which it can be summarized as follows:(1)
The proposed method can easily deal with the CR cases, which makes up the inadequacy that only excitation randomness is considered in the random analysis of footbridges in the past. In a future study, in addition to the randomness of excitation, the randomness of structural parameters also needs to be considered because the comparison between the SR and CR cases suggests that the influence brought by the randomness of structural parameters on the probability results is important.(2)
In this study, the serviceability analysis of the London Millennium Bridge is conducted in a convenient manner rather than the previously used power spectrum method. It is found that there is a nonlinear relationship between the number of pedestrians and the RMS acceleration. When the number of pedestrians is small, the RMS acceleration increases linearly with the number of pedestrians, but it will rapidly increase when the number of pedestrians is relatively large. It is also found the RMS acceleration will not exceed the threshold that causes pedestrian uncomfortable if the number of pedestrians is less than 150.(3)
Instead of using the numerical MC method with enormous computation burden or the conventional analytical method based on the backward Kolmogorov equation, the dynamic reliability of the London Millennium Bridge is successfully analyzed by making advantage of the PDF results of structural response. This demonstrates that the complex dynamic reliability analysis of nonlinear stochastic vibration of the footbridge can be dealt with an effective and relatively simple method.(4)
Based on the random stability analysis, it is found that when the number of pedestrians grows to a certain value, the single peak in the zero-value region disappears, while a volcano-shaped peak around the nonzero region arises, suggesting a large vibration may likely happen. The critical number of pedestrians triggering the divergency of vibration is in a good agreement with the existent reports.


## 6. Discussions

The proposed method is not limited to solving such problems and can be used for other nonlinear stochastic models. Hence, the proposed method is worthy of promotion. Nevertheless, there are also some limitations in this study. For the sake of simplicity, the intersubject variability among crowds is not considered in this paper, and the synchronized pedestrians are assumed to be identical. Moreover, although the proposed method can determine the stability of lateral vibration via the density flow or joint PDF, this approach belongs to the numerical analysis method, and only a rough interval of the critical condition can be obtained with multiple trials. Alternatively, the stochastic average method or equivalent nonlinear approximate method can be applied to obtain an analytical expression of the critical condition triggering the vibration instability. However, this type of method can only consider random excitation by ignoring the randomness in structural parameters. Another strategy called the energy criterion method can identify the stability/instability in the nonlinear stochastic system through the difference between the input and output energies. This method can also be used within the frame of the PDE method to obtain the probability of stability/instability. This related research will be completed in our future work.



---

*Source: 2606395-2019-10-27.xml*