Dynamic Response of Bridge Subjected to Eccentrically Moving Flexible Vehicle: A Semianalytical Approach

Abstract

Dynamic response of a single span bridge subjected to moving flexible vehicles has been studied using a semianalytical approach. The eccentricity of vehicle path giving rise to torsional motion of the bridge has been incorporated in the approach. The bridge surface irregularity has been considered as the nonhomogeneous process in spatial domain. A closed form expression has been derived to generate response samples corresponding to each input of roughness profile to form an ensemble. Thereafter, averaging across the ensemble has been carried out at each time step to determine mean and standard deviation of bridge and vehicle response. Further, dynamic amplification factor (DAF) of the bridge response has been obtained for several combinations of bridge-vehicle parameters. The study reveals that structural bending modes of vehicle can significantly reduce dynamic response of the bridge. The eccentricity of vehicle path and flexural/torsional rigidity ratios plays a significant role in dynamic amplification of bridge response.

1. Introduction

The dynamic effect resulting from the passage of vehicles is an important problem generally encountered in the bridge design. The irregularity or unevenness of the bridge pavement surface is the main cause of exciting the vehicle which in turn imposes a time varying load as it travels along the span of the bridge. Starting from the year 1922, various theoretical and experimental studies have been conducted to understand the dynamic behavior of bridge subjected to moving load. A review of literature on the said topic starting from basic formulation with moving mass has been published by Fryba [1] and Yang et al. [2] in their books with detailed discussion on the formulation and its limitations. Earlier researchers have considered vehicle as a moving mass on the bridge, either neglecting their inertia effect or incorporating the same. Although researchers have revealed various dynamic characteristics for practical applications, modern bridges of slender cross-section and larger span do not actually reflect true behavior when moving mass problems have been solved. The deformation of bridge can cause significant change in dynamic forces at the contact point of the vehicle wheel. Realizing these facts, numerous studies have been conducted considering bridge-vehicle a coupled system, with consideration of stiffness and damping parameters of suspension systems. Biggs [3], Fryba [4], and Wen [5] are some of the earlier authors who considered vehicle as single lumped mass system having only bounce motion or a rigid system with bounce and pitch motion. A three-dimensional heave-pitch-roll model has been investigated by Yadav and Upadhyay [6] to find the response of railway tracks on elastic subgrade. Vehicle models with seven and twelve degrees of freedom were developed by Wang et al. [7] according to H20-44 and HS20-44 which are major design vehicles in the American Association of State Highway and Transportation Officials (AASHTO). However, analyses of vehicle motions were confined to rigid modes only. Kou and de Wolf [8], Cheung et al. [9], and Marchesiello et al. [10] have studied the vibrational behavior of multispan continuous bridge by adoption of beam or isotropic plate model. The effect of moving mass or rigid standard vehicle models on the bridge response was examined. Yin et al. [11] investigated lateral vibration of high pier bridge subjected to moving vehicle. In addition to employment of three-dimensional rigid vehicle model, special attention was given to the modeling of tyre with a three-dimensional linear spring and a rectangular contact surface. Other notable works in bridge-vehicle interaction dynamics that include various dynamic properties of rigid vehicle have been reported by Chen and Cai [12], Law and Zhu [13], Zhang et al. [14], and Green and Cebon [15]. Marchesiello et al. [16] studied the interaction of multispan continuous bridges modeled by isotropic plates with multi-degree of freedom vehicles moving at constant speed. Theoretical modes of plate vibration have been obtained with the help of Rayleigh-Ritz approximate method. Prestressed concrete single track railway bridges, consisting of several box girders, were studied by Chu et al. [17] in which they considered each girder to share equal loads and modelled as a beam. Rail irregularities are modeled using power spectral density function. Various works on bridge-vehicle interaction on multispan bridges have been reported by Huang et al. [18], Wang [19], and Ichikawa et al. [20]. Numerical techniques for solving the finite element model with rigid vehicle have been adopted. Other important works which include complexity in bridge model are that of horizontally curved girder [21, 22]. Analytical methods were employed to solve bridge dynamic problem subject to moving mass taking effect of centrifugal force. Effect of vehicle braking has been included in study of vehicle-bridge interaction problem by Kishan and Traill-Nash [23] and Radley and Traill-Nash [24].

Interests in research on bridge-vehicle interaction dynamics are still growing because of several complexities in modeling and uncertainties in excitation. In the last decade, researchers have undertaken more and more complex problems and attempted to solve by newer methodology taking into consideration support flexibility and nonuniform cross-section [25]. A vehicle-track-bridge interaction element considering vehicle’s pitching effect has been developed by Lou [26]. Experimental results and their comparison with finite element model analysis of vehicle-bridge interaction problem have been presented by Brady et al. [27]. Rigid model of vehicle has been considered in the study to determine set of critical velocity associated with peaks of dynamic amplification factor. Multiple resonance response of railway bridge has been investigated by Yau and Yang [28]. A vehicle-bridge interaction process has been simulated with MATLAB Simulink by Harris et al. [29] to find out bridge-friendly damping control strategy with a tractor semitrailer vehicle model. Very recently, some interesting studies of bridge-vehicle or road-vehicle coupled dynamics include stochastic numerics and discrete integration schemes for digital simulation of road-vehicle system by Wedig [30], application of spectral stochastic finite element by Wu and Law [31], and use of spectral matrix operator for direct solution of stochastic coupled differential equations by Kozar and Malic [32]. Those studies have significantly improved the understanding of complex problem in vehicle-bridge interactions. However, vehicle model has been assumed as lumped masses with rigid link connected by suspension elements exhibiting various discrete degrees of freedom. In most of the studies, numerical integration and Monte Carlo simulation techniques have been used. In the modern days, characteristics of vehicles have greatly changed due to incorporation of larger pay load and for increasing demand of traffic. In the past, vehicles have been modeled by a rigid 2D or 3D system having degrees of freedom in bounce, pitch, and roll. However, due to long and slender vehicle plying frequently over the bridges, there is a need to consider flexibility in the vehicle model and to examine the effect of flexible modes of vehicle on the dynamics of bridge. This aspect of bridge-vehicle dynamic interaction is not adequately addressed in the literature. With this in mind, present study has been conducted to find out the response statistics of single span bridge due to movement of flexible vehicle along an eccentric path. By the term “flexible vehicle,” it is understood here that flexural modes have been included in addition to rigid body modes. The bridge has been considered to be under independent transverse bending as well as under torsional excitation arising out of the eccentric path of the vehicle. Nonhomogeneous profile of deck surface has been incorporated in the study by considering a deterministic mean surface superimposed by zero mean random process. Such formulation can take care of defects in surface finishing, construction joints, potholes, bump, approach slab settlement, and so forth. Bridge-vehicle system equations are expressed in state-space form and decoupled at each time step using complex eigenvalue analysis. Closed form expression has been derived for state vectors for each sample input of bridge profile to form the ensemble of response. Finally, averaging across the ensemble has been carried out to determine mean and standard deviation of the response quantities. In present approach numerical integration can be avoided to generate response samples, which can save considerable amount of computational time. The results obtained from present approach have been validated by numerical simulation and available experimental results from the literature. The effect of vehicle velocity on the response and combined effect of several vehicle-bridge parameters on dynamic amplification factor (DAF) have been examined.

2. Mathematical Model

The bridge-vehicle model has been shown in Figure 1. The bridge has been modeled as a uniform beam with simply supported end conditions. The mass, stiffness, and damping are assumed to be uniform along the span of bridge. Due to eccentricity of the vehicle path, the bridge is subjected to flexure as well as torsion. The bridge deck is uneven which has been realized as nonhomogeneous process in spatial domain. This is represented by a function ℎ(𝑥).

Figure 1 

Bridge-vehicle model.

2.1. Equation of Motion of Vehicle

Vehicle body has been idealized as Euler-Bernoulli beam of length 𝑙_𝑣. The behavior of suspension systems consisting of spring and dashpot is assumed as linear. The governing differential equation of motion of the vehicle deflection can be expressed as 𝐸_𝑣𝐼_𝑣𝜕⁴𝑧(𝑢,𝑡)𝜕𝑢⁴+𝐶_𝑣𝜕𝑧(𝑢,𝑡)𝜕𝑡+𝑚_𝑣𝜕²𝑧(𝑢,𝑡)𝜕𝑡²=𝑓_𝑣(𝑢,𝑡) (1) in which 𝑚_𝑣, 𝐸_𝑣𝐼_𝑣, and 𝐶_𝑣 denote the mass per unit length, flexural rigidity, and viscous damping per unit length of the vehicle body and 𝑧(𝑢,𝑡) represents vertical deflection of the vehicle body measured at location 𝑢 from the reference point (taken at the left end of the vehicle) at time instant 𝑡. The impressed vertical force on the vehicle body is given by 𝑓_𝑣()=𝑘𝑢,𝑡_𝑣1𝑧(𝑢,𝑡)−𝑧₁(𝑡)+𝑐_𝑣1𝑧̇(𝑢,𝑡)−𝑧̇₁(𝑡)×𝛿𝑢−𝑢₁+𝑘_𝑣2𝑧(𝑢,𝑡)−𝑧₂(𝑡)+𝑐_𝑣2𝑧̇(𝑢,𝑡)−𝑧̇₂(𝑡)×𝛿𝑢−𝑢₂, (2) where 𝑢₁ and 𝑢₂ represent the location of the attachment point of vehicle suspension from the reference point; 𝑧₁ and 𝑧₂ denote the vertical displacement of front and rear wheel masses respectively. 𝑘_𝑣1 and 𝑘_𝑣2 are the front and rear vehicle suspension stiffness, respectively; 𝑐_𝑣1 and 𝑐_𝑣2 represent damping for vehicle front and rear suspension, respectively, 𝛿 represents Dirac delta function with the property ^∞_−∞𝑓(𝑥)𝛿𝑥−𝑥₁𝑑𝑥=𝑓𝑥₁. (3) The equation of motion for the front unsprung mass is given by 𝑚₁𝑧(̈𝑡)₁+𝑘_𝑡1𝑧₁(𝑡)−𝑦𝑥₁,𝑡−ℎ𝑥₁+𝑘_𝑣1𝑧₁(𝑡)−𝑧𝑢₁,𝑡+𝑐_𝑣1𝑧̇₁(𝑡)−𝑧̇𝑢₁,𝑡+𝑐_𝑡1𝑧̇₁(𝑡)−𝐷𝑦𝑥𝐷𝑡¹,𝑡+ℎ𝑥₁=0. (4) The equation of motion for the rear unsprung mass is given by 𝑚₂𝑧̈₂(𝑡)+𝑘_𝑡2𝑧₂(𝑡)−𝑦𝑥₂,𝑡−ℎ𝑥₂+𝑘_𝑣2𝑧₂(𝑡)−𝑧𝑢₂,𝑡+𝑐_𝑣2𝑧̇₂(𝑡)−𝑧̇𝑢₂,𝑡+𝑐_𝑡2𝑧̇₂(𝑡)−𝐷𝑦𝑥𝐷𝑡²,𝑡+ℎ𝑥₂=0, (5) where 𝑚₁ and 𝑚₂ are front and rear wheel mass, respectively; 𝑘_𝑡1, 𝑘_𝑡2 are front and rear suspension stiffness respectively; and 𝑐_𝑡1, 𝑐_𝑡2 are front and rear suspension damping, respectively. ℎ(𝑥₁) and ℎ(𝑥₂) represent the nonhomogeneous deck profile under the front and rear wheels, respectively. 𝑦(𝑥₁,𝑡) and 𝑦(𝑥₂,𝑡) are bridge displacements under front and rear wheels, respectively, at any instant of time 𝑡. 𝑧(𝑢₁,𝑡) and 𝑧(𝑢₂,𝑡) represent vehicle body deflection at the front and rear wheels position at any instant of time 𝑡. 𝑢₁ and 𝑢₂ are the location of wheel from the end of the vehicle body. Coriolis forces that arise due to rolling of wheel on the deflected profile of the bridge have been considered in the equation of motion using total derivative operator 𝐷/𝐷𝑡 (with 𝐷𝑦/𝐷𝑡=(𝜕𝑦/𝜕𝑥)(𝜕𝑥/𝜕𝑡)+𝜕𝑦/𝜕𝑡) [1, 33].

