Analysis of Train–Track–Bridge Coupling Vibration Characteristics for Heavy-Haul Railway Based on Virtual Work Principle

Abstract

This paper introduces an innovative model for heavy-haul train–track–bridge interaction, utilizing a coupling matrix representation based on the virtual work principle. This model establishes the relationship between the wheel–rail contact surface and the bridge–rail interface concerning internal forces and geometric constraints. In this coupled system’s motion equation, the degrees of freedom (DOFs) of the wheelsets in a heavy-haul train lacking primary suspension are interdependent. Additionally, the vertical and nodding DOFs of the bogie frame are linked with the rail element. A practical application, a Yellow River Bridge with a heavy-haul railway line, is used to examine the accuracy of the proposed model with regard to discrepancy between the simulated and measured displacement ranging from 1% to 11%. A comprehensive parametric analysis is conducted, exploring the impacts of track irregularities of varying wavelengths, axle load lifting, and the degradation of bridge stiffness and damping on the dynamic responses of the coupled system. The results reveal that the bridge’s dynamic responses are particularly sensitive to track irregularities within the wavelength range of 1 to 20 m, especially those within 1 to 10 m. The vertical displacement of the bridge demonstrates a nearly linear increase with heavier axle loads of the heavy-haul trains and the reduction in bridge stiffness. However, there is no significant rise in vertical acceleration under these conditions.

1. Introduction

With the growing demand for lot cargo transport, the operation of heavier, longer and faster trains must achieve a harmonious balance between carrying capacity and efficiency [1,2]. The planned heavy-haul railway line is disrupted by topographical obstacles such as rugged mountains, deep valleys, rivers, and other existing railway lines. The construction of bridges plays a crucial role in ensuring the uninterrupted operation of the railway line. The heavy and cyclic axle loads from trains exert high-intensity and short-time dynamic impacts on both the track and the bridge. Much attention is paid to the analysis of train–track–bridge coupling vibration characteristics caused by the dynamic impacts of heavy-haul trains [3].

Heavy-haul trains involve the significant amplitude and frequency of loading cycles, resulting in the excessive deformation and degradation of tracks and bridges [4,5,6]. The track system, consisting of rails, sleepers, and ballast, serves as a conduit for transferring the dynamic load of heavy-haul trains to the bridge deck. At the same time, the vibration of bridges and tracks exerts an adverse impact on running trains, thereby compromising the stability and safety of train operation [7]. The dynamic responses of the heavy-haul train, track, and bridge subsystems can be obtained through the establishment of a train–track–bridge coupling system. Therefore, the establishment of a train–track–bridge coupling system for heavy-haul railways holds immense significance in terms of the design, operation, and maintenance of railway bridges.

So far, numerical simulation has evolved from the simplest moving load model to the train–bridge coupling model, excluding the track subsystem, and ultimately to the refined train–track–bridge coupling model [8,9,10,11]. The wheel–rail contacts and the rail–bridge interaction play an essential role in determining the dynamic responses of trains, tracks, and bridges. In the iterative method, the subsystems of train and track are coupled through the balance of forces and the compatibility of displacement acting on the wheel–rail contact point [12]. Wang et al. [13] presented an iterative approach for solving the equations of motion pertaining to the vehicle and the track, respectively. The iteration persists until the discrepancy in the interaction force between the wheel and the rail falls below the specified tolerance. Zhang et al. [14] developed an intersystem iterative model, which solves the vehicle and bridge motions throughout the entire simulation. The updated dynamic responses form a new excitation to interact with each other’s system until the specified error threshold is reached. Another way to solve the dynamic responses in the time domain is to consider the train–track–bridge as a globally interconnected large-scale system for coupling [15,16]. Lou et al. [17] regarded the wheel–rail contact force as the internal force, and utilized the principle of total potential energy to derive the motion equation for the train–track–bridge system. Chen et al. [18] approached the vehicle and substructure as a holistic system, incorporating interconnected matrices of stiffness, damping, and mass through the energy variational principle and wheel–rail contact geometry. Compared with the aforementioned two methods, the iterative process of wheel–rail force may not easily or slowly converge due to the drastic dynamic contact changes between wheelsets and rails for heavy-haul trains. The coupling method employs a step-by-step integration method to obtain the dynamic responses of the train, track and bridge simultaneously, which avoids the numerical diffusion problem in the iterative process.

The majority of the relevant numerical models used in the existing literature are derived and analyzed for high-speed railway bridges and high-speed trains [19,20,21]. The main concern in high-speed railway bridge design is the severe vibration caused by vehicle–bridge resonance [22,23]. However, the vibration source complexity of heavy-haul railways surpasses that of high-speed railways. The vibration of heavy-haul railway bridges is primarily attributed to factors such as high axle loads [24], diverse car body suspension modes [25], serious rail wear [26], and ballast softening [27]. The static axle load of high-speed trains typically does not exceed 170 kN, while for heavy-haul trains it may range from 250 kN to 350 kN. The increase in train loading and train formation is likely to have a significant impact on substructures, resulting in dynamic amplifications, fatigue damage, and ballast settlement. Shi et al. [24] established a type of bridge–embankment transition model to predict the acceleration of the transition zone below a heavy-haul railway line under the existing loads. Wang et al. [26] investigated the competition and constrained relationship between fatigue crack damage and side wear for a heavy-haul railway rail. Feng et al. [28] evaluated the mechanical properties and weld fatigue behavior of orthotropic steel bridge decks in heavy-haul railway bridges. Emrah et al. [29] described a detailed investigation of the dynamic behavior of heavy-haul railway bridges through a moving load and moving mass model under various parameter variables, including bridge span, normalized train length, normalized train mass, bridge deck stiffness and train speed. In summary, the dynamic responses of heavy-haul railways need to be obtained through a comprehensive coupling analysis of train–track–bridge systems. However, there is currently a dearth of research in the field of heavy-haul railway bridges.

This study aims to establish and validate an innovative coupling model for heavy-haul trains, tracks, and bridges. This model accurately simulates the dynamic behavior of complex heavy-haul trains, considering that most heavy-haul trains lack primary suspension, unlike high-speed trains. In these heavy-haul trains, the bogie is connected to the rail element through vertical and nodding motions, necessitating the derivation of the wheelsets’ motion from the bogie frame. Utilizing the principle of virtual work, motion equations for each subsystem are derived. Notably, this approach avoids the use of wheel–rail and bridge–rail interfaces as system boundaries, eliminating the need for iterative solutions. The model’s validity is confirmed through its application to a continuous rigid frame bridge with a heavy-haul railway line. Furthermore, the study includes a parametric analysis, considering track irregularities of varying wavelengths, axle load variations, and the degradation of bridge stiffness and damping, thus comprehensively evaluating the system’s dynamic responses.

2. Modeling of Heavy-Haul Train–Track–Bridge Coupling System

2.1. Models of Vehicle, Rail and Bridge

The dynamic characteristics in the axial and vertical directions are exclusively considered in this study, as they are the primary factors contributing to bridge vibration caused by heavy-haul trains, while disregarding transverse deformation. A typical heavy-haul train–track–bridge coupling model is schematically shown in Figure 1.

This model is constructed based on the following five assumptions:

(1)

Rails and bridges exhibit elastic, homogeneous and isotropic behaviors, which are simulated by plane Euler–Bernoulli beam;

(2)

The wheel–rail contact region is a small elliptical area, and the linear Hertz contact theory is employed to simulate the wheel–rail contact relationship;

(3)

The heavy-haul train is modeled as a four-axle multi-rigid-body system;

(4)

The heavy-haul trains maintain a consistent velocity along the bridge span direction, regardless of the longitudinal connection and vibration between the car bodies;

(5)

Ballast is modeled as a continuous spring and damping system to provide track support, and the left and right rails are merged into a single strand in the model.

2.2. Equation of Motion for the Bridge by Virtual Work Principle

Firstly, it is assumed that the length of the bridge element is denoted by l_b, and the Young’s modulus and density per unit length of the bridge are denoted by E_b and ρ_b, respectively. The symbols k_rb and c_rb denote the spring stiffness and damping coefficients of the railway ballast, respectively. The bridge element is characterized by three degrees of freedom (DOFs) at each node, namely, a longitudinal displacement u_b, a vertical displacement v_b and a rotation θ_b about an axis normal to the plane of paper.

