Abstract
This paper investigates the dynamic response of a train-bridge system subjected to earthquakes, and the running safety indices of the train on the bridge under earthquake are studied. Taking a long span cable-stayed bridge across the Huangpu River as an example, a full three-dimensional finite element model of the train-bridge system was established, in which the soil-bridge and rail-train interactions were considered. Parallel computing based on contact balance was utilized to deal with this large-scale numerical simulation problem. The dynamic nonlinear analysis was performed on a Hummingbird supercomputer using the finite element code LS-DYNA 971. The results show that the acceleration responses of the train subjected to an earthquake are much greater than the ones without earthquake input, and the running safety of a moving train is affected by both the earthquake intensity and the running speed of the train. The running safety of the moving train can be evaluated by the threshold curve between earthquake intensity and train speed. The proposed modeling strategies and the simulated results can give a reference prediction of the dynamic behaviour of the train-bridge subjected to an earthquake.
1. Introduction
As bridges are indispensable structures for crossing rivers and bays, more and more bridges have been built throughout the world. In China, there are more than 5800 natural rivers. With the social development and economic promotion, numerous bridges carrying both highways and railways have been constructed in China. There has been a wealth of research focusing on train-bridge interaction problems. Diana et al. proposed an analysis of the dynamic behaviour of railway vehicles running on a bridge, and the train’s safety comfort and bridge structure fatigue resistance have been investigated [1]. A spatial dynamic analysis model for a train-bridge system under random excitation has been established by Xia et al. and the whole histories of the train and the responses of bridges have been calculated. Further, they developed a three-dimensional dynamic interaction model for a linear induction motor train and elevated bridge system, from which the dynamic response of the vehicle and bridge system was obtained [2, 3]. Kim et al. proposed a three-dimensional means of analysis for bridge-vehicle interaction, and the analytical dynamic wheel loads and acceleration responses of heavy vehicles and responses of the bridge have been compared with data from field tests to verify the validity of the proposed procedure [4]. Zhang et al. investigated the interactions between a train and a long span cable-stayed bridge using the finite element method [5]. Youcef et al. proposed a modal superposition approach to derive the dynamic response of coupled vehicle-bridge system [6]. A wheel-rail contact formulation for analyzing the nonlinear train-structure interaction that takes into account the wheel and rail geometry was proposed by Montenegro et al. [7].
According to previous research (Dong et al.), the major factor that poses a threat to the safety of bridges is the seismic load [8]. As seismic resistance is the principal concern when designing a bridge, various references tackled the dynamic response of a bridge during an earthquake [9–13]. However, once an earthquake occurs, it is possible that the seismic resistance is good enough to keep the bridge safe but may not be safe enough for moving trains on the bridge.
The running safety of moving trains over bridges during an earthquake receives increasing attention due to the wide use of bridges as railway supporting structures and the increasing speed of trains. Miyamoto et al. evaluated the running safety of trains using a simple three-dimensional car model as the track was vibrated in lateral and vertical directions, but it was assumed that the train is stationary on the track [14]. The dynamic stability of trains moving over bridges shaken by different earthquakes has been studied by Yang and Wu, and various regions of stability have been established for the train experiencing four ground motions [15]. Han et al. investigated the dynamic response of a cable-stayed bridge subjected to running trains and earthquakes, and the running safety has been evaluated [16]. Tanabe et al. described a simple and efficient mechanical model for interaction between wheel and rail. The railway structure including the track is modeled with various finite elements [17]. Combined transient dynamic response of the train and railway structure is obtained by solving the nonlinear equations of motion using a modal method.
In the last decade, a series of experimental and numerical studies were carried out. Kim and Kawatani investigated the seismic response of steel monorail bridges using three-dimensional dynamic response analysis, and the acceleration response of a monorail train has also been calculated to investigate the effect of structural types of bridges on the train’s dynamic response during earthquakes [18]. Xia et al. presented a dynamic model of a coupled train-bridge system, in which the nonuniform characteristics of the seismic wave input from different foundations have been considered [19]. Critical train speeds have been proposed for running safety on high-speed railway bridges during earthquakes. Tanabe et al. used a finite element model to simulate the dynamic interaction of a high-speed train and railway structure during an earthquake [20]. Numerical examples including the comparison with experiments of an actual Shinkansen car on a shaking table have been demonstrated. Yau and Frýba studied the vibration of a suspension bridge due to moving loads of equidistant, identical forces shaken by vertical support motions caused by an earthquake [21]. On this basis, Frýba and Yau carried out extensive investigations on dynamic responses of bridges subjected to moving loads and earthquake excitation [22, 23]. Konstantakopoulos et al. proposed a mathematical model for a combined cable system of bridges, and the dynamic behaviour of suspended bridges subjected to separately or simultaneously acting moving loads and earthquake’s loading was simulated [24, 25]. Ju investigated the derailment of trains moving on multispan bridges using finite element analyses, and he suggested that large pier stiffness should be used to ensure the safety of moving trains during earthquakes [26]. Du et al. presented a finite element method-based framework for dynamic analysis of coupled bridge-train systems under nonuniform seismic ground motion [27]. The results from the case study demonstrate that the spatial variation of seismic ground motion affects the dynamic responses of the bridge-train system. In spite of these works, the running safety of moving trains over bridges during an earthquake is still not clearly understood, and further studies are needed.
In this paper, a recently built long span cable-stayed bridge over the Huangpu River was used as the case study. A full three-dimensional finite element model of a train-bridge system subjected to an earthquake was established, in which the soil-bridge and rail-train interactions were considered. The dynamic behaviours of train-bridge have been simulated, and the running safety of a moving train on a bridge has been discussed in detail. Due to the large scale of the nonlinear dynamic system, a domain decomposition method called Contact Balance Bisection (CBB) was used in this paper, and the model was solved by the Hummingbird supercomputer adopting LS-DYNA.
2. Model Description
2.1. Model of the Cable-Stayed Bridge with Foundation
The study concerns a road-cum-rail bridge called Second Minpu Bridge with double decks situated over the Huangpu River in the city of Shanghai, China. The bridge is supported by a single-pylon and double plane of stay cables and has an overall length of 436 m with a main span of 251 m. Figure 1 shows the design of the bridge. According to the design drawings of the bridge, a full three-dimensional finite element model is constructed using beam, shell, and solid elements. The FE model contains 1373922 nodes and 1380634 elements. The numerical model of the bridge includes all the details, such as tower, cable, truss, deck, rail, and foundation. Figure 2 shows the full three-dimensional FE model of the cable-stayed bridge system with foundation. The modeling method of major components can be summarized as follows.
(a) Cable-stayed bridge and foundation
(b) Main tower
(c) Main piles
(d) Auxiliary piles
(e) Main span of main truss
(f) Side span of main truss
The tower is the key structure on the cable-stayed bridge and mainly consists of a column and up and down cross beams. The piers as well as the tower are hollow concrete structures and modeled by solid elements with eight or six nodes. The FE model of tower and piers is shown in Figure 2. The Colorado cap material model is adopted for the concrete of the tower and piers, and the parameters are listed in Table 1.
The model of the main truss and bridge deck is built using shell elements with three or four nodes. The elastic material is assigned with an elastic modulus; Poisson’s ratio and density are 210 GPa, 0.3 and 7850 kg/m3, respectively. The track including crossties and rail lay on a layer of concrete connected with the deck by weld stud connectors.
The cables are made of high-strength galvanized steel. The equivalent elastic modulus method is the most frequently used method to model the cable. However, this method is suitable for cables with small deformation. For long cables of this bridge, multisegment cable elements are applied to model the cables. Depending on their length, each cable is divided into 30 to 50 elements. The elastic modulus of the cable is 196 GPa. Poisson’s ratio and density are 0.3 and 7850 kg/m3, respectively.
The support piles of the main tower and piers are steel tube piles and concrete piles, respectively. All the piles are built by resultant beam elements with nodes coinciding with the corresponding nodes of pile caps and soil elements. In order to eliminate waves propagating outward from the boundary of foundation, the range of the foundation should be large enough. Hence, the size of the soil foundation is 730 m × 260 m × 90 m. According to the field survey data, a total number of 9 layers are adopted to discretize the soil stratum. The boundary at the bottom of the foundation is assumed to be horizontal and the vertical directions are assumed to be fixed. Viscoelastic artificial boundaries are applied along the four lateral boundaries of the mesh and the boundary conditions are those described in previous research [28–30]. The Drucker-Prager model is adopted for the foundation and the parameters are listed in Table 2.
2.2. Model of the Train
The train is composed of four carriages and each carriage consists of a car body, two bogies, four wheel sets, and four spring-dashpot systems, as shown in Figure 3. The dynamic analysis model of a single carriage is shown in Figure 4, and the main dimensions of the train are indicated in Table 3. The train is modeled by a multibody rigid-flexible model, and the following assumptions are used in the modeling process:(1)The carriages, bogies, and wheel sets are regarded as rigid bodies and the elastic deformations of these components are neglected. The accuracy and efficiency of this assumption have previously been validated by Xia and Zhang [31].(2)Five degrees of freedom are considered in the modeling of the carriage and bogies. They are lateral, vertical, roll, yaw, and pitch displacement. For the wheel set, the pitch displacement is neglected. The four-axle two-bogie vehicle is considered with 27 degrees of freedom in total, as shown in Figure 4.(3)The connections between the carriage, bogies, and wheel sets are represented by viscous dashpots and linear springs.
2.3. Wheel-Rail Interface
The dynamic interaction between the bridge and the train was represented by a wheel-rail contact model. The rails of the track are modeled by two parallel lines of solid elements. In this study, the wheel-rail contact style is treated as being automatic node-to-surface type. The wheel node sets represent the contact patch between wheel and rail, and rail nodes lying in the wheel-rail interface are referred to as master nodes. A penalty function is used to constrain the wheel nodes to slide along the track. A linear force-deflection relationship is assumed in compression, and no tensile force is generated in the process. If there is no penetration, nothing is done. If it does penetrate, a normal contact force can be calculated bywhere denotes normal contact force acting on the wheel and denotes normal contact force acting on the rail; is the stiffness of spring; is the damping coefficient; is the amount of penetration. It is important to evaluate the stiffness coefficient , which is usually determined by the experience of researcher. In this study, the stiffness coefficient is calculated bywhere is the bulk modulus of the element; is the element area; is the element volume; and is the scale factor.
The frictional force (tangential contact force) is defined as where is the nonlinear coefficient of friction; is the viscous damping coefficient; is the relative sliding rate.
The roughness of the rail is a principal source of vibration. According to the American railway standard, the vertical and lateral irregularities of the track can be calculated as follows:where and are the vertical and lateral power spectrum densities; is the spatial frequency; is the coefficient associated with the line grade (from one to six); and are the cutoff frequencies; and is a constant. In this study, class 6 track roughness spectrum parameters are adopted, where and are 0.8245 rad/m and 0.8209 rad/m; is 0.0339 cm2·rad/m; and is 0.25.
2.4. Seismic Loads
A first 20-second Shanghai artificial wave specified by Shanghai’s local code is adopted as excitation, as shown in Figure 5. According to the seismic design code of railway engineering [32], three earthquake levels are specified according to their probabilities of occurrence, including frequent, moderate, and severe. As for railway bridges, it is possible that the bridge itself may remain safe during even a severe earthquake but may not be safe enough for the trains to move over it. In this study, a severe seismic excitation with an occurrence probability of 3% in 50 years is used for seismic designing. The seismic wave is provided by the seismic safety assessment report, and the maximum acceleration in the analysis is normalized as 0.075 g. Regarding the direction of excitation, only the shear wave has been taken into account. The acceleration time history is applied to all the nodes at the bottom of the foundation. All the nodes, at a certain time, have the same acceleration. In other words, the propagating velocity of the seismic waves is infinite.
3. Numerical Methods
In this study, the total degrees of freedom in the full three-dimensional finite element model are greater than 5,000,000. In addition, a high computational cost required for contact detection and computation is very large. It is a challenge for the simulation as well as the data exchange and storage. To solve this large-scale nonlinear problem, explicit time integration and the parallel domain decomposition method are adopted.
3.1. Dynamic Finite Element Method
The finite element method has proved to be an effective method to analyze the dynamic response of a structure during an earthquake. The momentum equations are solved in both the spatial domain and the time domain with the finite element approach. In the spatial domain, the equation of motion for the nonlinear case at time is given bywhere and are the mass and damping matrices, respectively; and are the acceleration and velocity vectors, respectively; is the internal force vector; accounts for the external and body force vector.
In the time domain, (5) can be integrated by the explicit central difference integration rule and is rewritten as follows:
The velocities and displacements are updated in each time step as follows:where .
Compared with implicit schemes, the explicit integration scheme improves the computational efficiency and requires less memory by using a lumped-mass matrix. However, one of the disadvantages of the explicit integration scheme is that the method is conditionally stable. For stability, the time step should be less than the critical time step , which is determined by the characteristic length of the element and its material properties. The critical time step can be represented aswhere is the characteristic length of the element and is the wave speed in the element.
In large-scale problems, the small time step has an adverse effect on the computer runtime. The shortcoming can be overcome by parallel computation, where advantage can be taken of the fact that explicit integration methods are ideally suited for parallel processing of data.
3.2. Parallel Algorithm
The domain decomposition method is the most critical part in the parallel algorithm. The domain decomposition confines the computation and communication to a local space, which is the process of subdividing the spatial domain into subdomains that can be assigned to individual processors for the subsequent parallel solution process. The displacements and forces are computed in each processor, and the information on the boundary of each subdomain is exchanged through a message passing interface (MPI). There are two partitioning methods: the so-called node-cut partitioning and element-cut partitioning, as shown in Figure 6. The elements are assigned uniquely to subdomains in node-cut partitioning. The nodes through which the cut is led are shared by the subdomains. On the other hand, the cut is led across the edges or faces in element-cut partitioning. The nodes are still assigned uniquely to subdomains, and the shared elements are duplicated for each subdomain adjacent to the boundary.
Domain decomposition can be implemented in several ways. The recursive coordinate bisection (RCB) algorithm, which is the default decomposition algorithm of the LS-DYNA program, is generally used. However, it only takes into account the geometrical information of the finite element model. The load balance on each domain is not taken into consideration. The complicated contact taking place during the process requires a large amount of computation, including contact detection and computation. Besides, the highly nonlinear computation with elements moving in and out of contact is also very time-consuming. Taking the structure features of the analysis model and the characteristics of contact computation into account, the decomposition method Contact Balance Bisection (CBB) is put forward in this paper. The properties of this method are as follows:(1)The contact detection and computation are as balanced in each subdomain as possible.(2)The elements and nodes are balanced in each subdomain.(3)The communication between two subdomains is minimized.(4)The convergent velocity of every subdomain is generally increased.
The solution procedure of the explicit finite element modeling parallel algorithm based on CBB is illustrated in Figure 7.
4. Results and Discussion
4.1. Initial Equilibrium Configuration of the Bridge
As a highly statically indeterminate structure, the initial equilibrium state of the cable-stayed bridge is related to the construction stage information. There may be differences between the designed initial configuration and the final configuration. In the study, we supposed that the errors are comparatively small and within the engineering allowance. In order to establish the initial condition for subsequent dynamic simulation, the initial equilibrium state of cable-stayed bridge under gravity was determined by reiterative calculation with finite element method. Two aspects should be taken into consideration: the initial tension forces in cables and the initial configuration of the bridge. The initial tension forces should be identical to the designed value and the initial deformed equilibrium configuration should also be identical to the designed geometric alignment of the bridge deck. Hence, a series of assumed cable tension forces were manipulated as the specified input quantity of the cable element and analysis under a dead load was performed to compare the calculated deck alignment with the designed value. The cable tension forces were then adjusted until the best match was achieved. It is a trial-and-error process which needs a great deal of time and manpower. By means of LS-DYNA second development, a program is developed so that the initial equilibrium states of the bridge are adjusted automatically. Figure 8 shows the cable forces of the simulation results and designed value. It can be seen that the difference between them is less than 2%. Figure 9 shows the comparison of the calculated heights of the upper chords and designed value. It is clear from the figure that they are basically the same and the differences are less than 5%.
4.2. Dynamic Responses of the Bridge
The dynamic responses of the train-bridge system are divided into two main parts: the dynamic response of the bridge and the dynamic response of the train. In this paper, only the uniform seismic excitation is applied. Figure 10 shows the lateral and vertical displacement time histories of the cable-stayed bridge at the middle main-span and side-span during the earthquake when the train is moving on the bridge at a speed of 80 km/h. The maximum lateral displacement of the bridge at the middle main-span and side-span is almost the same, 0.0297 m and 0.0294 m, respectively. According to the lengths of the main-span and side-span, the maximum vertical displacement at the middle main-span is greater than that at the middle side-span. The maximum vertical displacement at the middle main-span is 0.0225 m which is quite small compared with the length (249 m) of the main-span.
Time histories of the lateral and vertical acceleration responses of the bridge are shown in Figure 11. It can be seen that the maximum lateral acceleration of the bridge at the middle main-span and side-span is 3.18 m/s2 and 3.15 m/s2, respectively. Compared with the lateral acceleration responses, the vertical acceleration is quite small with maximums of 0.78 m/s2 and 0.69 m/s2 at the middle main-span and side-span, respectively.
4.3. Dynamic Responses of the Train
Figure 12 shows the lateral and vertical acceleration time histories of the first car body during the earthquake when the speed of the train is 80 km/h. It can be seen that the seismic input has a great effect on both lateral acceleration and vertical acceleration of the first car body. Without the earthquake, the maximum lateral acceleration and vertical acceleration of the first car body are 0.41 c and 0.26 m/s2, respectively. Under the earthquake’s input, the maximum lateral acceleration and vertical acceleration of the first car body increase to 0.83 m/s2 and 0.66 m/s2, respectively.
The maximum responses of lateral and vertical acceleration of the first car body are shown in Figure 13, as a function of train speed. It is clear by reference to the figure that the maximum lateral acceleration and vertical acceleration of the first car body increase with increasing train speed under different earthquake intensities. The influence of the train speed on the lateral acceleration is greater than that on vertical acceleration. When the maximum acceleration is 0.075 g and the train speed increases from 70 km/h to 120 km/h, the lateral acceleration and vertical acceleration increase from 0.78 m/s2 and 0.65 m/s2 to 1.5 m/s2 and 1.05 m/s2 by 92.3% and 61.5%, respectively.
4.4. Running Safety of the Train
The process of a wheel flange climbing over the rail is very complicated. Many researchers have focused on the study of the wheel-rail interaction process. However, there are still many problems that need to be addressed, especially for running safety assessment. Generally, there are two major variables which are the key parameters to govern the derailment during the process, the lateral wheel-rail force , and the vertical wheel-rail force . The relevant factors, derailment factor and offload ratio, are generally used to evaluate the running safety of moving railway vehicles [33, 34]. In this study, we also use these two factors to give evaluation indices for the running safety of the train. The derailment factor is defined as the ratio of the lateral wheel-rail force to the vertical wheel-rail force: . The offload ratio is defined as the ratio of the offload vertical wheel-rail force to the average wheel-rail force wheel set: . The expressions and allowable values of these indices are given in Chinese codes as follows:
The lateral and vertical wheel-rail forces are determined by the penalty method, which is a very sophisticated approach for the simulation of contact phenomena. Figure 14 shows the time histories of vertical and lateral wheel-rail forces relating to the first wheel set without earthquake input, and the interaction wheel-rail forces under seismic input are shown in Figure 15. It can been seen that the influence of the earthquake on the lateral wheel-rail force is more noticeable than that on the vertical one.
The derailment and offload factors of the first wheel set of the first carriage are shown in Figures 16 and 17, when the light rail train travels on the bridge at the speed of 80 km/h, without an earthquake and with uniform seismic input.
(a)
(b)
(a)
(b)
It is clear by reference to Figures 17 and 18 that the running safety indices of the moving train are strongly affected by the earthquake input. The maximum derailment and offload factors without an earthquake are 0.19 and 0.216, respectively. With the earthquake, the maximum running safety indices are both amplified. For the derailment factor, the value increases from 0.19 to 0.391 by 105.8%; the offload factor increases from 0.216 to 0.41 by 89.8%. Moreover, when the cable-stayed bridge is excited by an earthquake with maximum acceleration of 0.075 g, the derailment and offload factor of the light rail train moving at 80 km/h are both within the corresponding safety allowance range.
In order to further investigate the running safety of the moving train subjected to an earthquake, a parameter analysis is performed. Here, train speed and peak ground acceleration are selected as the major factors. With regard to the velocity of the train, a typical operating and limited velocity range varying from 70 km/h to 120 km/h is considered. Figure 18 shows the maximum responses of the derailment factor and offload factor versus train speed, when the train runs across the bridge with a different earthquake input. The horizontal short dash lines represent the allowance in Chinese code. It can be seen that the derailment and offload factors increase with the increase of train speed. When the train speed reaches 90 km/h, the derailment factor during an earthquake with a maximum acceleration of 0.2 g exceeds the allowance value. Nevertheless, when the train speed reaches 80 km/h, the offload factor during the earthquake with a maximum acceleration of 0.2 g exceeds the allowance value. Moreover, it can be found that the train can move safely at the speed of 70 km/h for all maximum acceleration considered.
The maximum responses of the derailment factor and offload factor versus earthquake intensity (maximum acceleration) are shown in Figure 19. The influence of the earthquake on the running safety of the train is similar to that of the train speed. It can be seen that the derailment and offload factors increase with the increasing of the earthquake intensity in most cases. When the maximum acceleration reaches 0.125 g, the derailment factor of the train running at 120 km/h exceeds the allowance value, and when maximum acceleration reaches 0.1 g, the offload factor of the train running at 110 km/h exceeds the allowance value. Moreover, it can be found that the train can move safely during an earthquake with maximum acceleration of 0.075 g for all speeds considered.
It is important to assess the stability of the running train over the bridge with the earthquake. In order to determine the threshold curves for running safety of the train, the train speed and the earthquake intensity are taken as the key parameters as mentioned. The threshold speed can be determined as follows: keeping the train speed or the earthquake intensity as a constant, the derailment and offload factors are calculated by changing one of the key parameters. The corresponding train speed at which any safety factor exceeds the allowances is the threshold speed. Figure 20 shows the train speed—maximum acceleration—running safety indices plot.
The relationship between the earthquake intensity and train speed can be obtained by plotting the simulated results in one coordinated system, as shown in Figure 21. Generally, it can be observed that the greater the earthquake intensity, the lower the critical speed for running safety. In other words, the higher the train speed, the lower the earthquake intensity that the moving train can suffer. In this study, the critical train speed decreases with the increase of the earthquake’s intensity. Moreover, the maximum acceleration of the earthquake is about 0.17 g for the train at a speed of 80 km/h, which is the designed operational speed of the train.
5. Conclusion
The full 3D dynamic analysis of the train-bridge coupled system subjected to the earthquake shaking was proposed using explicit finite element method combined with parallel computing. Taking a long cable-stayed bridge as a case study, the dynamic responses of the bridge and train were computed, and the running safety indices were analyzed. Based on the simulated results from this study, several conclusions can be drawn as follows:(1)The proposed modeling strategies and the simulated results could give a reference prediction of the dynamic behaviour of the train-bridge during an earthquake. The simulated initial equilibrium configuration of the bridge is generally identical with the design plan. The acceleration responses of the train subjected to earthquake are much greater than the ones without earthquake input, and the lateral acceleration and vertical acceleration of the train increase with the increasing train speed.(2)The running safety indices of the train are strongly affected by the earthquake input. The derailment factor and offload factor obviously increase as the earthquake takes place. In this case, the train at a speed of designed operational value (80 km/h) could move across the bridge safely during an earthquake with the maximum acceleration of 0.075 g.(3)The running safety of the moving train is affected by both the earthquake intensity and the running speed of the train. The higher the train speed, the lower the earthquake intensity that the moving train can suffer. The running safety of the moving train could be evaluated by the threshold curve between earthquake intensity and train speed.
Conflict of Interests
All the authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The project was supported by the National 863 High Technology Research and Development Program of China (2012AA01A307) and the National Natural Science Foundation of China (51475287).