2.2. Equation of Motion of Bridge

It is assumed that, for symmetrical cross-section (symmetrical about vertical axis), bending and torsion of the bridge would be independent under vertically applied live load. Thus governing differential equation of motion of the bridge in flexure can be expressed as 𝐸_𝑏𝐼_𝑏𝜕⁴𝑦(𝑥,𝑡)𝜕𝑥⁴+𝐶_𝑏𝜕𝑦(𝑥,𝑡)𝜕𝑡+𝑚_𝑏𝜕²𝑦(𝑥,𝑡)𝜕𝑡²=𝑓_𝑏(𝑥,𝑡) (6) in which 𝑚_𝑏, 𝐸_𝑏𝐼_𝑏, and 𝐶_𝑏 represent the mass per unit length, flexural rigidity, and viscous damping per unit length of bridge. The impressed vertical force 𝑓_𝑏(𝑥,𝑡) on the bridge due to vehicle interaction is given by 𝑓_𝑏()𝑘𝑥,𝑡=−_𝑡1𝑧₁(𝑡)−𝑦(𝑥,𝑡)−ℎ(𝑥)+𝑐_𝑡1𝑧̇₁(𝑡)−𝐷𝑦()(𝑥𝐷𝑡𝑥,𝑡+ℎ)×𝛿𝑥−𝑥₁−𝑘_𝑡2𝑧₂(𝑡)−𝑦(𝑥,𝑡)−ℎ(𝑥)+𝑐_𝑡2𝑧̇₂(𝑡)−𝐷𝑦()(𝑥𝛿𝐷𝑡𝑥,𝑡+ℎ)𝑥−𝑥₂−𝑚₁+12𝑚_𝑣𝑙_𝑣𝑔𝛿𝑥−𝑥₁−𝑚₂+12𝑚_𝑣𝑙_𝑣𝑔×𝛿𝑥−𝑥₂+𝑚₁𝐷²𝐷𝑡²𝑦(𝑥,𝑡)+ℎ(𝑥)𝛿𝑥−𝑥₁+𝑚₂𝐷²𝐷𝑡²𝑦(𝑥,𝑡)+ℎ(𝑥)𝛿𝑥−𝑥₂, (7) where 𝑔 is the acceleration due to gravity. The governing differential equation of the bridge in torsion can be written as 𝐺_𝑏𝐽_𝑏𝜕²𝛾(𝑥,𝑡)𝜕𝑥²−𝐶_𝑏𝑇𝜕𝛾(𝑥,𝑡)𝜕𝑡−𝐼_𝑏𝜕²𝛾(𝑥,𝑡)𝜕𝑡²=𝑓_𝑇(𝑥,𝑡), (8) in which 𝐼_𝑏, 𝐺_𝑏𝐽_𝑏, 𝐶_𝑏𝑡, and 𝛾(𝑥,𝑡) represent the mass moment of inertia per unit length, torsional rigidity, distributed viscous damping to rotational motion, and torsional function of bridge, respectively. 𝐽_𝑏 is torsional constant, and 𝐺_𝑏 is the shear modulus of beam material. 𝑓_𝑇(𝑥,𝑡) is the torque produced in the bridge cross-section due to eccentric loading which is given by 𝑓_𝑇()𝑘𝑥,𝑡=−_𝑡1𝑧₁(𝑡)−𝑦(𝑥,𝑡)−ℎ(𝑥)+𝑐_𝑡1𝑧̇₁(𝑡)−𝐷𝑦()(𝑥𝐷𝑡𝑥,𝑡+ℎ)×𝑒_𝑥𝛿𝑥−𝑥₁−𝑘_𝑡2𝑧₂(𝑡)−𝑦(𝑥,𝑡)−ℎ(𝑥)+𝑐_𝑡2𝑧̇₂(𝑡)−𝐷𝑦()(𝑥𝐷𝑡𝑥,𝑡+ℎ)×𝑒_𝑥𝛿𝑥−𝑥₂−𝑚₁+12𝑚_𝑣𝑙_𝑣𝑔𝑒_𝑥𝛿𝑥−𝑥₁−𝑚₂+12𝑚_𝑣𝑙_𝑣𝑔𝑒_𝑥𝛿𝑥−𝑥₂+𝑚₁𝐷²𝐷𝑡²𝑦(𝑥,𝑡)+ℎ(𝑥)𝑒_𝑥𝛿𝑥−𝑥₁+𝑚₂𝐷²𝐷𝑡²𝑦(𝑥,𝑡)+ℎ(𝑥)𝑒_𝑥𝛿𝑥−𝑥_2𝑥𝛿𝑥−𝑥₂. (9) The parameter 𝑒_𝑥 in (9) denotes the eccentricity of vehicle wheels from the centre line of bridge deck.

2.3. Bridge Deck Roughness

In the present study we introduce a roughness, which is nonhomogeneous in space even though vehicle velocity is constant, by adopting the following relation: ℎ(𝑥)=ℎ_𝑚(𝑥)+^𝑁_𝑠=1𝜍_𝑠cos2𝜋Ω_𝑠𝑥+𝜃_𝑠, (10) where ℎ_𝑚(𝑥) is a deterministic mean which represents construction defects, expansion joints, created pot holes, approach slab settlement, expansion joints, development of corrugation, and so forth; 𝜍_𝑠 is the amplitude of cosine wave; and Ω_𝑠 is the spatial frequency (rad/m) within the interval Ω_𝐿Ω_𝑈 in which power spectral density is defined. Ω_𝐿 and Ω_𝑈 are lower and upper cut-off frequencies. The deck roughness is a Gaussian process [34] with a random phase angle 𝜃_𝑠 uniformly distributed from 0 to 2π. 𝑁 is the number of terms used to build up the road surface roughness. The parameters 𝜍_𝑠 and Ω_𝑠 are computed as 𝜍_𝑠=2𝑆Ω_𝑠ΩΔΩ,_𝑠=Ω_𝐿+1𝑠−2ΩΔΩ,ΔΩ=^𝑈−Ω_𝐿𝑁, (11) in which 𝑆(Ω_𝑠) is the power spectral density function (m³/rad) taken from [35] modifying the same with addition of one term in denominator so that the function exists when Ω→0: 𝑆(Ω)=𝑆Ω₀×Ω²₀Ω_𝑠²+Ω²_𝐿. (12) In the above equation, Ω₀=1/2𝜋 rad/m has been taken. The spatial frequency Ω (rad/m) and temporal frequency 𝜔 (rad/s) for the surface profile are related to the vehicle speed 𝑉 (m/s) as 𝜔=Ω𝑉. In the present study, vehicle forward velocity has been assumed as constant.

2.4. Discretization of Flexible Vehicle Equation of Motion

As mentioned earlier vehicle body has been modeled as free-free beam which has two rigid modes and 𝑛_𝑣 number of elastic modes. It can be shown that when the translation of the mass centroid and the rotational motion about the mass centroid are considered, the two motions are orthogonal with respect to each other and with respect to the elastic modes [36]. Thus total displacement of these rigid body degrees of freedom and elastic modes can be described by 𝑧(𝑢,𝑡)=^𝑛𝑣_𝑗=−1𝜙_𝑣𝑗(𝑢)𝜂_𝑗(𝑡), (13) where 𝜙_𝑣𝑗(𝑢) is the vehicle mode shapes; the subscript 𝑣 denotes vehicle; 𝜂_𝑗(𝑡) is the time dependent generalized coordinate; 𝑗 is the mode number; 𝑗=−1,0 are taken to denote rigid body translatory and pitching mode; 𝑗=1,2,3,…,𝑛_𝑣 represent elastic mode sequence of free-free beam; and 𝑛_𝑣 is the number of significant modes of flexible vehicle body. Thus two rigid modes can be written as 𝜙₋₁𝜙=1,₀=𝑢−𝐷₂, (14) where 𝐷₂ is a distance of vehicle centre of gravity from the trailing edge as given in Figure 1.

The elastic bending modes of free-free beam for 𝑗=1,2,3,… are given by [37]: 𝜙_𝑣𝑗=sin𝛼_𝑗𝑢+sinh𝛼_𝑗𝑢+𝛽_𝑗cos𝛼_𝑗𝑢+cosh𝛼_𝑗𝑢𝛽_𝑗=cos𝛼_𝑗𝑙_𝑣+cosh𝛼_𝑗𝑙_𝑣sin𝛼_𝑗𝑙_𝑣+sinh𝛼_𝑗𝑙_𝑣, (15) where corresponding nondimensional frequency parameters 𝛼_𝑗𝑙_𝑣 can be related to circular natural frequency as 𝜔_𝑣𝑗=𝛼²_𝑗𝐸_𝑣𝐼_𝑣𝑚_𝑣𝑙⁴_𝑣. (16) Substituting (13) in (1) and multiplying both sides of the equation by 𝜙_𝑣𝑘(𝑢) and then integrating with respect to 𝑢 from 0 to 𝑙_𝑣 along with orthogonality conditions, the equation of motion can be discretized as 𝜂̈_𝑗(𝑡)+2𝜉_𝑣𝑗𝜔_𝑣𝑗𝜂̇_𝑗(𝑡)+𝜔²_𝑣𝑗𝜂_𝑗(𝑡)=𝑄_𝑣𝑗(𝑡), (17) where 𝑗=−1,0,1,2,….

Generalized force 𝑄_𝑣𝑗(𝑡) in the 𝑗th mode acting on the vehicle is given as 𝑄_𝑣𝑗(𝑡)=1𝑀_𝑣𝑗^𝑙𝑉₀𝑓_𝑣(𝑢,𝑡)𝜙_𝑗(𝑢)𝑑𝑢, (18) in which generalized mass 𝑀_𝑣𝑗 in the 𝑗th mode is given by 𝑀_𝑣𝑗=^𝑙𝑣₀𝑚_𝑣𝜙²_𝑣𝑗(𝑢)𝑑𝑢. (19) Making use of (2) and (13) in (17) and integrating the expression using the property of Dirac delta function, one has the following expression for generalized force: 𝑄_𝑣𝑗(𝑡)=1𝑀_𝑣𝑗𝑘_𝑣1𝑧₁(𝑡)−^𝑛𝑣_𝑗=−1𝜙_𝑗𝑢₁𝜂(𝑡)𝜙_𝑗𝑢₁+𝑐_𝑣1𝑧̇₁(𝑡)−^𝑛𝑣_𝑗=−1𝜙_𝑗𝑢₁𝜂̇(𝑡)𝜙_𝑗𝑢₁+𝑘_𝑣2𝑧₂(𝑡)−^𝑛𝑣_𝑗=−1𝜙_𝑗𝑢₂𝜂(𝑡)𝜙_𝑗𝑢₂+𝑐_𝑣2𝑧̇₂(𝑡)−^𝑛𝑣_𝑗=−1𝜙_𝑗𝑢₂𝜂̇(𝑡)𝜙_𝑗𝑢₂. (20) It may be mentioned that infinite number of modes are possible in continuous system considered in the present study. However, for practical implementation only first 𝑛_𝑣 modes of vehicle body have been included.