Assuming that the cross-section of the bridge element is variable, the interpolation method can be employed to express the cross-sectional area and moment of inertia between the nodes based on the parameters of the section at both ends of the element. The cross-section area and moment of inertia at the coordinate x are denoted by A(x) and I(x), respectively. The beam element of the linear variable’s cross-section is shown in Figure 2.

In the next section, the equation of motion for the bridge is derived by considering the virtual work contributions of the inertia force, elastic deformation, interaction between the rail and bridge, as well as external forces. The virtual work of a bridge element with variable cross-section can be expressed by:

$$\begin{array}{l}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{E}_{\mathrm{b}}{A}_{\mathrm{b}}(x)u_{\mathrm{b}}^{\prime }\delta u_{\mathrm{b}}^{\prime }\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{E}_{\mathrm{b}}{I}_{\mathrm{b}}(x)v_{\mathrm{b}}^{″}\delta v_{\mathrm{b}}^{″}\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{\rho }_{\mathrm{b}}{A}_{\mathrm{b}}(x){\ddot{u}}_{\mathrm{b}}\delta u_{\mathrm{b}}\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{\rho }_{\mathrm{b}}{A}_{\mathrm{b}}(x){\ddot{v}}_{\mathrm{b}}\delta v_{\mathrm{b}}\mathrm{d}x+\\ {{\displaystyle \int }}_0^{l_{\mathrm{b}}}{c}_{\mathrm{b}}{\dot{u}}_{\mathrm{b}}\delta u_{\mathrm{b}}\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{c}_{\mathrm{b}}{\dot{v}}_{\mathrm{b}}\delta v_{\mathrm{b}}\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{k}_{\mathrm{rb}}\left(v_{\mathrm{b}}-v_{\mathrm{r}}\right)\delta v_{\mathrm{b}}\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{c}_{\mathrm{rb}}\left({\dot{v}}_{\mathrm{b}}-{\dot{v}}_{\mathrm{r}}\right)\delta v_{\mathrm{b}}\mathrm{d}x=\delta {\mathbf{d}}_{\mathrm{b}}{\mathbf{f}}_{\mathrm{b}}^{\mathrm{T}}\end{array}$$

(1)

with

$$\begin{gathered} \hspace{1em}{A}_{\mathrm{b}}(x)=A_i+\frac{A_j-A_i}{l_{\mathrm{b}}}x \\ {I}_{\mathrm{b}}(x)=I_i+\frac{I_j-I_i}{l_{\mathrm{b}}}x \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathbf{d}}_{\mathrm{b}}=\left\{\begin{array}{cccccc}{u}_{\mathrm{b}1}& v_{\mathrm{b}1}& {\theta }_{\mathrm{b}1}& u_{\mathrm{b}2}& v_{\mathrm{b}2}& {\theta }_{\mathrm{b}2}\end{array}\right\} \end{gathered}$$

where ${\mathbf{d}}_{\mathrm{b}}$ is the node motion vector at both ends of the bridge element; ${\mathbf{f}}_{\mathrm{b}}$ is external loading force vector acting on the bridge element; the superscripts “$\cdot$” and “$\cdot \cdot$” are the first and second derivatives with respect to time, and the superscripts “′” and “″” are the first and second derivatives with respect to coordinates.

2.3. Equation of Motion for the Rail by Virtual Work Principle

Firstly, it is assumed that the length of the rail element is denoted by l_r, while the Young’s modulus and density per unit length of the rail are denoted by E_r and ρ_r, respectively. The symbols A_r and I_r denote the cross-section area and moment of inertia of the rail, respectively. The rail element is characterized by three DOFs at each node, namely, a longitudinal displacement u_r, a vertical displacement v_r and a rotation θ_r about an axis normal to the plane of paper.

The geometric compatibility condition of the relative vertical displacement between the wheel–rail contact at the kth wheelset and the tth moment can be expressed by:

$$y_k^s=v_{\mathrm{w}k}(t)-v_{\mathrm{r}}(x_k,t)-I{R}_y(x_k)$$

(2)

where $v_{wk}(t)$ is the absolute vertical displacement of the kth wheelset at the tth moment; $v_{\mathrm{r}}$ is the absolute vertical displacement of the rail at the contact point $x_k$ and the moment t; $I{R}_y$ is the vertical geometric irregularity of the rail at the contact point $x_k$.

Considering that there is no tension in the wheel–rail contact, it is unnecessary to model wheel–rail contact and separation separately. Based on the assumption (2), the linear elastic stiffness coefficient of wheel–rail Hertz contact is expressed by $a_{\mathrm{w}k}{k}_{\mathrm{vr}}$. The contact coefficient $a_{\mathrm{w}k}$ should adhere to the following relationship:

$$\left\{\begin{array}{c}\begin{array}{c}{a}_{\mathrm{w}k}=1\hspace{1em}{y}_k^s>0\\ a_{\mathrm{w}k}=0\hspace{1em}{y}_k^s<0\end{array}\end{array}\right.$$

(3)

where the subscript “wk” denotes the kth wheelset.

In the next section, the equation of motion for the rail is derived by considering the virtual work contributions of the inertia force, the elastic deformation, the interaction between the rail and bridge, the contact force between the rail and wheelset, as well as external forces. The virtual work of the rail element can be expressed by:

$$\begin{array}{l}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{E}_{\mathrm{r}}{A}_{\mathrm{r}}{u}_{\mathrm{r}}^{\prime }\delta u_{\mathrm{r}}^{\prime }\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{E}_{\mathrm{r}}{I}_{\mathrm{r}}{v}_{\mathrm{r}}^{″}\delta v_{\mathrm{r}}^{″}\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{\rho }_{\mathrm{r}}{A}_{\mathrm{r}}{\ddot{u}}_{\mathrm{r}}\delta u_{\mathrm{r}}\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{\rho }_{\mathrm{r}}{A}_{\mathrm{r}}{\ddot{v}}_{\mathrm{r}}\delta v_{\mathrm{r}}\mathrm{d}x+\\ {{\displaystyle \int }}_0^{l_{\mathrm{r}}}{k}_{\mathrm{rb}}\left(v_{\mathrm{r}}-v_{\mathrm{b}}\right)\delta v_{\mathrm{r}}\mathrm{d}x+{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{c}_{\mathrm{rb}}\left({\dot{v}}_{\mathrm{r}}-{\dot{v}}_{\mathrm{b}}\right)\delta v_{\mathrm{r}}\mathrm{d}x+{\displaystyle \sum _{k=0}^{n_{\mathrm{w}}}{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}\left(v_{\mathrm{r}}-v_{\mathrm{w}k}-I{R}_y\right)\delta v_{\mathrm{r}}}=\delta {\mathbf{d}}_{\mathrm{r}}{\mathbf{f}}_{\mathrm{r}}^T\end{array}$$

(4)

with

$${\mathbf{d}}_{\mathrm{r}}=\left\{\begin{array}{cccccc}{u}_{\mathrm{r}1}& v_{\mathrm{r}1}& {\theta }_{\mathrm{r}1}& u_{\mathrm{r}2}& v_{\mathrm{r}2}& {\theta }_{\mathrm{r}2}\end{array}\right\}$$

where ${\mathbf{d}}_{\mathrm{r}}$ is the node motion vector at both ends of the rail element; ${\mathbf{f}}_{\mathrm{r}}$ is the external loading force vector acting on the rail element; $n_{\mathrm{w}}$ is the number of wheelsets on the rail element.

2.4. Equation of Motions for a Heavy-Haul Train by Virtual Work Principle

The heavy-haul train simulation employs a multi-rigid-body dynamic model, which consists of a car body, a front bogie, a rear bogie, a secondary suspension, and four wheelsets. The car body and the bogie frames of the ith heavy-haul railway vehicle possess two DOFs, respectively, namely, vertical motion ($v_{\mathrm{c}}$, $v_{\mathrm{t}1}$ and $v_{\mathrm{t}2}$) and nodding motion (${\theta }_{\mathrm{c}}$, ${\theta }_{\mathrm{t}1}$, and ${\theta }_{\mathrm{t}2}$), where the subscript “c” denotes car body, “t1” denotes front bogie, and “t2” denotes rear bogie. Since there is no primary suspension in the majority of heavy-haul trains, the DOFs of the wheelsets are no longer independent and can be derived from the motion of bogie. The positive direction is defined as downward for vertical motion and clockwise for nodding motion. The wheelset is coupled with the rail element through the Hertz contact spring.

The equation of motion for the car body is derived by considering the virtual work contributions of the inertia force acting on the car body, the interaction between the car body and the bogies, as well as gravity. The virtual work of the car body can be expressed by:

$$\begin{array}{l}{m}_{\mathrm{c}}{\ddot{v}}_{\mathrm{c}}\delta v_{\mathrm{c}}+I_{\mathrm{c}}{\ddot{\theta }}_{\mathrm{c}}\delta {\theta }_{\mathrm{c}}+k_2(v_{\mathrm{c}}-v_{\mathrm{t}1})\delta v_{\mathrm{c}}+k_2(v_{\mathrm{c}}-v_{\mathrm{t}2})\delta v_{\mathrm{c}}+k_2{L}_{\mathrm{c}}(L_{\mathrm{c}}{\theta }_{\mathrm{c}}-v_{\mathrm{t}1})\delta {\theta }_{\mathrm{c}}+\\ k_2{L}_{\mathrm{c}}(L_{\mathrm{c}}{\theta }_{\mathrm{c}}+v_{\mathrm{t}2})\delta {\theta }_{\mathrm{c}}+c_2({\dot{v}}_{\mathrm{c}}-{\dot{v}}_{\mathrm{t}1})\delta v_{\mathrm{c}}+c_2({\dot{v}}_{\mathrm{c}}-{\dot{v}}_{\mathrm{t}2})\delta v_{\mathrm{c}}+c_2{L}_{\mathrm{c}}(L_{\mathrm{c}}{\dot{\theta }}_{\mathrm{c}}-{\dot{v}}_{\mathrm{t}1})\delta {\theta }_{\mathrm{c}}+\\ c_2{L}_{\mathrm{c}}(L_{\mathrm{c}}{\dot{\theta }}_{\mathrm{c}}+{\dot{v}}_{\mathrm{t}2})\delta {\theta }_{\mathrm{c}}=m_{\mathrm{c}}g\delta v_{\mathrm{c}}\end{array}$$

(5)

where $m_{\mathrm{c}}$ and $I_{\mathrm{c}}$ are the mass and rotational inertia of the car body; $k_2$ and $c_2$ are the spring stiffness and damping coefficient of secondary suspension; $L_{\mathrm{c}}$ is half of the longitudinal distance between the center of gravity of the front bogie and the rear bogie.

The equation of motion for the front bogie is derived by considering the virtual work contributions of the inertia force acting on the front bogie, the interaction between the car body and the front bogie, the contact between the wheelsets and the rail, as well as gravity. The virtual work of the front bogie can be expressed by:

$$\begin{array}{l}(m_{\mathrm{t}}+2m_{\mathrm{w}}){\ddot{v}}_{\mathrm{t}1}\delta v_{\mathrm{t}1}+(I_{\mathrm{t}}{\ddot{\theta }}_{\mathrm{t}1}+2m_{\mathrm{w}}{l}_{\mathrm{t}}^2){\theta }_{\mathrm{t}1}\delta {\theta }_{\mathrm{t}1}+k_2(v_{\mathrm{t}1}-v_{\mathrm{c}}-L_{\mathrm{c}}{\theta }_{\mathrm{c}})\delta v_{\mathrm{t}1}+c_2({\dot{v}}_{\mathrm{t}1}-{\dot{v}}_{\mathrm{c}}-L_{\mathrm{c}}{\dot{\theta }}_{\mathrm{c}})\delta v_{\mathrm{t}1}+\\ a_{\mathrm{w}1}{k}_{\mathrm{wr}}(v_{\mathrm{t}1}+l_{\mathrm{t}}{\theta }_{\mathrm{t}1}-v_{\mathrm{r}}-I{R}_{y1})\delta v_{\mathrm{t}1}+a_{\mathrm{w}1}{k}_{\mathrm{wr}}{l}_{\mathrm{t}}(v_{\mathrm{t}1}+l_{\mathrm{t}}{\theta }_{\mathrm{t}1}-v_{\mathrm{r}}-I{R}_{y1})\delta {\theta }_{\mathrm{t}1}+a_{\mathrm{w}2}{k}_{\mathrm{wr}}(v_{\mathrm{t}1}-l_{\mathrm{t}}{\theta }_{\mathrm{t}1}-v_{\mathrm{r}}-I{R}_{y2})\delta v_{\mathrm{t}1}-\\ a_{\mathrm{w}2}{k}_{\mathrm{wr}}{l}_{\mathrm{t}}(v_{\mathrm{t}1}-l_{\mathrm{t}}{\theta }_{\mathrm{t}1}-v_{\mathrm{r}}-I{R}_{y2})\delta {\theta }_{\mathrm{t}1}=(m_{\mathrm{t}}+2m_{\mathrm{w}})g\delta v_{\mathrm{t}1}\end{array}$$

(6)

where $m_{\mathrm{t}}$ and $I_{\mathrm{t}}$ are the mass and rotational inertia of the bogie frame; $m_{\mathrm{w}}$ is the mass of the wheelset; $l_{\mathrm{t}}$ is half of the bogie axle’s base; the subscripts “w1” and “w2” are the front and rear wheelsets of the front bogie.

The equation of motion for the rear bogie is derived by considering the virtual work contributions of the inertia force acting on the rear bogie, the interaction between the car body and the rear bogie, the contact between the wheelsets and the rail, as well as gravity. The virtual work of the rear bogie can be expressed by:

$$\begin{array}{l}(m_{\mathrm{t}}+2m_{\mathrm{w}}){\ddot{v}}_{\mathrm{t}2}\delta v_{\mathrm{t}2}+(I_{\mathrm{t}}{\ddot{\theta }}_{\mathrm{t}2}+2m_{\mathrm{w}}{l}_{\mathrm{t}}^2){\theta }_{\mathrm{t}2}\delta {\theta }_{\mathrm{t}2}+k_2(v_{\mathrm{t}2}-v_{\mathrm{c}}+L_{\mathrm{c}}{\theta }_{\mathrm{c}})\delta v_{\mathrm{t}2}+c_2({\dot{v}}_{\mathrm{t}2}-{\dot{v}}_{\mathrm{c}}+L_{\mathrm{c}}{\dot{\theta }}_{\mathrm{c}})\delta v_{\mathrm{t}2}+\\ a_{\mathrm{w}3}{k}_{\mathrm{wr}}(v_{\mathrm{t}2}+l_{\mathrm{t}}{\theta }_{\mathrm{t}2}-v_{\mathrm{r}}-I{R}_{y3})\delta v_{\mathrm{t}2}+a_{\mathrm{w}3}{k}_{\mathrm{wr}}{l}_{\mathrm{t}}(v_{\mathrm{t}2}+l_{\mathrm{t}}{\theta }_{\mathrm{t}2}-v_{\mathrm{r}}-I{R}_{y3})\delta {\theta }_{\mathrm{t}2}+a_{\mathrm{w}4}{k}_{\mathrm{wr}}(v_{\mathrm{t}2}-l_{\mathrm{t}}{\theta }_{\mathrm{t}2}-v_{\mathrm{r}}-I{R}_{y4})\delta v_{\mathrm{t}2}-\\ a_{\mathrm{w}4}{k}_{\mathrm{wr}}{l}_{\mathrm{t}}(v_{\mathrm{t}2}-l_{\mathrm{t}}{\theta }_{\mathrm{t}2}-v_{\mathrm{r}}-I{R}_{y4})\delta {\theta }_{\mathrm{t}2}=(m_{\mathrm{t}}+2m_{\mathrm{w}})g\delta v_{\mathrm{t}2}\end{array}$$

(7)

where the subscripts “w3” and “w4” are the front and rear wheelsets of the rear bogie.

The equations of motion for wheelsets are not independent. The displacements, velocities, and accelerations of the wheelset are dependent on the corresponding displacements, velocities, and accelerations of the bogie.

$$v_{\mathrm{w}1}=v_{\mathrm{t}1}+l_{\mathrm{t}}{\theta }_{\mathrm{t}1},\hspace{1em}{\dot{v}}_{\mathrm{w}1}={\dot{v}}_{\mathrm{t}1}+l_{\mathrm{t}}{\dot{\theta }}_{\mathrm{t}1},\hspace{1em}{\ddot{v}}_{\mathrm{w}1}={\ddot{v}}_{\mathrm{t}1}+l_{\mathrm{t}}{\ddot{\theta }}_{\mathrm{t}1}$$

(8)

$$v_{\mathrm{w}2}=v_{\mathrm{t}1}-l_{\mathrm{t}}{\theta }_{\mathrm{t}1},\hspace{1em}{\dot{v}}_{\mathrm{w}2}={\dot{v}}_{\mathrm{t}1}-l_{\mathrm{t}}{\dot{\theta }}_{\mathrm{t}1},\hspace{1em}{\ddot{v}}_{\mathrm{w}2}={\ddot{v}}_{\mathrm{t}1}-l_{\mathrm{t}}{\ddot{\theta }}_{\mathrm{t}1}$$

(9)

$$v_{\mathrm{w}3}=v_{\mathrm{t}2}+l_{\mathrm{t}}{\theta }_{\mathrm{t}2},\hspace{1em}{\dot{v}}_{\mathrm{w}3}={\dot{v}}_{\mathrm{t}2}+l_{\mathrm{t}}{\dot{\theta }}_{\mathrm{t}2},\hspace{1em}{\ddot{v}}_{\mathrm{w}3}={\ddot{v}}_{\mathrm{t}2}+l_{\mathrm{t}}{\ddot{\theta }}_{\mathrm{t}2}$$

(10)

$$v_{w4}=v_{\mathrm{t}2}-l_{\mathrm{t}}{\theta }_{\mathrm{t}2},\hspace{1em}{\dot{v}}_{w4}={\dot{v}}_{\mathrm{t}2}-l_{\mathrm{t}}{\dot{\theta }}_{\mathrm{t}2},\hspace{1em}{\ddot{v}}_{\mathrm{w}4}={\ddot{v}}_{\mathrm{t}2}-l_{\mathrm{t}}{\ddot{\theta }}_{\mathrm{t}2}$$

(11)

3. Assembly and Solution of Heavy-Haul Train–Track–Bridge Coupling Matrix

The heavy-haul train, track and bridge constitute an integrated system in which the vibration characteristics vary with the passage of time during train crossings on the railway. Based on the equation of motion derived from the virtual work principle in Section 2 and following the rule of “seat by number” for matrix formation, the dynamic equations of motion for the heavy-haul train–track–bridge coupling system can be expressed in matrix form, that is:

$$\left[\begin{array}{ccc}{\mathbf{M}}_{\mathrm{bb}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{M}}_{\mathrm{rr}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{M}}_{\mathrm{vv}}\end{array}\right]\left\{\begin{array}{c}{\ddot{\mathbf{X}}}_{\mathrm{b}}\\ {\ddot{\mathbf{X}}}_{\mathrm{r}}\\ {\ddot{\mathbf{X}}}_{\mathrm{v}}\end{array}\right\}+\left[\begin{array}{ccc}{\mathbf{C}}_{\mathrm{bb}}& {\mathbf{C}}_{\mathrm{br}}& \mathbf{0}\\ {\mathbf{C}}_{\mathrm{rb}}& {\mathbf{C}}_{\mathrm{rr}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{C}}_{\mathrm{vv}}\end{array}\right]\left\{\begin{array}{c}{\dot{\mathbf{X}}}_{\mathrm{b}}\\ {\dot{\mathbf{X}}}_{\mathrm{r}}\\ {\dot{\mathbf{X}}}_{\mathrm{v}}\end{array}\right\}+\left[\begin{array}{ccc}{\mathbf{K}}_{\mathrm{bb}}& {\mathbf{K}}_{\mathrm{br}}& \mathbf{0}\\ {\mathbf{K}}_{\mathrm{rb}}& {\mathbf{K}}_{\mathrm{rr}}& {\mathbf{K}}_{\mathrm{rv}}\\ \mathbf{0}& {\mathbf{K}}_{\mathrm{vr}}& {\mathbf{K}}_{\mathrm{vv}}\end{array}\right]\left\{\begin{array}{c}{\mathbf{X}}_{\mathrm{b}}\\ {\mathbf{X}}_{\mathrm{r}}\\ {\mathbf{X}}_{\mathrm{v}}\end{array}\right\}=\left\{\begin{array}{c}{\mathbf{F}}_{\mathrm{b}}\\ {\mathbf{F}}_{\mathrm{r}}\\ {\mathbf{F}}_{\mathrm{v}}\end{array}\right\}$$

(12)

where M, C and K are the submatrices of mass, damping and stiffness respectively; $\mathbf{X}$, $\dot{\mathbf{X}}$ and $\ddot{\mathbf{X}}$ are the vectors of displacement, velocity and acceleration, respectively; F is the vector of force; the subscripts “b”, “r” and “v” are the bridge, rail and heavy-haul vehicle, respectively.

3.1. Submatrices of Bridge and Bridge-Rail Coupling

The symbols “${\mathbf{N}}_{u,j}$” and “${\mathbf{N}}_{v,j}$” represent axial and vertical interpolation functions of the bridge element with the same order 1 × 6, respectively. The former is the primary Hermitian interpolation function and the latter is the cubic Hermitian interpolation function [30], namely,

$$v_{w4}=v_{\mathrm{t}2}-l_{\mathrm{t}}{\theta }_{\mathrm{t}2},\hspace{1em}{\dot{v}}_{w4}={\dot{v}}_{\mathrm{t}2}-l_{\mathrm{t}}{\dot{\theta }}_{\mathrm{t}2},\hspace{1em}{\ddot{v}}_{\mathrm{w}4}={\ddot{v}}_{\mathrm{t}2}-l_{\mathrm{t}}{\ddot{\theta }}_{\mathrm{t}2}$$

(13)

$${\mathbf{N}}_{u,j}=\left\{\begin{array}{cccccc}{N}_1& 0& 0& N_2& 0& 0\end{array}\right\}$$

(14)

with

$$\begin{array}{l}{N}_1=1-\left(x/l_j\right),\hspace{1em}{N}_2=x/l_j,\hspace{1em}{N}_3=1-3{\left(x/l_j\right)}^2+2{\left(x/l_j\right)}^3\\ N_4=x{(1-x/l_j)}^2,\hspace{1em}{N}_5=3{\left(x/l_j\right)}^2-2{\left(x/l_j\right)}^3,\hspace{1em}{N}_6=x\left[{\left(x/l_j\right)}^2-x/l_j\right];\end{array}$$

where x is the horizontal distance from the left node of the element; $l_j(j=\mathrm{b},\mathrm{r})$ is the length of the beam element.

The axial and vertical displacements of the bridge element related to the nodal DOFs can be respectively expressed by:

$$v_{\mathrm{b}}\left(x\right)={\mathbf{N}}_{u,\mathrm{b}}{\mathbf{d}}_{\mathrm{b}}^{\mathrm{T}}$$

(15)

$$u_{\mathrm{b}}\left(x\right)={\mathbf{N}}_{v,\mathrm{b}}{\mathbf{d}}_{\mathrm{b}}^{\mathrm{T}}$$

(16)

The axial and vertical displacements of the rail element related to the nodal DOFs can be respectively expressed by:

$$v_{\mathrm{r}}\left(x\right)={\mathbf{N}}_{u,\mathrm{r}}{\mathbf{d}}_{\mathrm{r}}^{\mathrm{T}}$$

(17)

$$u_{\mathrm{r}}\left(x\right)={\mathbf{N}}_{v,\mathrm{r}}{\mathbf{d}}_{\mathrm{r}}^{\mathrm{T}}$$

(18)

Substituting the displacement fields, i.e., Equations (15)–(18), into the virtual work Equation (1) yields the bridge and bridge–rail coupling submatrices. The bridge submatrix is marked with the subscript “bb”. The mass submatrix ${\mathbf{M}}_{\mathrm{bb}}$ of the bridge with the order $n_{\mathrm{b}}\times n_{\mathrm{b}}$ can be expressed by:

$${\mathbf{M}}_{\mathrm{bb}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{M}}_{\mathrm{b}1}& {\mathbf{M}}_{\mathrm{b}2}& \cdots & {\mathbf{M}}_{\mathrm{b}{n}_{\mathrm{b}}}\end{array}\right\}$$

(19)

with

$${\mathbf{M}}_{\mathrm{b}i}={\rho }_{\mathrm{b}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{A}_{\mathrm{b}}(x){\mathbf{N}}_{u,\mathrm{b}}^{\mathrm{T}}{\mathbf{N}}_{u,\mathrm{b}}\mathrm{d}x+{\rho }_{\mathrm{b}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{A}_{\mathrm{b}}(x){\mathbf{N}}_{v,\mathrm{b}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{b}}\mathrm{d}x$$

The stiffness submatrix ${\mathbf{K}}_{\mathrm{bb}}$ of the bridge with the order $n_{\mathrm{b}}\times n_{\mathrm{b}}$ can be expressed by:

$${\mathbf{K}}_{\mathrm{bb}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{K}}_{\mathrm{b}1}& {\mathbf{K}}_{\mathrm{b}2}& \cdots & {\mathbf{K}}_{\mathrm{b}{n}_{\mathrm{b}}}\end{array}\right\}$$

(20)

with

$${\mathbf{K}}_{\mathrm{b}i}=E_{\mathrm{b}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{A}_{\mathrm{b}}(x){{\mathbf{N}}^{\prime }}_{v,\mathrm{b}}^{\mathrm{T}}{{\mathbf{N}}^{\prime }}_{v,\mathrm{b}}\mathrm{d}x+E_{\mathrm{b}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{I}_{\mathrm{b}}(x){{\mathbf{N}}^{″}}_{v,\mathrm{b}}^{\mathrm{T}}{{\mathbf{N}}^{″}}_{v,\mathrm{b}}\mathrm{d}x+k_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{\mathbf{N}}_{v,\mathrm{b}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{b}}\mathrm{d}x$$

The damping submatrix ${\mathbf{C}}_{\mathrm{bb}}$ of the bridge with the order $n_{\mathrm{b}}\times n_{\mathrm{b}}$ can be expressed by:

$${\mathbf{C}}_{\mathrm{bb}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{C}}_{\mathrm{b}1}& {\mathbf{C}}_{\mathrm{b}2}& \cdots & {\mathbf{C}}_{\mathrm{b}{n}_{\mathrm{b}}}\end{array}\right\}$$

(21)

with

$${\mathbf{C}}_{\mathrm{b}i}=c_{\mathrm{b}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{\mathbf{N}}_{u,\mathrm{b}}^{\mathrm{T}}{\mathbf{N}}_{u,\mathrm{b}}\mathrm{d}x+c_{\mathrm{b}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{\mathbf{N}}_{v,\mathrm{b}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{b}}\mathrm{d}x+c_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{\mathbf{N}}_{v,\mathrm{b}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{b}}\mathrm{d}x$$

The bridge–rail coupling submatrices are marked with the subscripts “br” and “rb”. The stiffness submatrices ${\mathbf{K}}_{\mathrm{br}}$ with the order $n_{\mathrm{b}}\times n_{\mathrm{r}}$ and ${\mathbf{K}}_{\mathrm{rb}}$ with the order $n_{\mathrm{r}}\times n_{\mathrm{b}}$ can be expressed by:

$${\mathbf{K}}_{\mathrm{br}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{K}}_{\mathrm{br}1}& {\mathbf{K}}_{\mathrm{br}2}& \cdots & {\mathbf{K}}_{\mathrm{br}{n}_{\mathrm{b}}}\end{array}\right\}$$

(22)

$${\mathbf{K}}_{\mathrm{rb}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{K}}_{\mathrm{rb}1}& {\mathbf{K}}_{\mathrm{rb}2}& \cdots & {\mathbf{K}}_{\mathrm{br}{n}_{\mathrm{r}}}\end{array}\right\}$$

(23)

with

$$\begin{gathered} {\mathbf{K}}_{\mathrm{br}i}=-k_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{\mathbf{N}}_{v,\mathrm{b}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{r}}\mathrm{d}x \\ {\mathbf{K}}_{\mathrm{rb}i}=-k_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{\mathbf{N}}_{v,\mathrm{r}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{b}}\mathrm{d}x \end{gathered}$$

The damping submatrices ${\mathbf{C}}_{\mathrm{br}}$ with the order $n_{\mathrm{b}}\times n_{\mathrm{r}}$ and ${\mathbf{C}}_{\mathrm{rb}}$ with the order $n_{\mathrm{r}}\times n_{\mathrm{b}}$ can be expressed by:

$${\mathbf{C}}_{\mathrm{br}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{C}}_{\mathrm{br}1}& {\mathbf{C}}_{\mathrm{br}2}& \cdots & {\mathbf{C}}_{\mathrm{br}{n}_{\mathrm{b}}}\end{array}\right\}$$

(24)

$${\mathbf{C}}_{\mathrm{rb}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{C}}_{\mathrm{rb}1}& {\mathbf{C}}_{\mathrm{rb}2}& \cdots & {\mathbf{C}}_{\mathrm{br}{n}_{\mathrm{r}}}\end{array}\right\}$$

(25)

with

$$\begin{gathered} {\mathbf{C}}_{\mathrm{br}i}=-c_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{\mathbf{N}}_{v,\mathrm{b}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{r}}\mathrm{d}x \\ {\mathbf{C}}_{\mathrm{rb}i}=-c_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{\mathbf{N}}_{v,\mathrm{r}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{b}}\mathrm{d}x \end{gathered}$$

3.2. Submatrix of Rail

Substituting the displacement fields, i.e., Equations (15)–(18), into the virtual work Equation (4) yields the rail submatrix. The rail submatrix is marked with the subscript “rr”. The mass submatrix ${\mathbf{M}}_{\mathrm{rr}}$ of the rail with the order $n_{\mathrm{r}}\times n_{\mathrm{r}}$ can be expressed by:

$${\mathbf{M}}_{\mathrm{rr}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{M}}_{\mathrm{r}1}& {\mathbf{M}}_{\mathrm{r}2}& \cdots & {\mathbf{M}}_{\mathrm{r}{n}_{\mathrm{r}}}\end{array}\right\}$$