2.5. Discretization of Bridge Equation of Motion

Using first 𝑛_𝑏 number of bridge bending modes, the bridge deflection in flexure be written as 𝑦(𝑥,𝑡)=^𝑛𝑏_𝑘=1𝜙_𝑏𝑘(𝑥)𝑞_𝑘(𝑡), (21) where subscript 𝑏 represents bridge, 𝜙_𝑏𝑘(𝑥) is the flexural mode of the beam for simply supported boundary condition corresponding to natural frequency, and 𝜔_𝑏𝑘 and 𝑞_𝑘(𝑡) are generalized coordinates in 𝑘th mode [37].

Now, substituting (21) in (6) and multiplying both sides of the equation by 𝜙_𝑏𝑗(𝑥) and then integrating with respect to 𝑥 from 0 to 𝐿 with the use of orthogonality conditions, the equation of motion can be discretized in normal coordinates as 𝑞̈_𝑘(𝑡)+2𝜉_𝑏𝑘𝜔_𝑏𝑘𝑞̇_𝑘(𝑡)+𝜔²_𝑏𝑘𝑞_𝑘(𝑡)=𝑄_𝑏𝑘(𝑡), (22) where 𝑘=1,2,3,…,𝑛_𝑏.

The generalized force 𝑄_𝑏𝑘(𝑡) in the 𝑘th mode of bridge in flexure is given as 𝑄_𝑏𝑘(𝑡)=1𝑀_𝑏𝑘^𝐿₀𝑓_𝑏(𝑥,𝑡)𝜙_𝑏𝑘(𝑥)𝑑𝑥, (23) in which generalized mass 𝑀_𝑏𝑘 in the 𝑘th mode is given by 𝑀_𝑏𝑘=^𝐿₀𝑚_𝑏𝜙²_𝑏𝑘(𝑥)𝑑𝑥. (24) The generalized force in the 𝑘th of mode of bridge transverse vibration has been worked out as 𝑄_𝑏𝑘(𝑡)1=−𝑀_𝑏𝑘𝑘_𝑡1𝑧₁(𝑡)−^𝑛𝑏_𝑘=1𝜙_𝑏𝑘𝑥₁𝑞_𝑘(𝑡)−ℎ𝑥₁×𝜙_𝑏𝑘𝑥₁+𝑘_𝑡2𝑧₂(𝑡)−^𝑛𝑏_𝑘=1𝜙_𝑏𝑘𝑥₂𝑞_𝑘(𝑡)−ℎ𝑥₂×𝜙_𝑏𝑘𝑥₂+𝑐_𝑡1𝑧̇₁(𝑡)−𝑉^𝑛𝑏_𝑘=1𝜙^′_𝑏𝑘𝑥₁𝑞_𝑘(𝑡)−𝑉^𝑛𝑏_𝑘=1𝜙_𝑏𝑘𝑥₁𝑞̇_𝑘(𝑡)−𝑉ℎ^′𝑥₁×𝜙_𝑏𝑘𝑥₁+𝑐_𝑡2𝑧̇₂(𝑡)−𝑉^𝑛𝑏_𝑘=1𝜙^′_𝑏𝑘𝑥₂𝑞_𝑘(𝑡)−𝑉^𝑛𝑏_𝑘=1𝜙_𝑏𝑘𝑥₂𝑞̇_𝑘(𝑡)−𝑉ℎ^′𝑥₂×𝜙_𝑏𝑘𝑥₂−𝑚₁+12𝑚_𝑣𝑙_𝑣𝑔𝜙_𝑏𝑘𝑥₁−𝑚₂+12𝑚_𝑣𝑙_𝑣𝑔𝜙_𝑏𝑘𝑥₂×𝑚₁^𝑛𝑏_𝑘=1𝜙_𝑏𝑘𝑥₁𝑞̈_𝑘(𝑡)+2𝑉^𝑛𝑏_𝑘=1𝜙^′_𝑏𝑘𝑥₁𝑞̇_𝑘(𝑡)+𝑉^2𝑛𝑏_𝑘=1𝜙^′′_𝑏𝑘𝑥₁𝑞_𝑘(𝑡)+4𝑉ℎ^′′𝑥₁×𝜙_𝑏𝑘𝑥₁+𝑚₂^𝑛𝑏_𝑘=1𝜙_𝑏𝑘𝑥₂𝑞̈_𝑘(𝑡)+2𝑉^𝑛𝑏_𝑘=1𝜙^′_𝑏𝑘𝑥₂𝑞̇_𝑘(𝑡)+𝑉^2𝑛𝑏_𝑘=1𝜙^′′_𝑏𝑘𝑥₂𝑞_𝑘(𝑡)+4𝑉ℎ^′′𝑥₂×𝜙_𝑏𝑘𝑥₂, (25) in which (·) denotes time derivative. Repeating the similar steps, the discretized bridge equation for torsion in normal coordinates can be expressed as 𝛾̈_𝑙(𝑡)+2𝜉_𝑇𝑙𝜔_𝑇𝑙𝛾̇_𝑙(𝑡)+𝜔²_𝑇𝑙𝛾_𝑙(𝑡)=𝑄_𝑇𝑙(𝑡)𝑙=1,2,3,…,𝑛_𝑇, (26) where 𝑛_𝑇 represents number of bridge torsional modes considered and 𝜔_𝑇𝑙 and 𝜉_𝑇𝑙 are the natural frequency and modal damping coefficient of the 𝑙th mode in torsion, respectively. The generalized torque in the 𝑙th mode is given by 𝑄_𝑇𝑙(𝑡)=1𝑀_𝑇𝑙^𝐿₀𝑓_𝑇(𝑥,𝑡)𝜙_𝑇𝑙(𝑥)𝑑𝑥. (27) The torsional natural frequency 𝜔_𝑇𝑙 and corresponding mode 𝜙_𝑇𝑙 for the given simply supported boundary conditions for no warping restrains have been taken from [37]. The generalized mass moment of inertia 𝑀_𝑇𝑙 in the  𝑙th mode is given by 𝑀_𝑇𝑙=^𝐿₀𝐼_𝑏𝜙²_𝑇𝑙(𝑥)𝑑𝑥. (28) The generalized torque in the 𝑙th mode can be expressed as 𝑄_𝑇𝑙(𝑡)𝑒=−^𝑥𝑀_𝑇𝑙𝑘_𝑡1𝑧₁(𝑡)−^𝑛𝑇_𝑙=1𝜙_𝑇𝑙𝑥₁𝛾_𝑙(𝑡)−ℎ𝑥₁×𝜙_𝑇𝑙𝑥₁+𝑘_𝑡2𝑧₂(𝑡)−^𝑛𝑇_𝑙=1𝜙_𝑇𝑙𝑥₂𝛾_𝑙(𝑡)−ℎ𝑥₂×𝜙_𝑇𝑙𝑥₂+𝑐_𝑡1𝑧̇₁(𝑡)−𝑉^𝑛𝑇_𝑙=1𝜙^′_𝑇𝑙𝑥₁𝛾_𝑙(𝑡)−𝑉^𝑛𝑇_𝑙=1𝜙_𝑇𝑙𝑥₁𝛾̇_𝑙(𝑡)−𝑉ℎ^′𝑥₁×𝜙_𝑇𝑙𝑥₁+𝑐_𝑡2𝑧̇₂(𝑡)−𝑉^𝑛𝑇_𝑙=1𝜙^′_𝑇𝑙𝑥₂𝛾_𝑙(𝑡)−𝑉^𝑛𝑇_𝑙=1𝜙_𝑇𝑙𝑥₂𝛾̇_𝑙(𝑡)−𝑉ℎ^′𝑥₂×𝜙_𝑇𝑙𝑥₂−𝑚₁+12𝑚_𝑣𝑙_𝑣𝑔𝜙_𝑇𝑙𝑥₁−𝑚₂+12𝑚_𝑣𝑙_𝑣×𝑔𝜙_𝑇𝑙𝑥₂×𝑚₁^𝑛𝑇_𝑙=1𝜙_𝑇𝑙𝑥₁𝛾̈_𝑙(𝑡)+2𝑉^𝑛𝑇_𝑙=1𝜙^′_𝑇𝑙𝑥₁𝛾̇_𝑙(𝑡)+𝑉^2𝑛𝑇_𝑙=1𝜙^′′_𝑇𝑙𝑥₁𝛾_𝑙(𝑡)+4𝑉ℎ^′′𝑥₁×𝜙_𝑇𝑙𝑥₁+𝑚₂^𝑛𝑇_𝑙=1𝜙_𝑇𝑙𝑥₂𝛾̈_𝑙(𝑡)+2𝑉^𝑛𝑇_𝑙=1𝜙^′_𝑇𝑘𝑥₂𝛾̇_𝑙(𝑡)+𝑉^2𝑛𝑇_𝑘=1𝜙^′′_𝑇𝑘𝑥₂𝛾_𝑙(𝑡)+4𝑉ℎ^′′𝑥₂×𝜙_𝑇𝑙𝑥₂. (29)

3. Method of Solution

The system of (4), (5), (17), (22), and (26) is coupled second order ordinary differential equations. In general for continuous system like the ones (vehicle and bridge) presented in the paper infinite number of modes exist. However, for practical applications, modes have to be truncated to a finite size. Let 𝑛_𝑣, 𝑛_𝑏, and 𝑛_𝑇 be number of significant modes of vehicle motion, bridge flexural, and torsional vibration, respectively. The number of coupled equations becomes 𝑛=2+𝑛_𝑣+𝑛_𝑏+𝑛_𝑇. The system equations can be expressed in matrix notation as [𝑀]{𝑟̈(𝑡)}+[𝐶]{𝑟̇(𝑡)}+[𝐾]{𝑟(𝑡)}={𝐹(𝑡)}, (30) where {𝑟(𝑡)}={𝜂₁(𝑡), 𝜂₂(𝑡),…,𝜂_𝑛𝑣(𝑡), 𝑧₁(𝑡), 𝑧₂(𝑡),𝑞₁(𝑡),𝑞₂(𝑡),…,𝑞_𝑛𝑏(𝑡),𝛾₁(𝑡), 𝛾₂(𝑡),…,𝛾_𝑛𝑇(𝑡)}^𝑇 is the response vector, {𝐹(𝑡)} is the generalized force vector and, [𝑀], [𝐂], and [𝐾] are system mass, damping, and stiffness matrix, respectively.

The system equation can be cast into a 2𝑛-dimensional first order differential equation of the following form: 𝑝̇(𝑡)+[𝐴]𝑝(𝑡)={𝑃(𝑡)}, (31) where 𝑝(𝑡=𝑟̇(𝑡))𝑟(𝑡),{𝑃(𝑡=[𝑀])}⁻¹{𝐹(𝑡)}{0},[𝐴]=[𝑀]⁻¹[𝐶]−[𝐼][𝑀]⁻¹[𝐾][0]. (32) This form is suitable for bridge-vehicle interaction problems, since suspension damping is not small and diagonalization of damping matrix as in case of Rayleigh’s damping may not be fully convincing. [𝐼] is an identity matrix and 𝟎 a null vector or matrix. Let the eigenvalues of the state matrix 𝐀 be 𝛼₁,𝛼₂,𝛼₃,…,𝛼_2𝑛 and the corresponding eigenvectors {𝐮}¹, {𝐮}², {𝐮}³,…,{𝐮}^2𝑛. The modal matrix is defined as [𝑈]={𝑢}¹{𝑢}²,…,{𝑢}^2𝑛. (33) The eigenvalues are assumed to be distinct. Thus one has [𝑈]⁻¹[𝐴][𝑈]=diagonal𝛼₁,𝛼₂,𝛼₃,…,𝛼_2𝑛. (34) Let us assume 𝑝=[𝑈]{𝑣}. (35) Using linear transformation in (35), the state-space equations in decoupled form can be written as 𝑣̇_𝑗(𝑡)+𝛼_𝑗𝑣_𝑗(𝑡)=𝑅_𝑗(𝑡),𝑗=1,2,3,…,2𝑛, (36) with 𝑅_𝑗=^𝑛_𝑘=1𝑢^′𝑛_{𝑗𝑠𝑠=1}𝑚^′_𝑠𝑘𝐹_𝑘(𝑡), (37)𝑢^′_𝑗𝑠 denotes elements in the inverse of the matrix [𝑈] and 𝑚^′_𝑠𝑘 the elements in the inverse of matrix [𝑀].

The general solution of (36) in Stieltjes integral [38] form may be expressed as 𝑣_𝑗(𝑡)=𝑋_0𝑗exp−𝛼_𝑗𝑡+^∞_−∞𝐻_𝑗(𝜔,𝑡)𝑑𝑆𝐹_𝑘(𝜔), (38) where 𝑋_0𝑗 are constants of integration to be determined from the initial conditions. 𝐻_𝑗(𝜔,𝑡) is the transient frequency response function given by [38] 𝐻_𝑗()=1𝜔,𝑡−𝑖𝜔+𝛼_𝑗exp(−𝑖𝜔𝑡)−exp−𝛼_𝑗𝑡−𝑡₀. (39) It can be seen that as 𝑡₀→−∞, for the positive real part of 𝛼_𝑗, 𝐻_𝑗(𝜔,𝑡) approaches the limiting value (1/−𝑖𝜔+𝛼_𝑗) exp(−𝑖𝜔𝑡). This is the characteristics of stable dynamic system. Using (37) and (39) in (38), the response in generalized normal coordinate may be expressed as 𝑟_𝑚(𝑡)=^2𝑛_𝑗=1𝑢_{𝑚+𝑛,𝑗}𝑋_0𝑗exp−𝛼_𝑗𝑡+^2𝑛_𝑗=1𝑢^𝑛_{𝑚+𝑛,𝑗𝑠=1}𝑢^′𝑛_{𝑗𝑠𝑘=1}𝑚^′_𝑠𝑘^∞_−∞𝐻_𝑗(𝜔,𝑡)𝑑𝑆𝐹_𝑘(𝜔)𝑚=1,2,3,…,𝑛. (40) The first term of (40) represents homogeneous solution of the system equation due to initial condition and the second term of the equation may be written using the limiting value of transient frequency response function 𝑟_𝑃𝑚(𝑡)=^2𝑛_𝑗=1𝑢^𝑛_{𝑚+𝑛,𝑗𝑠=1}𝑢^′_𝑗𝑠×^𝑛_𝑘=1𝑚^′_𝑠𝑘^∞_−∞exp(−𝑖𝜔𝑡)−𝑖𝜔+𝛼_𝑗𝐹_𝑘(𝜔)𝑑𝜔. (41) Using Fourier integral [39] and changing the order of integration, (41) can be written as 𝑟_𝑃𝑚(𝑡)=^2𝑛_𝑗=1𝑢^𝑛_{𝑚+𝑛,𝑗𝑠=1}𝑢^′𝑛_{𝑗𝑠𝑘=1}𝑚^′_𝑠𝑘×^∞_−∞𝐹_𝑘(𝜏1)2𝜋^∞_−∞exp[𝑖𝜔(𝜏−𝑡)]−𝑖𝜔+𝛼_𝑗𝑑𝜔𝑑𝜏. (42) Using Cauchy’s residue theorem [39], the particular solution can be expressed as 𝑟_𝑃𝑚(𝑡)=^2𝑛_𝑗=1𝑢^𝑛_{𝑚+𝑛,𝑗𝑠=1}𝑢^′𝑛_{𝑗𝑠𝑘=1}𝑚^′_𝑠𝑘𝐼_𝑗𝑘, (43) where 𝐼_𝑗𝑘=^𝑡₀exp𝛼_𝑗(𝜏−𝑡)𝐹_𝑘(𝜏)𝑑𝜏, (44) in which 𝑡 is the bridge loading time. 𝐼_𝑗𝑘 can be split up for convenience as 𝐼_𝑗𝑘=𝐼_𝑗𝑤+𝐼_𝑗𝑥+𝐼_𝑗𝑦+𝐼_𝑗𝑧, (45) where 𝑤=1,2,3,…,𝑛_𝑣; 𝑥=𝑛_𝑣+1, 𝑛_𝑣+2; 𝑦=𝑛_𝑣+3, 𝑛_𝑣+4,…,𝑛_𝑣+2+𝑛_𝑏;  𝑧=𝑛_𝑣+𝑛_𝑏+3, 𝑛_𝑣+𝑛_𝑏+4,…,𝑛_𝑣+𝑛_𝑏+𝑛_𝑇+2.

Closed form expressions for the components of the above integral which generates each of the response samples have been developed and given in the appendix. The response samples thus form complete ensemble of the process. Averaging across the ensemble at each time step yields mean 𝜇_𝑌 (𝑡_𝑘) and standard deviation 𝜎_𝑌(𝑡_𝑘) of a response process 𝑌.

4. Results and Discussions

4.1. Model Validation

The results obtained by present analytical method have been first compared with numerically simulated results and published results from the literature. Newmark integration scheme [40] has been adopted for numerical simulation. Time step for integration purpose has been selected as 0.005 sec. Two-axle vehicle model data and bridge of single span with simple supports have been used for comparison purpose. For this purpose, bridge and vehicle parameters have been selected from study conducted by Deng and Cai [41] in which they have also presented experimental data. Following are the inputs taken in both analytical model and numerical integration scheme:

span (𝐿): 16.764 m, mass per unit length (𝑚_𝑏): 1.08×10⁴ kg/m, flexural rigidity (𝐸_𝑏𝐼_𝑏): 2.48×10¹² Nm², torsional rigidity (𝐺_𝑏𝐽_𝑏): 2.48×10¹² Nm², mass moment of inertia (𝐼_𝑏): 1.87×10⁵ kgm²/m; 𝐸_𝑏=32 GPa, deck = 25.12 GPa; vehicle mass (𝑚_𝑣): 2.48×10⁴ kg, mass of front and rear wheel (𝑚_𝑤1,𝑚_𝑤2): 725.4 kg each, mass moment of inertia (𝐼_𝑣): 1.72×10⁵ kgm², front and rear suspension stiffness (𝑘_𝑣1,𝑘_𝑣2): 7.28×10⁵ N/m, 1.97×10⁶ N/m, front and rear suspension damping (𝑐_𝑣1,𝑐_𝑣2): 2189.6 Ns/m, 7181.8 Ns/m, and front and rear tyre stiffness (𝑘_𝑡1,𝑘_𝑡2): 1.97×10⁶ N/m, 4.74×10⁶ N/m. In Figure 2, the analytically obtained response sample has been compared with numerically simulated response sample and experimental results found in [41]. The comparison of response sample exhibits close agreement between analytical and numerically simulated results. The experimental data in the initial portion has shown considerable deviation from theoretical results, which, however, becomes closer in the later part of the response time history. The peak magnitude of the displacement in analytical, numerical, and experimental results agrees well.

Figure 2 

Comparison of analytical, numerical, and experimental results.

4.2. Parametric Study

The following system data have been adopted to generate numerical results and to conduct parametric study: a RC slab-girder bridge of span (𝐿): 20 m; three longitudinal girders along the span and three cross-girders at midspan and at supports are provided in bridge; the lane width: 8.6 m, deck thickness: 200 mm, and concrete characteristic strength 25 N/mm². The cross-section of the bridge is shown in Figure 3. A finite element (FE) model of bridge in SAP2000 commercial software is first developed using the above details of the bridge so as to match the fundamental natural frequency and first modal damping ratio of the simply supported beam model of T-beam bridge. The sectional properties of FE model are then used in the present analytical program of the bridge. The following physical parameters are finally selected for the beam model to represent a T-beam RC concrete bridge:

Figure 3 

Cross-section of T-beam bridge (all dimensions are in meter).

mass (𝑚_𝑏): 11.15×10³ kg/m, flexural rigidity (𝐸_𝑏𝐼_𝑏): 3.7×10¹⁰ N-m², and torsional rigidity (𝐺_𝑏𝐽_𝑏): 1.695×10¹⁰ N-m². The sectional parameters to be given in the beam model of bridge have been found from the FE model of bridge in SAP 2000 commercial software to match the fundamental natural frequency and first modal damping of both the beam and the FE model of bridge.

Vehicle parameters are as follows: a long vehicle carrying heavy load often crossing the bridge has been chosen to illustrate the present approach. The standards of vehicle are different from the live load prescribed by bridge code. In the present study, we use a vehicle type TATA 3516C-EX as shown in Figures 4 and 5. This has been idealized as Euler-Bernoulli beam adopted in present formulation. Following are the important physical parameters pertaining to vehicle: length (𝑙_𝑣): 12 m, flexural rigidity (𝐸_𝑣𝐼_𝑣): 5.3×10⁶ N-m², mass per unit length (𝑚_𝑣): 1500 kg/m, front and rear wheel masses (𝑚_𝑤1,  𝑚_𝑤2): 800 kg each, suspension stiffness front and rear (𝑘_𝑣1,𝑘_𝑣2): 3.6×10⁷ N/m, suspension damping front and rear (𝑐_𝑣1,𝑐_𝑣2): 7.2×10⁴ N-sec/m, front and rear tyre stiffness: (𝑘_𝑡1,𝑘_𝑡2): 0.9×10⁷ N/m, and front and rear tyre damping (𝑘_𝑡1,𝑘_𝑡2): 0.7×10⁴ N-sec/m.

Figure 4 

Longitudinal section of vehicle (all dimensions are in meter).

Figure 5 

Vehicle cross-section (all dimensions are in meter).

4.2.1. Bridge Response Statistics

Bridge mean responses at midspan displacement (μ_𝑑𝑦), velocity (μ_𝑣𝑦), and acceleration (μ_𝑎𝑦) subjected to three vehicle speed have been shown in Figures 6, 7, and 8. The roughness of the bridge deck has been assumed to be of poor category as per ISO classification corresponding to roughness coefficient 𝜍_𝑠=16×10⁻⁶ m²/(m/cycle)  [42]. The mean peak central deflection is found to increase with the vehicle speed. The corresponding standard deviation of each response displacement (𝜎_𝑑𝑧), velocity (𝜎_𝑣𝑧), and acceleration (𝜎_𝑎𝑧) has been shown in Figures 9, 10, and 11. It may be noted that bridge response frequency is modified by change in vehicle velocity and obtained behavior has appropriate physical sense that higher velocity increases the temporal frequency of road excitation. This is obvious from left shift of the peak when velocity increases from lower to higher. The peak magnitude is also higher in case of vehicle travelling with higher velocity. No definite pattern is discernible in the standard deviation. Although, the magnitudes of standard deviation are insignificant for practical purpose, magnitudes of the coefficient of variation of peak displacement, velocity, and accelerations are found to be 0.11, 0.183, and 0.105, respectively. It may be mentioned that dynamic loads do not lead to major bridge damage, except in resonance, but they contribute to continuous degradation of bridge increasing necessity of bridge maintenance.