(26)

with

$${\mathbf{M}}_{\mathrm{r}i}={\rho }_{\mathrm{r}}{A}_{\mathrm{r}}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{\mathbf{N}}_{u,\mathrm{r}}^{\mathrm{T}}{\mathbf{N}}_{u,\mathrm{r}}\mathrm{d}x+{\rho }_{\mathrm{r}}{A}_{\mathrm{r}}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{\mathbf{N}}_{v,\mathrm{r}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{r}}\mathrm{d}x$$

The stiffness submatrix ${\mathbf{K}}_{\mathrm{rr}}$ of the rail with the order $n_{\mathrm{r}}\times n_{\mathrm{r}}$ can be expressed by:

$${\mathbf{K}}_{\mathrm{rr}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{K}}_{\mathrm{r}1}& {\mathbf{K}}_{\mathrm{r}2}& \cdots & {\mathbf{K}}_{\mathrm{r}{n}_{\mathrm{r}}}\end{array}\right\}$$

(27)

with

$${\mathbf{K}}_{\mathrm{r}i}=E_{\mathrm{r}}{A}_{\mathrm{r}}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{{\mathbf{N}}^{\prime }}_{v,\mathrm{r}}^{\mathrm{T}}{{\mathbf{N}}^{\prime }}_{v,\mathrm{r}}\mathrm{d}x+E_{\mathrm{r}}{I}_{\mathrm{r}}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{{\mathbf{N}}^{″}}_{v,\mathrm{r}}^{\mathrm{T}}{{\mathbf{N}}^{″}}_{v,\mathrm{r}}\mathrm{d}x+k_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{\mathbf{N}}_{v,\mathrm{r}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{r}}\mathrm{d}x+{\displaystyle \sum _{k=0}^{n_{\mathrm{w}}}{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}{\mathbf{N}}_{v,\mathrm{r}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{r}}}$$