Figure 6 

Mean displacement of bridge at the midspan.

Figure 7 

Mean velocity of bridge at the midspan.

Figure 8 

Mean acceleration of bridge at the midspan.

Figure 9 

Standard deviation of bridge midspan deflection.

Figure 10 

Standard deviation of bridge midspan velocity.

Figure 11 

Standard deviation of bridge midspan acceleration.

4.2.2. Vehicle Response Statistics

The mean vehicle body responses (at centroid)—displacement (𝜇_𝑑𝑧), velocity (𝜇_𝑣𝑧), and acceleration (𝜇_𝑎𝑧) time history—including residual response have been shown in Figures 12, 13, and 14. The corresponding standard deviation of responses displacement (𝜎_𝑑𝑧), velocity (𝜎_𝑣𝑧), and acceleration (𝜎_𝑎𝑧) is shown in Figures 15, 16, and 17. It has been seen that the peak response is increased with increase in vehicle speed. Standard deviation has no definite trend in the present case for vehicle response.

Figure 12 

Mean displacement of vehicle C.G.

Figure 13 

Mean velocity of vehicle at C.G.

Figure 14 

Mean acceleration of vehicle at C.G.

Figure 15 

Standard deviation of vehicle displacement at C.G.

Figure 16 

Standard deviation of vehicle velocity at C.G.

Figure 17 

Standard deviation of vehicle acceleration at C.G.

4.2.3. Effect of Vehicle Flexibility on Mean and Standard Deviation

The contribution of significant number of structural modes of vehicle to bridge response has been examined in a bar diagram presented in Figures 18 and 19 by comparing maximum mean and standard deviation, respectively. It may be seen that inclusion of bending modes of the vehicle reduces the response magnitude of the bridge. The same figure also reveals that when flexible modes of the vehicle are considered, summation of first five modes for finding the response is adequate. Figures 20 and 21 show the comparison of mean and standard deviation of imposed force on the bridge due to passage of a rigid vehicle with those of flexible vehicle. It reveals that when structural bending modes are considered, the imposed force on the bridge is less. This may be attributed to the fact that part total strain energy has been utilized in bending of elastic vehicle body as compared to rigid beam, which reduces imposed force on the bridge.

Figure 18 

Maximum mean displacement of bridge at midspan, (A) rigid vehicle, (B) flexible vehicle only with first structural mode, (C) flexible vehicle with first three structural modes, and (D) flexible vehicle with first five structural modes.

Figure 19 

Maximum standard deviation of displacement at midspan, (A) rigid vehicle, (B) flexible vehicle only with first structural mode, (C) flexible vehicle with first three structural modes, and (D) flexible vehicle with first five structural modes.

Figure 20 

Mean of imposed force time history of the bridge due to passage of rigid and flexible vehicle.

Figure 21 

Standard deviation of imposed force time history of the bridge due to passage of rigid and flexible vehicle.

The comparison of mean response of bridge at midspan with different values of vehicle flexural rigidity has been presented in Figure 22. Corresponding standard deviation is presented in Figure 23. It may be noted that vehicle with lower flexural rigidity generates less response in the bridge. In generating numerical results under this section vehicle velocity is taken as 60 km/h.

Figure 22 

Effect of flexibility on mean response of bridge, (A) 𝐸𝑣𝐼𝑣=5.3×106 N-m2, (B) 𝐸𝑣𝐼𝑣=6.4×107 N-m2, and (C) 𝐸𝑣𝐼𝑣=8.2×1010 N-m2.

Figure 23 

Effect of flexibility on standard deviation of bridge response, (A) 𝐸𝑣𝐼𝑣=5.3×106 N-m2, (B) 𝐸𝑣𝐼𝑣=6.4×107 N-m2, and (C) 𝐸𝑣𝐼𝑣=8.2×1010 N-m2.

4.2.4. Effect of Eccentricity of Vehicle Path on Mean and Standard Deviation

The effect of eccentricity on the bridge response has been studied by varying the loading position from the centre of the bridge width. In the result shown in Figures 24 and 25, 6% to 8% increment in the bridge dynamic responses has been observed as the eccentricity varies from 0.5 m to 1.5 m. The vehicle stiffness for comparing the results in Figures 24 and 25 is taken uniform (𝐸_𝑣𝐼_𝑣=5.3×10⁶ N·m²). In generating numerical results under this section vehicle velocity is taken as 60 km/h.

Figure 24 

Mean of bridge midspan displacement due to different eccentricity of vehicle path (𝑒𝑥), (A) 𝑒𝑥=0.5 m, (B) 𝑒𝑥=1 m, and (C) 𝑒𝑥=1.5 m.

Figure 25 

Standard deviation of bridge midspan displacement due to different eccentricity of vehicle path (𝑒𝑥), (A) 𝑒𝑥=0.5 m, (B) 𝑒𝑥=1 m, and (C) 𝑒𝑥=1.5 m.

4.3. Dynamic Amplification Factor (DAF)

Considering any response variable 𝑌 as random process, the maximum dynamic response can be written as 𝑌_dynamic=𝜇_𝑌𝑥_𝑘,𝑡_𝑘+𝜎_𝑌𝑥_𝑘,𝑡_𝑘, (46) where 𝑌_dynamic denotes the maximum response due to fluctuating load imposed on bridge due to vibratory motion of the vehicle excited by road unevenness.

Thus, dynamic amplification factor (DAF) in this study is defined as 𝑌DAF=^static+𝑌_dynamic𝑌_static, (47) where 𝑌_static refers to the response of the bridge at the midspan location for adverse position of static wheel loads.

4.3.1. Effect of Bridge Torsional Rigidity on Dynamic Amplification Factor (DAF)

The effect of torsional rigidity of bridge has been studied by obtaining DAF for different ratios of torsional rigidity (𝐺_𝑏𝐽_𝑏) to flexural rigidity (𝐸_𝑏𝐼_𝑏) of bridge with various values of vehicle flexural rigidity (𝐸_𝑣𝐼_𝑣) and presented in Figure 26. The results reveal that dynamic amplification factor decreases by an amount of 23% to 34% when the ratio of torsional rigidity to flexural rigidity of bridge increases from 0.002 to 0.01. This indicates that torsionally stiffer bridge produces less dynamic amplification of static live load response.

Figure 26 

Dynamic amplification factor (DAF) with the ratio of bridge torsional rigidity to flexural rigidity (𝐺𝑏𝐽𝑏/𝐸𝑏𝐼𝑏) for different vehicle flexural rigidity. (A) 𝐸𝑣𝐼𝑣=5.3×106 N-m2, (B) 𝐸𝑣𝐼𝑣=6.4×107 N-m2, and (C) 𝐸𝑣𝐼𝑣=8.2×1010 N-m2.

4.3.2. Effect of Vehicle Mass and Speed on DAF

Since dynamic amplification factor depends on several variables, we choose to represent it by surface plot rather than a two-dimensional plot. The effect of vehicle mass as well as speed on the dynamic amplification factor (DAF) has been shown in Figure 27. Individual effect on earlier classical studies [43] has shown that effect of increasing mass has reducing effect on DAF. However, a surface plot shown in Figure 18 shows that combined effect has an increasing trend. Earlier authors argued that reduction of dynamic amplification factor with increasing vehicle weight is due to increase of static deflection. However, the coupled dynamic interaction may be the cause of increased inertia force in addition to suspension force imposed by moving mass at greater velocity, which increases dynamic component of the response.

Figure 27 

Dynamic amplification factor with change in vehicle speed and vehicle mass.

4.3.3. Effect of Bridge Surface Roughness and Speed

The surface roughness and speed of the vehicle are two most influential factors that can cause increased dynamic amplification factor and rapid degradation of the bridge. Bridge dynamic amplification factor has been found by changing bridge surface roughness coefficient for good condition to poor condition as mentioned in ISO specification [42] with change in vehicle speed. Figure 28 shows that poor condition of road induces more deflection on the bridge when vehicle moves over it and also can be catalyzed by vehicle speed.

Figure 28 

Dynamic amplification factor with vehicle speed and bridge surface roughness.

4.3.4. Effect of Bridge Span and Velocity

The span of the bridge is an important factor which decides the impact factor in most of the bridge design codes. Combined effect of bridge span and speed of the vehicle on DAF is not fully known. Figure 29 shows the DAF with variation of bridge span and velocity. It has been found that when the bridge span increases from 14 m to 30 m dynamic amplification factor decreases by an amount of 12% to 15%, irrespective of increase of vehicle velocity. Although increasing span shows a decreasing trend in DAF similar to the earlier studies, the increment found in the present case is not very significant for the range of span 14–30 m

Figure 29 

Dynamic amplification factor with change in bridge span and vehicle speed.

5. Conclusions

In the present study, an analytical solution for coupled vehicle-bridge interaction problem considering eccentrically moving flexible vehicle as well as random deck surface roughness has been obtained. Response samples are generated from the expressions derived in the present study considering nonhomogeneity of deck roughness. The results obtained from analytical approach have been validated with numerical simulation and experimental data available in the literature. Individual and combined effects of several bridge-vehicle parameters have been considered to find out the response statistics. The increased effect of vehicle speed has been found more significant in changing the frequency of imposed oscillation rather than noticeable increase in response peaks. Combined effect of increasing vehicle weight and speed has been found to increase the dynamic amplification factor. The study reveals that flexibility of long vehicle is an important consideration in obtaining bridge response and, due to change in load carrying vehicle configuration, it is now imperative to address this issue in bridge design codes. Torsion of the bridge is activated by eccentric movement of vehicle, and dynamic amplification factor is largely dependent on the ratio of torsional rigidity to flexural rigidity of bridge cross-section.

Appendix

Expression of the Integral 𝐼_𝑗𝑘 for Generating Response Sample