The damping submatrix ${\mathbf{C}}_{\mathrm{rr}}$ of the rail with the order $n_{\mathrm{r}}\times n_{\mathrm{r}}$ can be expressed by:

$${\mathbf{C}}_{\mathrm{rr}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{C}}_{\mathrm{r}1}& {\mathbf{C}}_{\mathrm{r}2}& \cdots & {\mathbf{C}}_{\mathrm{r}{n}_{\mathrm{r}}}\end{array}\right\}$$

(28)

with

$${\mathbf{C}}_{\mathrm{r}i}=c_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{r}}}{\mathbf{N}}_{v,\mathrm{r}}^{\mathrm{T}}{\mathbf{N}}_{v,\mathrm{r}}\mathrm{d}x$$

3.3. Submatrices of Heavy-Haul Train and Train–Rail Coupling

The displacement vector of each heavy-haul train comprises six degrees of freedom, which can be expressed by:

$${\mathbf{d}}_{\mathrm{v}}=\left\{\begin{array}{cccccc}{v}_{\mathrm{c}}& {\theta }_{\mathrm{c}}& v_{\mathrm{t}1}& {\theta }_{\mathrm{t}1}& v_{\mathrm{t}2}& {\theta }_{\mathrm{t}2}\end{array}\right\}$$