The vector {𝐹(𝑡)} needed to perform the integration is given below: 𝐹_𝑗(𝑡)=0for𝑗=1,2,3,…,𝑛_𝑣=𝑘_{𝑡𝑗−𝑛𝑣}ℎ𝑥_{𝑗−𝑛𝑣}+𝑐_{𝑡𝑗−𝑛𝑣}ℎ̇𝑥_{𝑗−𝑛𝑣}For𝑗=𝑛_𝑣+1,𝑛_𝑣+2=𝑘_𝑡1ℎ𝑥₁𝜙_{𝑏𝑗−(𝑛𝑣+2)}𝑥₁+𝑘_𝑡2ℎ𝑥₂𝜙_{𝑏𝑗−(𝑛𝑣+2)}𝑥₂+𝑐_𝑡1ℎ̇𝑥₁𝜙_{𝑏𝑗−(𝑛𝑣+2)}𝑥₁+𝑐_𝑡2ℎ̇𝑥₂𝜙_{𝑏𝑗−(𝑛𝑣+2)}𝑥₂For𝑗=𝑛_𝑣+3,𝑛_𝑣+4,…,𝑛_𝑣+2+𝑛_𝑏=𝑒_𝑥𝑘_𝑡1ℎ𝑥₁𝜙_{𝑡𝑗−(𝑛𝑣+2+𝑛𝑏)}𝑥₁+𝑘_𝑡2ℎ𝑥₂𝜙_{𝑡𝑗−(𝑛𝑣+2+𝑛𝑏)}𝑥₂+𝑐_𝑡1ℎ̇𝑥₁𝜙_{𝑡𝑗−(𝑛𝑣+2+𝑛𝑏)}𝑥₁+𝑐_𝑡2ℎ̇^′𝑥₂𝜙_{𝑡𝑗−(𝑛𝑣+2+𝑛𝑏)}𝑥₂For𝑗=𝑛_𝑣+𝑛_𝑏𝑛+3,_𝑣+𝑛_𝑏+4,…,𝑛_𝑣+𝑛_𝑏+𝑛_𝑡+2. (A.1) The mean surface profile has been taken as shallow parabolic, the equation of which with respect to one end of the bridge is given by ℎ_𝑚(𝑥)=4ℎ₀𝐿²𝑥(𝐿−𝑥). (A.2) For one trial, a generation of a set of random phase angles 𝜃_𝑠 (𝑠=1,2,…,𝑁) is employed to express a Gaussian process is ℎ_𝑟(𝑥)=^𝑁_𝑠=1𝐴_𝑠cos2𝜋Ω_𝑠𝑥+𝜃_𝑠. (A.3) Further ℎ^′(𝑥)=𝑑ℎ/𝑑𝑥 and ℎ̇(𝑥)=𝑑ℎ/𝑑𝑡=𝑉(𝑑ℎ/𝑑𝑥) where 𝑉 is the speed of vehicle.

For convenience of computation, 𝐼_𝑗𝑘=𝐼_𝑗𝑤+𝐼_𝑗𝑥+𝐼_𝑗𝑦+𝐼_𝑗𝑧, (A.4) where 𝑤=1,2,3,…,𝑛_𝑣; 𝑥=𝑛_𝑣+1, 𝑛_𝑣+2; 𝑦=𝑛_𝑣+2+1, 𝑛_𝑣+2+2,…,𝑛_𝑣+2+𝑛_𝑏; 𝑧=𝑛_𝑣+2+𝑛_𝑏+1, 𝑛_𝑣+2+𝑛_𝑏+2,…,𝑛_𝑣+2+𝑛_𝑏+𝑛_𝑇.

The first term 𝐼_𝑗𝑤 (𝑤=1,2,…,𝑛_𝑣) of (A.4) can be written as 𝐼_𝑗1=𝐼_𝑗2=⋯=𝐼_𝑗𝑛𝑣=0. (A.5) The second term of (A.4) can be written as 𝐼_𝑗𝑥=𝐼_𝑗𝑥1+𝐼_𝑗𝑥2. (A.6) Splitting the first term of integral (A.6) again into two, we get 𝐼_𝑗𝑥1=𝐼_𝑗𝑥11+𝐼_𝑗𝑥12. (A.7) Introducing 𝐴_𝑠=2𝑆(𝑉Ω_𝑠)⁻² and 𝐵_𝑠=2𝜋{Ω_𝐿+(𝑠−0.5)(Ω_𝑈−Ω_𝐿)/𝑁}, the components can be written as 𝐼_𝑗𝑥11=𝑘_𝑡1^𝑁_𝑠=1𝐴_𝑠exp−𝛼_𝑗𝑡𝛼²_𝑗+𝐵_𝑠𝑉²−𝛼_𝑗cos𝜃_𝑠−𝐵_𝑠𝑉sin𝜃_𝑠+exp𝛼_𝑗𝑡×𝛼_𝑗cos𝐵_𝑠𝑉𝑡+𝜃_𝑠+𝐵_𝑠𝑉cos𝐵_𝑠𝑉𝑡+𝜃_𝑠+4ℎ₀𝑉𝛼³_𝑗𝐿²𝛼_𝑗𝐿𝛼_𝑗𝑡+exp−𝛼_𝑗𝑡−𝑉−1𝛼³_𝑗𝐿²𝛼_𝑗𝑡𝛼_𝑗𝑡−2−2exp−𝛼_𝑗𝑡,𝐼+2_𝑗𝑥12=𝑐_𝑡1^𝑁_𝑠=1𝑉𝐴_𝑠𝐵_𝑠exp−𝛼_𝑗𝑡𝛼²_𝑗+𝐵_𝑠𝑉²−𝛼_𝑗sin𝜃_𝑠−𝐵_𝑠𝑉cos𝜃_𝑠+exp𝛼_𝑗𝑡×𝛼_𝑗sin𝐵_𝑠𝑉𝑡+𝜃_𝑠−𝐵_𝑠𝑉sin𝐵_𝑠𝑉𝑡+𝜃_𝑠−4ℎ₀𝑉exp−𝛼_𝑗𝑡𝛼³_𝑗𝐿²×2𝑉exp−𝛼_𝑗𝑡𝛼_𝑗−1𝑡−1+1𝛼_𝑗𝐿exp𝛼_𝑗𝑡−2. (A.8) Similarly, splitting the second term of integral (A.6) into two, we get 𝐼_𝑗𝑥2=𝐼_𝑗𝑥21+𝐼_𝑗𝑥22𝐼_𝑗𝑥21=𝑘_𝑡2^𝑁_𝑠=1𝐴_𝑠exp−𝛼_𝑗𝑡𝛼²_𝑗+𝐵_𝑠𝑉²×−𝛼_𝑗cos𝐵_𝑠𝑙_𝑣−𝜃_𝑠−𝐵_𝑠𝑉sin𝐵_𝑠𝑙_𝑣−𝜃_𝑠+exp𝛼_𝑗𝑡𝛼_𝑗cos𝐵_𝑠(𝑙−𝑉𝑡)−𝜃_𝑠+𝐵_𝑠𝑉sin𝐵_𝑠𝑙_𝑣−𝑉𝑡+𝜃_𝑠+𝑉²2−exp𝛼_𝑗𝑡𝛼_𝑗𝑡𝛼_𝑗+𝑡−2+24ℎ₀𝛼³_𝑗𝐿²−𝛼²_𝑗𝑙_𝑣𝑙_𝑣+𝐿exp𝛼_𝑗𝑡−1+𝛼_𝑗𝑉2𝑙_𝑣+𝐿exp𝛼_𝑗𝑡𝛼_𝑗,𝐼𝑡−1+1_𝑗𝑥22=𝑐_𝑡2^𝑁_𝑠=1𝐴_𝑠𝐵_𝑠𝑉exp−𝛼_𝑗𝑡𝛼²_𝑗+𝐵_𝑠𝑉²𝛼_𝑗sin𝐵_𝑠𝑙_𝑣−𝜃_𝑠+𝐵_𝑠𝑉cos𝐵_𝑠𝑙_𝑣−𝜃_𝑠−exp𝛼_𝑗𝑡𝛼_𝑗sin𝐵_𝑠(𝑙−𝑉𝑡)−𝜃_𝑠+𝐵_𝑠𝑉cos𝐵_𝑠𝑙_𝑣−𝑉𝑡+𝜃_𝑠+4ℎ₀𝑉exp−𝛼_𝑗𝑡𝛼²_𝑗𝐿²×2𝑉exp𝛼_𝑗𝑡𝛼_𝑗𝑡−1+1−𝛼_𝑗2𝑙_𝑣+𝐿exp𝛼_𝑗𝑡−1. (A.9) The third term of (A.4) is written as 𝐼_𝑗𝑦=^𝑛𝑏_𝑘=1𝐼_{𝑗𝑛𝑣+2+𝑘}. (A.10) The 𝑘th term of the above series is composed of ten different terms as given below 𝐼_{𝑗𝑛𝑣+2+𝑘}=¹⁰_𝑎=1𝐼_{(𝑗𝑛𝑣+2+𝑘)𝑎}𝑘=1,2,3,…,𝑛_𝑏. (A.11) The components of (A.11) are as follows:

Introducing the following parameters: 𝑄_𝑠=𝐵_𝑠𝑉;𝑅_𝑘=𝑘𝜋𝑉𝐿;𝑆_𝑘=𝑘𝜋𝐿;𝑍_𝑠=𝐵_𝑠𝑙_𝑣,𝐼_{(𝑗𝑛𝑣+2+𝑘)1}=12𝑘^𝑁_𝑡1𝑠=1𝐴_𝑠2𝑅_𝑘exp−𝛼_𝑗𝑡×𝛼²_𝑗−𝑄²_𝑠+𝑅_𝑘²cos𝜃_𝑠+2𝛼_𝑗𝑄_𝑠sin𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹−𝛼_𝑗sin𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠+𝑅_𝑘−𝑄_𝑠×cos𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²⁻¹+𝛼_𝑗sin𝑄_𝑠+𝑅_𝑘𝑡+𝜃_𝑠−𝑅_𝑘+𝑄_𝑠×cos𝑄_𝑠+𝑅_𝑘𝑡+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹,𝐼_{(𝑗𝑛𝑣+2+𝑘)2}=12𝑘_𝑡14ℎ₀𝑉exp−𝛼_𝑗𝑡𝐿²𝛼²_𝑗+𝑄²_𝑠³×𝐿𝛼²_𝑗+𝑄²_𝑠exp𝛼_𝑗𝑡sin𝑄_𝑠𝑡×𝑄²_𝑠𝛼_𝑗𝑡+1+𝛼²_𝑗𝛼_𝑗𝑡−1−𝑄_𝑠exp𝛼_𝑗𝑡cos𝑄_𝑠𝑡×𝑄²_𝑠𝑡+𝛼_𝑗𝛼_𝑗𝑡−2−2𝛼_𝑗𝑄_𝑠+𝑉2𝑄³_𝑠+𝑄_𝑠exp𝛼_𝑗𝑡cos𝑄_𝑠𝑡×𝑄⁴_𝑠𝑡²+2𝑄²_𝑠𝛼²_𝑗𝑡²−2𝛼_𝑗𝑡−1+𝛼²_𝑗𝛼²_𝑗𝑡²−4𝛼_𝑗𝛼𝑡+6−exp_𝑗𝑡sin𝑄_𝑠𝑡×𝑄⁴_𝑠𝛼_𝑗𝑡+2+𝛼_𝑗𝑄²_𝑠𝛼²_𝑗𝑡²−3+𝛼³_𝑗𝛼²_𝑗𝑡²−2𝛼_𝑗𝑡+2−6𝛼²_𝑗𝑄_𝑠,𝐼_{(𝑗𝑛𝑣+2+𝑘)3}=12𝑘^𝑁_𝑡2𝑠=1𝐴_𝑠exp−𝛼_𝑗𝑡×𝑅_𝑘+𝑄_𝑠cos𝑆_𝑠−𝑍_𝑠−𝜃_𝑠+𝛼_𝑗sin𝑆_𝑘−𝑍_𝑠−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²⁻¹×𝑅_𝑘−𝑄_𝑠cos𝑆_𝑘−𝑍_𝑠+𝜃_𝑠+𝛼_𝑗sin𝑆_𝑘−𝑍_𝑠+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑠²⁻¹−𝛼_𝑗sin𝑄_𝑠−𝑅_𝑘𝑡+𝑆_𝑘−𝑍_𝑠+𝜃_𝑠−𝑅_𝑘+𝑄_𝑠cos𝑄_𝑠−𝑅_𝑘𝑡+𝑆_𝑘−𝑍_𝑠+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑠²⁻¹+𝛼_𝑗sin𝑆_𝑘+𝑍_𝑠−𝑄_𝑠+𝑅𝑡−𝜃_𝑠+𝑅_𝑘+𝑄_𝑠cos𝑆_𝑘+𝑍_𝑠+𝑄_𝑠+𝑅_𝑘𝑡−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹,𝐼_{(𝑗𝑛𝑣+2+𝑘)4}1=−2𝑘_𝑡24ℎ₀𝑉𝐿²𝛼²_𝑗+𝑉²³×sin𝑉𝑡−𝑙_𝑣𝛼_𝑗𝑙²_𝑣𝛼²_𝑗+𝑉²²−𝑙_𝑣𝛼²_𝑗+𝑉²×−𝛼³_𝑗𝐿−𝛼_𝑗𝐿𝑉²+2𝑉³𝛼_𝑗𝑡+1+2𝛼²_𝑗𝑉𝛼_𝑗𝑉𝑡−1+𝑉2𝛼³_𝑗+𝛼_𝑗𝑡𝛼²_𝑗+𝑉²²+2𝑡𝛼⁴_𝑗−𝑉⁴−6𝛼_𝑗𝑉²−𝐿𝛼²_𝑗+𝑉²×𝛼²_𝑗𝛼_𝑗𝑡−1+𝑉²𝛼_𝑗𝑙𝑡+1+𝑉cos_𝑣−𝑉𝑡𝑉⁴𝑙²_𝑣+𝑙_𝑣𝐿+2𝛼²_𝑗𝑡²−4𝛼_𝑗𝑡−2+𝛼²_𝑗𝑉²2𝑙²_𝑣+2𝑙_𝑣𝐿+4𝛼²_𝑗𝑡²+6+𝛼⁴_𝑗𝑙_𝑣𝑙_𝑣+𝐿−𝑉⁵𝑡2𝑙_𝑣×𝛼+𝐿_𝑗𝑡−1−𝛼³_𝑗𝑉2𝑙_𝑣×𝛼+𝐿_𝑗𝑡−2+𝑙²_𝑣𝑉⁶𝐼_{(𝑗𝑛𝑣+2+𝑘)5}=12𝑐^𝑁_𝑡1𝑠=1𝐴_𝑠𝑄_𝑠−2𝑅_𝑘exp−𝛼_𝑗𝑡×𝛼²_𝑗−𝑄²_𝑠+𝑅_𝑘²sin𝜃_𝑠−2𝛼_𝑗𝑄_𝑠cos𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅²×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹−𝛼_𝑗cos𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠+𝑅_𝑘−𝑄_𝑠×sin𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹+𝛼_𝑗cos𝑄_𝑠+𝑅_𝑘𝑡+𝜃_𝑠+𝑅_𝑘+𝑄_𝑠×sin𝑄_𝑠+𝑅_𝑘𝑡+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹𝐼_{(𝑗𝑛𝑣+2+𝑘)6}=−𝑐_𝑡12ℎ₀𝑉exp−𝛼_𝑗𝑡𝐿²𝛼²_𝑗+𝑉²³×exp𝛼_𝑗𝑡cos𝑄_𝑠𝑡×𝑄²_𝑠2𝛼_𝑗𝑡𝑉−𝛼_𝑗𝐿+2𝑉+𝛼²_𝑗2𝛼_𝑗𝑡𝑉−𝛼_𝑗𝐿−2𝑉−𝑄_𝑠exp𝛼_𝑗𝑡sin𝑄_𝑠𝑡×𝑄²_𝑠(𝐿−𝑉𝑡)+𝛼_𝑗𝛼_𝑗𝐿−2𝛼_𝑗,𝐼𝑡𝑉+4𝑉_{(𝑗𝑛𝑣+2+𝑘)7}=12𝑐^𝑁_𝑡2𝑠=1𝐴_𝑠𝑄_𝑠×exp−𝛼_𝑗𝑡𝑅_𝑘+𝑄_𝑠sin𝑆_𝑘+𝑍_𝑠−𝜃_𝑠−𝛼_𝑗cos𝑆_𝑘+𝑍_𝑠−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²⁻¹×𝑄_𝑠−𝑅_𝑘sin𝑆_𝑘−𝑍_𝑠+𝜃_𝑠+𝛼_𝑗cos𝑆_𝑘−𝑍_𝑠+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅²⁻¹+𝛼_𝑗cos𝑆_𝑘+𝑍_𝑠−𝑄_𝑠+𝑅_𝑘𝑡−𝜃_𝑠−𝑅_𝑘+𝑄_𝑠sin𝑆_𝑘+𝑍_𝑠−𝑄_𝑠+𝑅_𝑘𝑡−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹−𝛼_𝑗cos𝑆_𝑘−𝑍_𝑠+𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠+𝑄_𝑠−𝑅_𝑘sin𝑆_𝑘−𝑍_𝑠+𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²⁻¹,𝐼_{(𝑗𝑛𝑣+2+𝑘)8}=−𝑐_𝑡22ℎ₀𝑉exp−𝛼_𝑗𝑡𝐿²𝛼²_𝑗+𝑉²³×exp𝛼_𝑗𝑡sin𝑄_𝑠𝑡×𝑄²_𝑠2𝛼_𝑗𝑡𝑉−𝛼_𝑗𝐿+2𝑉+𝛼²_𝑗2𝛼_𝑗𝑡𝑉−𝛼_𝑗𝐿−2𝑉−𝑄_𝑠exp𝛼_𝑗𝑡cos𝑄_𝑠𝑡×𝑄²_𝑠(𝐿−𝑉𝑡)+𝛼_𝑗𝛼_𝑗𝐿−2𝛼_𝑗,𝐼𝑡𝑉+4𝑉_{(𝑗𝑛𝑣+2+𝑘)9}=−exp−𝛼_𝑗𝑡𝑔𝐿𝑙_𝑣𝑚_𝑣+2𝑚_𝑤12𝐿²𝛼²_𝑗+𝑛²𝜋²𝑉²×𝑘𝜋𝑉+exp𝛼_𝑗𝑡𝑘𝜋𝑉cos𝑘𝜋𝑉𝑡𝐿+𝐿𝛼_𝑗sin𝑘𝜋𝑉𝑡𝐿−exp−𝛼_𝑗𝑡𝑔𝐿𝑙_𝑣𝑚_𝑣+2𝑚_𝑤22𝐿²𝛼²_𝑗+𝑛²𝜋²𝑉²×𝑘𝜋𝑉cos𝑘𝜋𝐷𝐿−exp𝛼_𝑗𝑡cos𝑘𝜋𝐷−𝑉𝑡𝐿−𝐿𝛼_𝑗sin𝑘𝜋𝐷𝐿−exp𝛼_𝑗𝑡sin𝑘𝜋𝐷−𝑉𝑡𝐿,𝐼_{(𝑗𝑛𝑣+2+𝑘)10}=−2𝐿𝑉²𝑚_𝑤1+𝑚_𝑤2×^𝑁_𝑠=1𝐴_𝑠𝐵²_𝑠𝑘𝜋−𝐵_𝑠𝐿𝑉cos𝑄_𝑠−𝑅_𝑠𝑡+𝜃_𝑠+𝐿𝛼_𝑗sin𝑄_𝑠−𝑅_𝑠𝑡+𝜃_𝑠×𝐵²_𝑠𝐿²𝑉²−2𝑘𝜋𝐵_𝑠𝐿𝑉²+𝑘²𝜋²𝑉²+𝐿²𝛼²_𝑗⁻¹−𝑘𝜋+𝐵_𝑠𝐿𝑉cos𝑄_𝑠+𝑅_𝑠𝑡+𝜃_𝑠+𝐿𝛼_𝑗sin𝑄_𝑠+𝑅_𝑠𝑡+𝜃_𝑠×𝐵²_𝑠𝐿²𝑉²+2𝑘𝜋𝐵_𝑠𝐿𝑉²+𝑘²𝜋²𝑉²+𝐿²𝛼²_𝑗⁻¹−exp−𝛼_𝑗𝑡−𝑘𝜋+𝐵_𝑠𝐿𝑉cos𝜃_𝑠+𝐿𝛼_𝑗sin𝜃_𝑠×𝐵²_𝑠𝐿²𝑉²+2𝑘𝜋𝐵_𝑠𝐿𝑉²+𝑘²𝜋²𝑉²+𝐿²𝛼²_𝑗⁻¹+𝑘𝜋−𝐵_𝑠𝐿𝑉cos𝜃_𝑠−𝐿𝛼_𝑗sin𝜃_𝑠×𝑘𝜋−𝐵_𝑠𝐿²𝑉²+𝐿²𝛼²_𝑗⁻¹−8ℎ₀𝑉²𝑚_𝑤1+𝑚_𝑤2𝑘²𝜋²𝑉²+𝐿²𝛼²_𝑗××𝑘𝜋𝑉exp−𝛼_𝑗𝑡−cos𝑘𝜋𝑉𝑡𝐿+𝐿𝛼_𝑗sin𝑘𝜋𝑉𝑡𝐿. (A.12) The fourth term of (A.4) is written as 𝐼_𝑗𝑧=^𝑛𝑡_𝑘=1𝐼_{𝑗𝑛𝑣+2+𝑛𝑏+𝑘}. (A.13) The 𝑘th term of the above series, can be split up into ten parts as given below: 𝐼_{𝑗𝑛𝑣+2+𝑛𝑏+𝑙}=¹⁰_𝑎=1𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)𝑎}𝑘=1,2,3,…,𝑛_𝑇. (A.14) The components of (A.14) are as follows: 𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)1}1=−2𝑘_𝑡1𝑒_𝑥×^𝑁_𝑠=1𝐴_𝑠exp−𝛼_𝑗𝑡×𝛼_𝑗cos𝜃_𝑠+𝑄_𝑠−𝑅_𝑘sin𝜃_𝑠𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²+𝛼_𝑗cos𝜃_𝑠+𝑄_𝑠+𝑅_𝑘sin𝜃_𝑠𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²+𝛼_𝑗cos𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠+𝑄_𝑠−𝑅_𝑘sin𝑄_𝑠−𝑅_𝑘𝑡−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²⁻¹+𝛼_𝑗cos𝑄_𝑠+𝑅_𝑘𝑡+𝜃_𝑠+𝑄_𝑠+𝑅_𝑘sin𝑄_𝑠+𝑅_𝑘𝑡+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹,𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)2}=4𝑘_𝑡1𝑒_𝑥ℎ₀𝐿𝑉𝛼²_𝑗+𝑅_𝑘²³×𝛼²_𝑗+𝑅_𝑘²exp−𝛼_𝑗𝑡×𝛼²_𝑗+exp−𝛼_𝑗𝑡cos𝑅_𝑘𝑡×𝑅_𝑘²𝛼_𝑗𝑡+1+𝛼²_𝑗𝛼_𝑗𝑡−1+𝑅_𝑘exp𝛼_𝑗𝑡sin𝑅_𝑘𝑡×𝑅_𝑘²𝑡+𝛼_𝑗𝛼_𝑗𝑡−2−𝑅_𝑘²+2𝛼_𝑗exp−𝛼_𝑗𝑡𝛼²_𝑗+3𝑅_𝑘²+𝑅_𝑘sin𝑅_𝑘𝑡6𝛼²_𝑗+𝑡²𝛼²_𝑗+𝑅_𝑘²²+4𝛼_𝑗𝑡𝛼²_𝑗+𝑅_𝑘²−2𝑅_𝑘²+cos𝑅_𝑘𝑡2𝑡𝑅_𝑘⁴−𝛼⁴_𝑗+2𝛼_𝑗𝛼²_𝑗−3𝑅_𝑘²×𝛼_𝑗𝑡²𝛼²_𝑗+𝑅_𝑘²²,𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)3}1=−2𝑘_𝑡2𝑒_𝑥×^𝑁_𝑠=1𝐴_𝑠exp−𝛼_𝑗𝑡×𝛼_𝑗cos𝐵_𝑠𝑙_𝑣−𝑆_𝑘−𝜃_𝑠+𝑅_𝑘−𝑄_𝑠sin𝐵_𝑠𝑙_𝑣−𝑆_𝑘−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅²⁻¹+𝛼_𝑗cos𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠+𝑅_𝑘+𝑄_𝑠sin𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹+𝛼_𝑗cos𝑅_𝑘−𝑄_𝑠𝑡+𝐵_𝑠𝑙_𝑣−𝑆_𝑘−𝜃_𝑠+𝑅_𝑘−𝑄_𝑠sin−𝑄_𝑠+𝑅_𝑘𝑡+𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²⁻¹+𝛼_𝑗cos−𝑅_𝑘+𝑄_𝑠𝑡+𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠−𝑅_𝑘+𝑄_𝑠sin−𝑄_𝑠+𝑅_𝑘𝑡+𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹,𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)4}=4𝑘_𝑡2𝑒_𝑥ℎ₀𝐿²𝑙_𝑣𝑙_𝑣+𝐿exp−𝛼_𝑗𝑡×𝛼_𝑗cos𝑆_𝑘−𝑅_𝑘𝑡−𝑅_𝑘sin𝑆_𝑘−𝑅_𝑘𝑡+𝑅_𝑘sin(𝑆)−𝛼_𝑗cos𝑆_𝑘×𝛼²_𝑗+𝑅_𝑘²⁻¹+2𝐿𝑉exp−𝛼_𝑗𝑡𝛼²_𝑗+𝑅_𝑘²²×𝛼²_𝑗+𝑅_𝑘²cos𝑆_𝑘+exp𝛼_𝑗𝑡cos𝑆_𝑘−𝑅_𝑘𝑡×𝑅²𝛼_𝑗𝑡+1+𝛼²_𝑗𝛼_𝑗𝑡−1+𝑅_𝑘𝑅_𝑘²𝑡+𝛼_𝑗𝛼_𝑗𝑆𝑡−2×sin_𝑘−𝑅_𝑘𝑡−2𝛼_𝑗𝑅_𝑘sin𝑆_𝑘−𝑉²𝛼²_𝑗+𝑅_𝑘²³×𝑅_𝑘sin𝑆_𝑘−𝑅_𝑘𝑡6𝛼²_𝑗+𝑡²𝛼²_𝑗+𝑅_𝑘²²−4𝛼_𝑗𝑡𝛼²_𝑗+𝑅_𝑘²−2𝑅_𝑘²+2exp−𝛼_𝑗𝑡×𝛼³_𝑗−3𝛼_𝑗𝑅_𝑘²cos𝑆_𝑘+𝑅_𝑘𝑅_𝑘²−3𝛼²_𝑗×sin𝑆_𝑘−2𝑡𝑅_𝑘⁴−𝛼⁴_𝑗+2𝛼_𝑗𝛼²_𝑗−3𝑅_𝑘²+𝛼_𝑗𝑡²𝛼²_𝑗+𝑅_𝑘²²cos𝑆_𝑘−𝑅_𝑘𝑡,𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)5}1=−2𝑘_𝑡1𝑒_𝑥×^𝑁_𝑠=1𝐴_𝑠exp−𝛼_𝑗𝑡×𝛼_𝑗sin𝜃_𝑠+𝑄_𝑠−𝑅_𝑘cos𝜃_𝑠𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²+𝛼_𝑗sin𝜃_𝑠−𝑄_𝑠+𝑅_𝑘cos𝜃_𝑠𝛼²_𝑗+𝑄_𝑠+𝑅²+𝛼_𝑗sin𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠+𝑅_𝑘−𝑄_𝑠cos𝑄_𝑠−𝑅_𝑘𝑡+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅²⁻¹+𝛼_𝑗sin𝑄_𝑠+𝑅_𝑘𝑡+𝜃_𝑠−𝑄_𝑠+𝑅sin𝑄_𝑠+𝑅_𝑘𝑡+𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹,𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)6}=4𝑐_𝑡1𝑒_𝑥ℎ₀𝐿𝛼²_𝑗+𝑅_𝑘²³𝑉𝐿²×𝛼²_𝑗+𝑅_𝑘²𝛼_𝑗cos𝑅_𝑘𝑡+𝑅_𝑘sin𝑅_𝑘𝑡−𝛼_𝑗exp−𝛼_𝑗𝑡+2exp−𝛼_𝑗𝑡𝛼²_𝑗+exp𝛼_𝑗𝑡cos𝑅_𝑘𝑡×𝑅_𝑘²𝛼_𝑗𝑡+1+𝛼²_𝑗𝛼_𝑗𝑡−1+𝑅_𝑘exp𝛼_𝑗𝑡sin𝑅_𝑘𝑡×𝑅_𝑘²𝑡+𝛼_𝑗𝛼_𝑗𝑡−2−𝑅_𝑘²,𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)7}=4𝑐_𝑡2𝑒_𝑥ℎ₀𝑉𝛼²_𝑗+𝑅_𝑘²²×−𝑅_𝑘sin𝑆_𝑘−𝑅_𝑘𝑡×𝑙_𝑣𝛼²_𝑗+𝑅_𝑘²+𝐿𝛼²_𝑗+𝑅_𝑘²−2𝑉𝑅_𝑘²𝑡+𝛼_𝑗𝛼_𝑗𝑆𝑡−2+cos_𝑘−𝑅_𝑘𝑡×𝑅_𝑘²𝛼_𝑗𝑙_𝑣+𝐿−2𝑉𝛼_𝑗𝑡+1+𝛼²_𝑗𝛼_𝑗𝑙_𝑣+𝐿+2𝑉1−𝛼_𝑗𝑡+exp−𝛼_𝑗𝑡𝑅_𝑘sin𝑆_𝑘×𝑙_𝑣𝛼²_𝑗+𝑅_𝑘²+𝐿𝛼²_𝑗+𝑅_𝑘²+4𝛼_𝑗𝑉+cos𝑆_𝑘×𝑅_𝑘²2𝑉−𝛼_𝑗𝑙_𝑣+𝐿−𝛼²_𝑗𝛼_𝑗𝑙_𝑣,𝐼+𝐿+2𝑉_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)8}=12𝑐_𝑡2𝑒_𝑥×^𝑁_𝑠=1𝐴_𝑠𝑄_𝑠exp−𝛼_𝑗𝑡×𝛼_𝑗cos𝐵_𝑠𝑙_𝑣−𝑆_𝑘−𝜃_𝑠+𝑄_𝑠−𝑅_𝑘sin𝐵_𝑠𝑙_𝑣−𝑆_𝑘−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²⁻¹+𝛼_𝑗sin𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠+𝑅_𝑘+𝑄_𝑠cos𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹−𝛼_𝑗sin𝑅−𝑄_𝑠𝑡+𝐵_𝑠𝑙_𝑣−𝑆_𝑘−𝜃_𝑠+𝑄_𝑠−𝑅_𝑘cos𝑅_𝑘−𝑄_𝑠𝑡+𝐵_𝑠𝑙_𝑣−𝑆_𝑘−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠−𝑅_𝑘²⁻¹−𝛼_𝑗sin−𝑅+𝑄_𝑠𝑡+𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠+𝑅_𝑘+𝑄_𝑠cos−𝑄_𝑠+𝑅_𝑘𝑡+𝐵_𝑠𝑙_𝑣+𝑆_𝑘−𝜃_𝑠×𝛼²_𝑗+𝑄_𝑠+𝑅_𝑘²⁻¹,𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)9}=−𝑒_𝑥exp−𝛼_𝑗𝑡𝑔𝐿𝑙_𝑣𝑚_𝑣+2𝑚_𝑤12𝐿²𝛼²_𝑗+𝑛²𝜋²𝑉²×exp𝛼_𝑗𝑡𝑘𝜋𝑉sin𝑘𝜋𝑉𝑡𝐿−𝐿𝛼_𝑗+exp𝐿𝛼_𝑗𝛼_𝑗𝑡cos𝑘𝜋𝑉𝑡𝐿−𝑒_𝑥exp−𝛼_𝑗𝑡𝑔𝐿𝑙_𝑣𝑚_𝑣+2𝑚_𝑤22𝐿²𝛼²_𝑗+𝑛²𝜋²𝑉²×𝑘𝜋𝑉sin𝑘𝜋𝐷𝐿−exp𝛼_𝑗𝑡sin𝑘𝜋𝐷−𝑉𝑡𝐿−𝑒_𝑥𝐿𝛼_𝑗cos𝑘𝜋𝐷𝐿−exp𝛼_𝑗𝑡×cos𝑘𝜋𝐷−𝑉𝑡𝐿,𝐼_{(𝑗𝑛𝑣+2+𝑛𝑏+𝑙)10}=−2𝐿𝑉²𝑚_𝑤1+𝑚_𝑤2𝑒_𝑥×^𝑁_𝑠=1𝐴_𝑠𝐵²_𝑠𝑘𝜋−𝐵_𝑠𝐿𝑉cos𝑄_𝑠−𝑅_𝑠𝑡+𝜃_𝑠+𝐿𝛼_𝑗sin𝑄_𝑠−𝑅_𝑠𝑡+𝜃_𝑠×𝐵²_𝑠𝐿²𝑉²−2𝑘𝜋𝐵_𝑠𝐿𝑉²+𝑘²𝜋²𝑉²+𝐿²𝛼²_𝑗⁻¹−𝑘𝜋+𝐵_𝑠𝐿𝑉×cos𝑄_𝑠+𝑅_𝑠𝑡+𝜃_𝑠+𝐿𝛼_𝑗sin𝑄_𝑠+𝑅_𝑠𝑡+𝜃_𝑠×𝐵²_𝑠𝐿²𝑉²+2𝑘𝜋𝐵_𝑠𝐿𝑉²+𝑘²𝜋²𝑉²+𝐿²𝛼²_𝑗⁻¹−exp−𝛼_𝑗𝑡×−𝑘𝜋+𝐵_𝑠𝐿𝑉cos𝜃_𝑠+𝐿𝛼_𝑗sin𝜃_𝑠×𝐵²_𝑠𝐿²𝑉²+2𝑘𝜋𝐵_𝑠𝐿𝑉²+𝑘²𝜋²𝑉²+𝐿²𝛼²_𝑗⁻¹+𝑘𝜋−𝐵_𝑠𝐿𝑉cos𝜃_𝑠−𝐿𝛼_𝑗sin𝜃_𝑠×𝑘𝜋−𝐵_𝑠𝐿²𝑉²+𝐿²𝛼²_𝑗⁻¹−8ℎ₀𝑉²𝑒_𝑥𝑚_𝑤1+𝑚_𝑤2𝑘²𝜋²𝑉²+𝐿²𝛼²_𝑗×𝑘𝜋𝑉exp−𝛼_𝑗𝑡−cos𝑘𝜋𝑉𝑡𝐿+𝐿𝛼_𝑗sin𝑘𝜋𝑉𝑡𝐿. (A.15)

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.