(29)

The mass matrix ${\mathbf{M}}_{\mathrm{vv}}$, stiffness matrix ${\mathbf{K}}_{\mathrm{vv}}$ and damping matrix ${\mathbf{C}}_{\mathrm{vv}}$ of the heavy-haul train can be derived from the virtual work Equations (5)–(7). The mass matrix ${\mathbf{M}}_{\mathrm{vv}}$ with the order $n_{\mathrm{v}}\times n_{\mathrm{v}}$ can be expressed by:

$${\mathbf{M}}_{\mathrm{vv}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{M}}_{\mathrm{v}1}& {\mathbf{M}}_{\mathrm{v}2}& \cdots & {\mathbf{M}}_{\mathrm{v}{n}_{\mathrm{v}}}\end{array}\right\}$$

(30)

with

$${\mathbf{M}}_{\mathrm{v}i}=\mathrm{diag}\left\{\begin{array}{cccccc}{m}_{\mathrm{c}}& I_{\mathrm{c}}& m_{\mathrm{t}}+2m_{\mathrm{w}}& I_{\mathrm{t}}+2m_{\mathrm{w}}{l}_{\mathrm{t}}^2& m_{\mathrm{t}}+2m_{\mathrm{w}}& I_{\mathrm{t}}+2m_{\mathrm{w}}{l}_{\mathrm{t}}^2\end{array}\right\}$$

The stiffness matrix ${\mathbf{K}}_{\mathrm{vv}}$ with the order $n_{\mathrm{v}}\times n_{\mathrm{v}}$ can be expressed by:

$${\mathbf{K}}_{\mathrm{vv}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{K}}_{\mathrm{v}1}& {\mathbf{K}}_{\mathrm{v}2}& \cdots & {\mathbf{K}}_{\mathrm{v}{n}_{\mathrm{v}}}\end{array}\right\}$$

(31)

with

$$\begin{array}{l}{\mathbf{K}}_{\mathrm{v}i}=\left[\begin{array}{cccccc}2k_2& 0& -k_2& 0& -k_2& 0\\ & 2k_2{L}_{\mathrm{c}}^2& -k_2{L}_{\mathrm{c}}& 0& k_2{L}_{\mathrm{c}}& 0\\ & & k_2+{\displaystyle \sum _{k=1}^2{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}}& {\displaystyle \sum _{k=1}^2{\eta }_k{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}{l}_{\mathrm{t}}}& 0& 0\\ & & & {\displaystyle \sum _{k=1}^2{a}_{\mathrm{w}k}{k}_{\mathrm{w}r}{l}_{\mathrm{t}}^2}& 0& 0\\ & & symm.& & k_2+{\displaystyle \sum _{k=3}^4{a}_{\mathrm{w}k}{k}_{\mathrm{w}r}}& {\displaystyle \sum _{k=3}^4{\eta }_k{a}_{\mathrm{w}k}{k}_{\mathrm{w}r}{l}_{\mathrm{t}}}\\ & & & & & {\displaystyle \sum _{k=3}^4{a}_{\mathrm{w}k}{k}_{\mathrm{w}r}{l}_{\mathrm{t}}^2}\end{array}\right]\\ {\eta }_k=\left\{\begin{array}{c}1\hspace{1em}k=1,3\\ -1\hspace{1em}k=2,4\end{array}\right.\end{array}$$

The damping matrix ${\mathbf{C}}_{\mathrm{vv}}$ with the order $n_{\mathrm{v}}\times n_{\mathrm{v}}$ can be expressed by:

$${\mathbf{C}}_{\mathrm{vv}}=\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{C}}_{\mathrm{v}1}& {\mathbf{C}}_{\mathrm{v}2}& \cdots & {\mathbf{C}}_{\mathrm{v}{n}_{\mathrm{v}}}\end{array}\right\}$$

(32)

with

$${\mathbf{C}}_{\mathrm{v}i}=\left[\begin{array}{cccccc}2c_2& 0& -c_2& 0& -c_2& 0\\ & 2c_2{L}_{\mathrm{c}}^2& -c_2{L}_{\mathrm{c}}& 0& c_2{L}_{\mathrm{c}}& 0\\ & & c_2& 0& 0& 0\\ & & & 0& 0& 0\\ & & symm.& & c_2& 0\\ & & & & & 0\end{array}\right]$$

The load vector ${\mathbf{F}}_{\mathrm{v}}$ of the train with the order $n_{\mathrm{v}}\times 1$, which contains the vector of vehicle gravity ${\mathbf{F}}_{\mathrm{v},\mathrm{g}}$ and additional force ${\mathbf{F}}_{\mathrm{v},\mathrm{a}}$ induced by track irregularity, can be expressed by:

$${\mathbf{F}}_{\mathrm{v}}={\mathbf{F}}_{\mathrm{v},\mathrm{g}}+{\mathbf{F}}_{\mathrm{v},\mathrm{a}}$$

(33)

with

$$\begin{array}{l}{\mathbf{F}}_{\mathrm{v},\mathrm{g}}={\left\{\begin{array}{cccccccc}0& m_{\mathrm{c}}g& 0& (m_{\mathrm{t}}+2m_{\mathrm{w}})g& 0& (m_{\mathrm{t}}+2m_{\mathrm{w}})g& 0& 0\end{array}\right\}}^{\mathrm{T}}\\ {\mathbf{F}}_{\mathrm{v},\mathrm{a}}={\left\{\begin{array}{cccccccc}0& 0& 0& {\displaystyle \sum _{k=1}^2{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}I{R}_{yk}}& {\displaystyle \sum _{k=1}^2{\eta }_k{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}{l}_{\mathrm{t}}I{R}_{yk}}& {\displaystyle \sum _{k=3}^4{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}I{R}_{yk}}& {\displaystyle \sum _{k=3}^4{\eta }_k{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}{l}_{\mathrm{t}}I{R}_{yk}}& 0\end{array}\right\}}^{\mathrm{T}}\end{array}$$

The train–rail coupling submatrices are marked with the subscripts “vr” and “rv”. The stiffness submatrices ${\mathbf{K}}_{\mathrm{vr}}$ with the order $n_{\mathrm{v}}\times n_{\mathrm{r}}$ and ${\mathbf{K}}_{\mathrm{rv}}$ with the order $n_{\mathrm{r}}\times n_{\mathrm{v}}$ can be expressed by:

$${\mathbf{K}}_{\mathrm{vr}}={\mathbf{K}}_{\mathrm{rv}}^{\mathrm{T}}==\mathrm{diag}\left\{\begin{array}{cccc}{\mathbf{K}}_{\mathrm{vr}1}& {\mathbf{K}}_{\mathrm{vr}2}& \cdots & {\mathbf{K}}_{\mathrm{vr}{n}_{\mathrm{v}}}\end{array}\right\}$$

(34)

with

$${\mathbf{K}}_{\mathrm{vr}i}=\left\{\begin{array}{cccccc}0& 0& {\displaystyle \sum _{k=1}^2{a}_{\mathrm{w}k}{k}_{\mathrm{vr}}{\mathbf{N}}_{v,\mathrm{r}}}& {\displaystyle \sum _{k=1}^2{\eta }_k{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}{l}_t{\mathbf{N}}_{v,\mathrm{r}}}& {\displaystyle \sum _{k=3}^4{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}{\mathbf{N}}_{v,\mathrm{r}}}& {\displaystyle \sum _{k=3}^4{\eta }_k{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}{l}_{\mathrm{t}}{\mathbf{N}}_{v,\mathrm{r}}}\end{array}\right\}$$

The load vector ${\mathbf{F}}_{\mathrm{r}}$ of the rail induced by track irregularities with the order $n_{\mathrm{r}}\times 1$ can be expressed by:

$${\mathbf{F}}_{\mathrm{r}}={\left\{\begin{array}{ccc}0& {\displaystyle \sum _{k=0}^{n_{\mathrm{w}}}{a}_{\mathrm{w}k}{k}_{\mathrm{wr}}{\mathbf{N}}_{v,\mathrm{r}}I{R}_{yk}}& 0\end{array}\right\}}^{\mathrm{T}}(k=1~4)$$

(35)

3.4. Numerical Integration of Coupling Matrix

The passage of heavy-haul trains through the bridge (track) is discretized. The Newmark-β method is employed to solve the train–track–bridge coupling equation at each time step, which assumes that the displacement, velocity and acceleration at previous time $t_k$ have been obtained to calculate the current state of motion at $t_{k+1}$. The equation of motion at time $t_{k+1}$ can be expressed by:

$$\mathbf{M}{\ddot{\mathbf{X}}}_{k+1}+\mathbf{C}{\dot{\mathbf{X}}}_{k+1}+\mathbf{K}{\mathbf{X}}_{k+1}={\mathbf{F}}_{k+1}$$

(36)

By integrating the acceleration from time $t_k$ to $t_{k+1}$, the state of current motion at $t_{k+1}$ can be obtained:

$${\ddot{\mathbf{X}}}_{k+1}=\frac{1}{\beta \mathrm{\Delta }{t}^2}({\mathbf{X}}_{k+1}-{\mathbf{X}}_k)-\frac{1}{\beta \mathrm{\Delta }t}{\dot{\mathbf{X}}}_k-\left(\frac{1}{2\beta }-1\right){\ddot{\mathbf{X}}}_k$$

(37)

$${\dot{\mathbf{X}}}_{k+1}=\frac{\gamma }{\beta \mathrm{\Delta }t}({\mathbf{X}}_{k+1}-{\mathbf{X}}_k)-\left(1-\frac{\gamma }{\beta }\right){\dot{\mathbf{X}}}_k+\left(\frac{\gamma }{2\beta }-1\right){\ddot{\mathbf{X}}}_k\mathrm{\Delta }t$$

(38)

$${\mathbf{X}}_{k+1}={\widehat{\mathbf{K}}}^{-1}{\widehat{\mathbf{F}}}_{k+1}$$

(39)

with

$$\widehat{\mathbf{K}}=\mathbf{K}+\frac{1}{\beta \mathrm{\Delta }{t}^2}\mathbf{M}+\frac{1}{\beta \mathrm{\Delta }t}\mathbf{C}$$

$${\widehat{\mathbf{F}}}_{k+1}={\mathbf{F}}_{k+1}+\left[\frac{1}{\beta \mathrm{\Delta }{t}^2}{\mathbf{X}}_k+\frac{1}{\beta \mathrm{\Delta }t}{\dot{\mathbf{X}}}_k+\left(\frac{1}{2\beta }-1\right){\ddot{\mathbf{X}}}_k\right]\mathbf{M}+\left[\frac{\gamma }{\beta \mathrm{\Delta }t}{\mathbf{X}}_k+\left(\frac{\gamma }{\beta }-1\right){\dot{\mathbf{X}}}_k+\frac{\mathrm{\Delta }t}{2}\left(\frac{\gamma }{\beta }-2\right){\ddot{\mathbf{X}}}_k\right]\mathbf{C}$$

When γ = 1/2 and β = 1/4, the Newmark-β method is unconditionally stable. The dynamic responses of the heavy-haul train–track–bridge coupling system at discrete time intervals can be solved by Equations (37)–(39).

From the above analysis, it can be seen that the coupling between the track and the bridge remains time-invariant, while the coupling between the heavy-haul train and the track exhibits time-varying characteristics. The time-invariant submatrices are initially assembled in the solution, followed by the establishment of the integral heavy-haul train–track–bridge coupling matrix in each integration, taking into account the time step and the varying component of wheel–rail contact stiffness. Finally, the dynamic equations are solved by the Newmark-β method. The calculation flow of the heavy-haul train–track–bridge coupling system is shown in Figure 3, and the computational program is compiled by MATLAB programming language.

4. Case Study

4.1. Project Background

The Yellow River Bridge, located on a heavy-haul railway line, is the subject of investigation, as shown in Figure 4. The bridge is a single-track structure designed for a speed of 80 km/h, with a total length of 325.4 m. The bridge structure is a continuous rigid frame bridge featuring a single-hole box girder with a variable cross-section. The spans of the left, middle and right are 96.7 m, 132.0 m and 96.7 m, respectively. The ballasted track is laid on the bridge with a rail mass of 60 kg/m per unit length.

The side view model of the bridge is depicted in Figure 5a. There are four sections of variable cross-section with parabolic heights, each 55.0 m in length on both sides of Piers 2 and 3, as shown in Figure 5b. The variable section of the bridge is divided into 16 elements at unequal intervals based on the length of each construction section. The main beam consists of 17 key sections. The cross-section is shown in Figure 5c. The inertia moment and area parameters for 17 key sections are presented in

Table 1. Parameters of the key sections of Yellow River Bridge.

Section	1	2	3	4	5	6	7	8	9
Area (m2)	23.92	22.74	21.95	21.19	20.48	15.96	15.38	14.76	14.20
Inertia moment (m4)	247.47	210.91	188.35	168.38	150.74	120.48	106.24	92.18	80.50
Section	10	11	12	13	14	15	16	17
Area (m2)	13.62	13.13	12.71	12.36	12.09	10.65	10.55	10.52
Inertia moment (m4)	69.61	60.95	54.17	48.98	45.16	40.60	39.14	38.66

.

In order to effectively release the axial displacement of the bridge, a consolidation mode is adopted at the pier–beam connection, and the other two piers are provided with supports. The box girder is constructed using C55 concrete, which has a density of 26.5 kN/m³ and an elastic modulus of 35.5 GPa. The pier body is constructed using C45 concrete, which has an elastic modulus of 34.5 GPa and a Poisson’s ratio of 0.2.

The following parameters are assumed for the track: $E_{\mathrm{r}}=2.06\times {10}^{11}\mathrm{Pa}$, $I_{\mathrm{r}}=3.217\times {10}^{-5}{\mathrm{m}}^4$, $m_{\mathrm{r}}=60.64 \mathrm{kg}/\mathrm{m}$, $k_{\mathrm{rb}}=6\times {10}^7\mathrm{N}/\mathrm{m}$ and $c_{\mathrm{rb}}=7.5\times {10}^4\mathrm{N}\cdot \mathrm{s}/\mathrm{m}$. In addition, the following parameters are assumed for the heavy-haul train: $m_{\mathrm{c}}=9\times {10}^4\mathrm{kg}$, $I_{\mathrm{c}}=1.4075\times {10}^6\mathrm{kg}\cdot {\mathrm{m}}^2$, $m_{\mathrm{t}}=460 \mathrm{kg}$, $I_{\mathrm{t}}=175 \mathrm{kg}\cdot {\mathrm{m}}^2$, $L_{\mathrm{c}}=4.6 \mathrm{m}$, $l_{\mathrm{t}}=0.915 \mathrm{m}$, $L_{\mathrm{v}}=12 \mathrm{m}$ (total length of each vehicle), $k_2=1.9\times {10}^6\mathrm{N}/\mathrm{m}$, $c_2=3\times {10}^6\mathrm{N}\cdot \mathrm{s}/\mathrm{m}$ and $k_{\mathrm{wr}}=2.8\times {10}^8\mathrm{N}/\mathrm{m}$.

4.2. Dynamic Analysis on the Effect of Track Irregularities

The presence of geometric irregularities on the track is a significant factor contributing to the excitation of vibration in the heavy-haul train–track–bridge coupling system. The wavelength range of track irregularities is extensive, encompassing numerous harmonic components with varying amplitudes and wavelengths. The common wavelength range of the track irregularity is 0.01~200 m, in which the short and medium wave have great effects on the vibration of the coupling system. In order to comprehensively capture the vibrations of trains, tracks, and bridges caused by track geometries, track irregularities of different wavelengths must be included in the dynamic analysis.

The heavy-haul train–track–bridge coupling dynamic model, based on the virtual work principle, incorporates a time domain input as the system excitation. Inverse Fourier transform (IFFT) is used to transform the power spectral density function of the track’s random irregularity into a spatial sample that varies with the distance traveled along the track. The low interference spectrum is used to generate longitudinal irregularity with a wavelength range of 1 to 10 m, as shown in Figure 6. It can be seen that the simulation of the power spectral density is consistent with the analysis, indicating that the generated track irregularity sample is reasonable and effective.

To better verify the effectiveness of the proposed model, a vibration test of the middle span of the Yellow River Bridge was carried out. The displacement and acceleration obtained from the proposed model and the experimental test are compared in Figure 7. It can be seen that the calculated displacement and acceleration are consistent with the measured results in terms of trend. Among them, the calculated absolute maximum displacement and acceleration of the bridge are 16.96 mm and 1.56 m/s², respectively, while they are 17.13 mm and 1.95 m/s² in the field test. The discrepancy of displacement in amplitudes is limited, ranging from 1% to 11%. However, there exists a slightly larger disparity between the simulated and measured acceleration. One reason may be that the track irregularities in other wavelength ranges are sensitive to bridge acceleration. Other primary factors contributing to the inconsistency include the inconsistent axle loads, structural fatigue damages, uncertain measurement noises and environmental changes.

The vertical displacement and acceleration of the rail in the transition section are shown in Figure 8. Ballast can be clearly observed to effectively attenuate the free vibration caused by axle loads of the heavy-haul train. The heavy-haul train lacking primary suspension fails to effectively mitigate the impact of the nodding motion of the bogie on the track.

Figure 9 shows the contact force between the front wheelset and the rear wheelsets of both the front and rear bogies of a vehicle, from which one can notice that the variation in wheelset contact force is more violent in the rear bogie compared to the front bogie. The fluctuation of contact force of the rear wheelset significantly increases, indicating that the nodding motion of the front bogie has a significant effect on the contact force of the rear wheelset.

The maximum dynamic responses of the system in the four wavelength ranges, i.e., 1~10 m, 1~20 m, 1~30 m and 1~40 m, are shown in Figure 10. The track irregularities in the wavelength range of 1 to 20 m have a significant effect on the responses of the bridge, whereas wavelengths greater than 20 m show a negligible influence on the vibration of the bridge. The resonance of the car body occurs when the wavelength exceeds 20 m, resulting in a substantial amplification of both the acceleration amplitude of the car body and the wheel–rail contact force. The influence of the interaction between the heavy-haul train and the rail on the track vibration is significantly more pronounced compared to that of the interaction between the bridge and the track.

4.3. Dynamic Analysis on the Effect of Axle Load

Increasing the axle loads of a heavy-haul train can significantly improve the efficiency of transportation. In order to investigate the influence of axle load on the dynamic behavior of system, three kinds of axle loads of 25 t, 27.5 t and 30 t are studied. In this section, calculations still consider track irregularities with wavelengths ranging from 1 to 10 m.

Figure 11 and Figure 12 show the effects of mid-span displacement and acceleration of the bridge, as well as car body displacement and acceleration under track irregular excitation, respectively. From these results, it can be seen that the displacement response of the system increases linearly with the increase in axle load. When the freight volume of heavy-haul trains increases, the displacement can be reduced by strengthening measures to improve the stiffness of the structure. In Figure 11b, the increase of axle load has a minimal impact on bridge acceleration. The coupling vibration and the irregular excitation of the track are the main factors affecting bridge acceleration. As shown in Figure 12b, the increase in the mass of the car body can restrain the peak acceleration of the car body to a certain extent. The irregular track excitation has a significant effect on the acceleration of the car body if the freight volume of the vehicle is small. The utilization of a new heavy-haul train equipped with primary suspension can effectively enhance the vibration isolation performance against track irregular excitation, thereby mitigating the acceleration of the car body.

4.4. Dynamic Analysis on the Effect of Bridge Structural Degradation

Under the cyclic loading of heavy-haul trains, structural fatigue and performance degradation of railway bridges are inevitable. Rayleigh damping is adopted into the bridge model. Considering track irregularities with wavelengths ranging from 1 to 10 m as excitation, a degradation of overall stiffness and damping by 10% and 20%, respectively, is taken into account.

Figure 13 and Figure 14 plot the responses of the displacement and acceleration in the middle span of the bridge after the stiffness and damping of the bridge are degraded, respectively. From these results, it can be observed that the stiffness degradation of the bridge has the most significant influence on the displacement of the bridge. As shown in Figure 13a and Figure 14a, the stiffness of the bridge is reduced by 10% and 20%, respectively, while the mid-span displacement of the bridge is increased by 11% and 25%, respectively. However, the decrease in bridge damping can be disregarded when considering the displacement in the middle of the bridge span. As shown in Figure 13b and Figure 14b, the decrease in bridge stiffness and damping leads to an increase in the amplitude of acceleration. However, the magnitude of the change is not obvious. Therefore, the fatigue cracking of the bridge occurs due to the cyclic loading from heavy-haul trains, resulting in a reduction in bridge stiffness. The excessive deformation further deteriorates the mechanical performance of the bridge and exacerbates its dynamic behavior.

5. Conclusions

A heavy-haul train–track–bridge coupling system, based on the virtual work principle and represented by a coupling matrix, is presented in this paper. The research focuses on a continuous rigid frame bridge accommodating a heavy-haul railway line as the object of study. The validity and accuracy of the proposed model are verified through comparison with the model and experimental results. The following major conclusions are drawn from this study:

(a)

According to the vibration characteristics of the interaction between the heavy-haul train, track and bridge, the motion equation of the train–track–bridge coupling system is derived by the principle of virtual work, which is combined with the displacement field to obtain the expression of the coupling matrix of the motion equation of the whole system. The dynamic responses of the heavy-haul train, track and bridge are solved by a step-by-step integration method at the same time;

(b)

Notably, the contact force fluctuations on the rear wheelset of heavy-haul trains are considerably higher than those on the front wheelset within the same bogie. This discrepancy arises due to the absence of primary suspension isolating the nodding motion of the bogie. Additionally, the bridge’s vibration is particularly sensitive to track irregularities with wavelengths ranging from 1 to 20 m, especially within the 1 to 10 m range;

(c)

The displacement of both the bridge and the car body increases linearly as the load on the car body increases, while the acceleration of the bridge and car body have negligible influence;

(d)

The study incorporates the impact of performance degradation on bridge vibration by reducing the bridge’s stiffness and damping by 10% and 20%, respectively. This attenuation significantly affects the vertical displacement of the bridge, whereas the impact on vertical displacement due to damping attenuation can be disregarded. Both the stiffness and damping attenuation contribute to a subtle increase in the bridge’s acceleration.