An Effective Tangent Stiffness of Train–Track–Bridge Systems Based on Artiﬁcial Neural Network

: Dynamic response analysis of a train–track–bridge (TTB) system is a challenging task for researchers and engineers, partially due to the complicated nature of the wheel–rail interaction (WRI). When Newton’s method is used to solve implicit nonlinear ﬁnite element equations of a TTB system, consistent tangent stiffness (CTS) is essential to guarantee the quadratic convergence rate. However, the derivation and software implementation of CTS for the WRI element require signiﬁcant efforts. Artiﬁcial neural network (ANN) can directly obtain a potentially good tangent stiffness by a trained relationship between input nodal displacement/velocity and output tangent stiffness. In this paper, the backpropagation neural-network-based tangent stiffness (BPTS) of the WRI element is derived and implemented into a general ﬁnite element software, OpenSees, and veriﬁed by dynamic response analysis of a high-speed train running on a seven span simply supported beam bridge. The accuracy and efﬁciency are compared between the use of BPTS and CTS. The results demonstrate that BPTS can not only save the signiﬁcant efforts of deriving and software implementing CTS but also improve computational efﬁciency while ensuring good accuracy. discusses a simple working condition with a TTB system, and the load on the train is relatively simple, that is, a small lateral disturbance. The application of BPTS under complex working conditions will be studied in the future, such as the effect of track irregularity and wind load, in order to improve the generality of this method.


Introduction
Dynamic response analysis of a train-track-bridge (TTB) system is a challenging task for researchers and engineers, partially due to the strong nonlinearity and complex contact conditions of the wheel-rail interaction (WRI). In recent years, many researchers have proposed different WRI models based on various methods. Zhai et al. [1][2][3][4] developed a novel explicit two-step method and a new family of predictor-corrector integration algorithms used in the solving of dynamic problems, based on which they developed a computer program named TTISIM to predict the vertical and lateral dynamic responses of the vehicle-track coupled system. Montenegro et al. [5] proposed a wheel-rail contact formulation for analyzing the nonlinear train-structure interaction that takes into account wheel and rail geometry and used a WRI element to model the behavior of the contact interface based on Hertz's theory and Kalker's laws. Pombo [6] presented a general wheel-rail contact detection formulation and calculated the contact forces using Kalker linear theory, a heuristic nonlinear model, and the Polach formulation. Liu et al. [7] developed a 3D WRI element by wrapping the WRI into an element where the contact force is obtained by using Kalker linear theory, nonlinear Hertz theory, and the Polach formulation. The WRI element has been demonstrated to be efficient and accurate, able to simulate complicated TTB systems, and is therefore used in the dynamic responses analysis of the TTB systems in this paper.
When Newton's method is used to solve implicit nonlinear finite element equations, e.g., the TTB system, the consistent tangent stiffness (CTS) is essential to guarantee the quadratic convergence rate [8][9][10][11][12]. For the WRI element mentioned above, the computation of the CTS requires the differentiation of the internal forces with respect to the nodal quadratic convergence rate [8][9][10][11][12]. For the WRI element mentioned above, the computation of the CTS requires the differentiation of the internal forces with respect to the nodal displacement, following the complicated calculation process of the internal forces. Although the CTS is accurate and efficient and is able to guarantee the quadratic rate of convergence of Newton's method [13], its derivation and software implementation require significant efforts. Therefore, it is necessary to find an alternative method for computing the tangent stiffness that can reduce the efforts of derivation and software implementation while maintaining high computational efficiency.
Recently, artificial neural network (ANN) has become increasingly popular in solving a wide range of engineering problems. The ANN can directly obtain a potentially good tangent stiffness by a trained relationship between input nodal displacement/velocity and output tangent stiffness. There are relatively less formula derivation and programming work. The ANN model has been used to solve complex nonlinear problems. Jung et al. [14] presented a rate-dependent ANN constitutive formulation to analyze the time-dependent behaviors of concrete. Luo et al. [15] proposed a novel ANN-based model for fast predicting the backbone curves of nonlinear flexureand shear-critical columns. Huang [16] proposed an ANN-based method for the identification of a discrete-time nonlinear hysteretic system under strong earthquakes. Jung et al. and Kang et al. [17] proposed an ANN-based health monitoring model for predicting the displacement responses of gravity dams. The ANN has been applied to various fields as an accuracy and efficiency method; however, it has not yet been applied to solving the tangent stiffness of complex nonlinear TTB systems.
In this paper, a BP neural-network-based tangent stiffness (BPTS) method is presented for the WRI element proposed by Liu and co-workers [8]. The input of the BP neural network is the displacements/velocities of the wheel node and rail virtual node, while the output is the tangent stiffness matrix of these nodes. The element tangent stiffness is obtained straightforwardly based on the output tangent stiffness. The trained BPTS model is implemented in a general finite element software, OpenSees [18], to substitute the calculation of the consistent tangent stiffness (CTS) [13] and verified by dynamic response analysis of a high-speed train running on a seven-span simply supported beam bridge. The accuracy and efficiency are compared between the use of the BPTS and the CTS. The BPTS is demonstrated to not only be able to save significant efforts of deriving and software implementing the CTS but also to greatly improve computational efficiency while ensuring acceptable accuracy.

The Internal Resisting Force and Consistent Tangent Stiffness of the WRI Element
The WRI element consists of a wheel node and all the rail nodes (b1, b2 … bi …) that the wheel node may pass through, and the crucial part of the WRI element [8] is the socalled active portion, as can be seen in Figure 1. The active portion contains one wheel node and two rail nodes, marked as O w , b i , and b i+1 , respectively. Then, by projecting O w vertically onto the rail element, a so-called 'virtual' rail node O r can be obtained. The coordinates of O r can be obtained by linear interpolation from the coordinates of b i and b i+1 , which can be expressed as where X i j denotes the jth coordinate component of the point i, N L I (I = 1, 2) denotes the shape function of the linear interpolation, N H I (I = 1, 2, 3, 4) denotes the shape function of the Hermite interpolation, and N H I,1 (I = 1, 2, 3, 4) denotes the derivative of the Hermite interpolation shape functions with respect to x. In each time step, the wheel node moves forward, and the shape functions change. After the wheel node O w crosses b i+1 , the active position is updated immediately and consists of O w , b i+1, and b i+2 .
The element internal resisting force R consists of the forces on the nodes O w , b i, and b i+1 , i.e., R= [R w , R b i , R b i+1 ] and can be obtained based on the nodal displacement u and velocity v, i.e., R = R(u, v), where the nodal displacement u consists of the displacements of the wheel node and two rail nodes, i.e., u= [u w , u b i , u b i+1 ], while the velocity v is a function of u, i.e., v = v(u). It is worth mentioning that the function can be explicitly obtained based on the approximation of velocity using displacement when a time-stepping discretization (e.g., Newmark method) is employed. The element internal resisting force R can be calculated based on the contact force F between wheel node O w and virtual rail node O r , which consists of the normal contact force F n and the tangential contact force F t , i.e., where X denotes the coordinates obtained by element nodal displacement u. Suppose there is a virtual interpenetration between the rail and wheel profile, as can be found in Figure 2. The normal contact force F c n and tangential contact forces F c t1 and F c t2 are calculated using nonlinear Hertz theory, Kalker linear theory, and the Polach formulation [19].  After the contact force F is obtained, the resisting force on the wheel node, Rw and the resisting force on the virtual rail node Rr (in Figure 3) can be calculated respectively by

R T A FAT
where A denotes the transformation matrix of the coordinate system, Tw and Tr denote the transformation matrix that converts F to the wheel node and the virtual rail node, which is related to the vector from the wheel-rail contact point to the wheel node and that to the virtual rail node [7]. The superscript T denotes the transpose of the matrix. The resisting force on the two rail nodes,  After the contact force F is obtained, the resisting force on the wheel node, R w and the resisting force on the virtual rail node R r (in Figure 3) can be calculated respectively by where A denotes the transformation matrix of the coordinate system, T w and T r denote the transformation matrix that converts F to the wheel node and the virtual rail node, which is related to the vector from the wheel-rail contact point to the wheel node and that to the virtual rail node [7]. The superscript T denotes the transpose of the matrix. The resisting force on the two rail nodes, R b i and R b i+1 , are calculated using the force on the virtual node R r based on their displacement relation in Equation (1), i.e., using the virtual work principle, as shown in Figure 3.
After the contact force F is obtained, the resisting force on the wheel node, Rw and the resisting force on the virtual rail node Rr (in Figure 3) can be calculated respectively by

R T A FAT
where A denotes the transformation matrix of the coordinate system, Tw and Tr denote the transformation matrix that converts F to the wheel node and the virtual rail node, which is related to the vector from the wheel-rail contact point to the wheel node and that to the virtual rail node [7]. The superscript T denotes the transpose of the matrix. The resisting force on the two rail nodes, i b R and i b R 1 + , are calculated using the force on the virtual node Rr based on their displacement relation in Equation (1), i.e., using the virtual work principle, as shown in Figure 3. The consistent tangent stiffness (CTS) is important in solving implicit nonlinear finite element equations, which can guarantee the quadratic rate of convergence when the Newton-Raphson algorithm is used. The CTS is calculated by differentiating the element resisting force with nodal displacement using the chain rule of derivation, i.e., The consistent tangent stiffness (CTS) is important in solving implicit nonlinear finite element equations, which can guarantee the quadratic rate of convergence when the Newton-Raphson algorithm is used. The CTS is calculated by differentiating the element resisting force with nodal displacement using the chain rule of derivation, i.e.,

∂R ∂u
where ∂R w ∂u and ∂R r ∂u are calculated as where the subscript "i" can be w or r, denoting the force of wheel and rail, respectively. The derivation and implementation of Equation (5) require significant efforts, and the detailed derivation formula can be found in [13].

The BP-Based Tangent Stiffness (BPTS)
A backpropagation (BP) neural network is used for obtaining an alternative tangent stiffness. The BP neural network consists of an input layer, a hidden layer, and an output layer, as can be seen in Figure 4. Each layer is made up of a few neurons. The number of neurons in the input and output layers is usually determined by the number of input parameters and output parameters, while the number of neurons in the hidden layer is uncertain. Too few neurons may affect the accuracy of the predicted results, but too many may lead to overfitting. Therefore, multiple tests are often needed to select the best number of hidden layer neurons.
stiffness. The BP neural network consists of an input layer, a hidden layer, and an output layer, as can be seen in Figure 4. Each layer is made up of a few neurons. The number of neurons in the input and output layers is usually determined by the number of input parameters and output parameters, while the number of neurons in the hidden layer is uncertain. Too few neurons may affect the accuracy of the predicted results, but too many may lead to overfitting. Therefore, multiple tests are often needed to select the best number of hidden layer neurons. In the hidden layer, each neuron receives parameters from the input layer and then calculates as follows: In the hidden layer, each neuron receives parameters from the input layer and then calculates as follows: where n indexes the neurons in the hidden layer, and m denotes the number of inputs. h i denotes the ith neuron of the hidden layer, where i ranges from 1 to n. X j is the jth input, where j ranges from 1 to m. w ij is a weighting coefficient for X j 's contribution to h i , while b i is a bias for the ith neuron in the hidden layer. w ij can be assembled into an n m matrix. The f represents a nonlinear activation function, which is the crucial part of a neural network that builds the nonlinear mappings. There are many forms for f, but in this study, the tanh function was used. It can be expressed as where e denotes the basis of the natural logarithms. After each neuron in the hidden layer has been calculated, the neurons in the output layer can be expressed as where M denotes the number of outputs, and N is the number of neurons in the hidden layer. Y j represents the jth output, where j ranges from 1 to M. h k is the kth neuron of the hidden layer, where k ranges from 1 to N. v jk is a weighting coefficient for the contribution of h k to Y j , while d j denotes the bias for the jth output. V jk can be assembled into an M × N matrix. The weighting coefficients and biases (w ij , b i , v jk, and d j ) in Equations (6) and (8) must be obtained by training with a large number of samples. The BP neural network's performance is evaluated by comparing its predicted results, and the training samples are performed by using the mean squared error (MSE) of the discrepancies for quantification. The MSE is where Y j denotes the output calculated by Equation (8), Y j denotes the sample value, and R is the number of the output sample. The MSE is usually guaranteed within a narrow range (0.001 in this study) to get weighting coefficients and biases that ensure the network's performance. In the ANN model for the tangent stiffness of WRI, the inputs X are nodal displacement u and velocity v, while the outputs Y are ∂R w ∂u and ∂R r ∂u in Equation (4). The BPTS of the WRI element can be obtained by the following steps: (1) Selecting an appropriate example, running enough time steps, recording input X and output Y at each iteration of each time step. (2) Selecting an appropriate activation function in Equation (7) and the learning algorithm, determining the number of hidden layers and number of neurons in each layer.
(3) Training the BP neural network model to obtain the weight coefficients and deviations in Equations (6) and (8) (i.e., w ij , b i , v jk , and d j ). (4) Substituting the CTS with the BPTS, i.e., the mapping relationship in Equations (6) and (8), by software implementing the BPTS algorithm in the WRI element. In practice, the calculation to obtain input X and output Y in step 1 and the training process in step 3 can be performed simultaneously, i.e., at the beginning of the calculation, the CTS (or the tangent by the perturbation method that obtains the tangent by the perturbation of displacement [7]) is used, and the input and output obtained are used to train the ANN model. After the ANN model has been well trained, the BPTS is used instead of the CTS to continue the calculation process so as to realize the intelligent choice of the tangent stiffness.

The Model of Train-Rail-Bridge System
A seven-span simply supported beam bridge [20] is analyzed in this study, the model of which is shown in Figure 5a. The single-span length of the bridge is 32.6 m, and the total length is 228.2 m. Both ends are supported on the foundation, and the middle spans are supported by 8-m-high piers. The model in OpenSees is shown in Figure 5b. The track, girder, and piers are simulated by the elasticBeamColumn elements in OpenSees, the parameters of which can be found in Table 1. The simulation of the cross-section of the bridge in OpenSees can be referred to in Figure 6. The tolerance of the convergence is selected as 0.1 N. The analysis type in this paper is algorithm Newton.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 7 of 18 model has been well trained, the BPTS is used instead of the CTS to continue the calculation process so as to realize the intelligent choice of the tangent stiffness.

The Model of Train-Rail-Bridge System
A seven-span simply supported beam bridge [20] is analyzed in this study, the model of which is shown in Figure 5a. The single-span length of the bridge is 32.6 m, and the total length is 228.2 m. Both ends are supported on the foundation, and the middle spans are supported by 8-m-high piers. The model in OpenSees is shown in Figure 5b. The track, girder, and piers are simulated by the elasticBeamColumn elements in OpenSees, the parameters of which can be found in Table 1. The simulation of the cross-section of the bridge in OpenSees can be referred to in Figure 6. The tolerance of the convergence is selected as 0.1 N. The analysis type in this paper is algorithm Newton. The connection between the track and the track plate is simplified by the springdampers, and the rail nodes are connected with the track plate nodes through the twoNodeLink element in OpenSees. The length of the twoNodeLink element is 0.787 m, and the spring and damping are three-dimensional. The bridge bearing is also simulated with the twoNodeLink element, and the element length is 0.2 m. However, since the bridge is simply supported, the bearings of each span are different, as can be seen in Figure 7, where kf, cf represents the stiffness and damping of the fixed bearing, km; cm represents the stiffness and damping of the movable bearing. The parameters of the twoNodeLink element in OpenSees can be seen in Table 2. All the mass of the girder is concentrated on the girder nodes, and the girder node is connected with the track slab node and the upper node of the bearing through the rigidLink beam in OpenSees. The lower node of the bearing is also connected with the upper node of the pier through the      The connection between the track and the track plate is simplified by the springdampers, and the rail nodes are connected with the track plate nodes through the twoN-odeLink element in OpenSees. The length of the twoNodeLink element is 0.787 m, and the spring and damping are three-dimensional. The bridge bearing is also simulated with the twoNodeLink element, and the element length is 0.2 m. However, since the bridge is simply supported, the bearings of each span are different, as can be seen in Figure 7, where k f , c f represents the stiffness and damping of the fixed bearing, k m ; c m represents the stiffness and damping of the movable bearing. The parameters of the twoNodeLink element in OpenSees can be seen in Table 2. All the mass of the girder is concentrated on the girder nodes, and the girder node is connected with the track slab node and the upper node of the bearing through the rigidLink beam in OpenSees. The lower node of the bearing is also connected with the upper node of the pier through the rigidLink beam. The length of the rigidLink beam mentioned above is determined by the actual bridge size. The nodes at the bottom of each pier are fixed.

The Data for the BP Neural Network
As mentioned above, the WRI is the most critical part of the TTB system, while the rest, such as the train, is relatively simple and has a limited effect on the element reaction of the WRI element in terms of efficiency and convergence rate. Therefore, only the CTS of the WRI element is replaced by the BPTS. This will greatly reduce the time spent training the neural network. In this paper, the single wheelset is running on the bridge mentioned in Section 3.1 at a constant speed of 20 m/s, as can be seen in Figure 8. Before the dynamic analysis of the wheelset, a horizontal force is applied in the lateral direction such that an original misalignment of the wheelset is 1 mm. The wheelset herein used is the ML95 [6]; the parameters of the wheelset can be found in Table 3, where M w represents the mass of the wheelset. J xx , J yy , J zz represents inertia about the x, y, and z-axes, respectively. R 0 represents the rolling radius of the wheel. W r and W w represent the distance of the right and left rail and wheel nodes, respectively. E is Young's modulus of the rail and wheel and µ is Poisson's ratio. The analysis type in this paper of the OpenSees model is algorithm Newton. The Newmark-Beta method [21] is used for the time integration, with γ = 0.5 and β = 0.25 with a time step of 0.001 s. The sign for judging convergence in OpenSees is that the maximum unbalanced force at each time step is less than 0.1 N.
As mentioned in Section 2.2, for the ANN model for the tangent stiffness of WRI, the inputs X are the nodal displacement u and velocity v, while the outputs Y are the approximate of ∂R w ∂u and ∂R r ∂u obtained by perturbation [7] because its derivation and programming are quite simple. There are 36 inputs and 216 outputs (decomposed as 12 groups with 18 outputs in each group for fast training purposes). A hidden layer neural network is composed of 20 neurons. For training purposes, the above wheelset example runs 3000 time steps, with an average convergence of 6-7 iterations per time step, so as to obtain a total of 19,927 sample data sets.

The Data for the BP Neural Network
As mentioned above, the WRI is the most critical part of the TTB system, whil rest, such as the train, is relatively simple and has a limited effect on the element rea of the WRI element in terms of efficiency and convergence rate. Therefore, only the of the WRI element is replaced by the BPTS. This will greatly reduce the time s training the neural network. In this paper, the single wheelset is running on the br mentioned in Section 3.1 at a constant speed of 20 m/s, as can be seen in Figure 8. Be the dynamic analysis of the wheelset, a horizontal force is applied in the lateral dire such that an original misalignment of the wheelset is 1 mm. The wheelset herein us the ML95 [6]; the parameters of the wheelset can be found in Table 3

Verification of BPTS in the Wheelset Example
The input of the above wheelset example is used to test the performance of the obtained BP neural network, and the predicted results are compared with the actual results, i.e., by CTS. In order to make the results more visible, the data are sampled and plotted every 50 points. Additionally, only some data comparisons are listed to save space, as can be seen in Figures 9-11, respectively. The data in Figure 9 represent the derivative of the resisting force in the x-direction of the wheel node with respect to the translational displacement in the x-direction of the wheel node, marked as ∂R w1 /∂u 1 . The data in Figure 10 represent the derivative of the resisting force in the y-direction of the rail virtual node to the translational displacement in the y-direction of the wheel node, marked as ∂R r2 /∂u 2 . The data in Figure 11 represent the derivative of the moment of the rail virtual node around the x-direction to the rotational displacement of the rail node b i+1 around the x-direction, marked as ∂R r4 /∂u 16 . It can be observed that the predicted components of tangent stiffness agree well with the actual ones by CTS, although there are irregular jumps of the CTS with iterations. A mean absolute percentage error (MAPE) is used to measure the error between the predicted result and the actual result herein, i.e., where Y represents the actual result and Y represents the result obtained by the BP neural network. The MAPEs of the two curves are 1.14%, 0.34%, and 0.0042%, respectively, which shows that the neural network has a good training effect.    After the training is performed, the CTS is replaced by the tangent stiffness predicted by the BP neural network (i.e., BPTS) and implanted into the above wheelset to verify its responses. Figure 12 shows the comparison of the lateral displacement of the left wheel obtained by two tangent stiffness. Figure 13 shows the comparison of the acceleration of the left wheel obtained by two tangent stiffness. The wheelset response calculated by BPTS agrees well with that calculated by CTS.
To compare the convergence rate by using the two tangent stiffness, the convergence data, i.e., the norm of the unbalance force at each iteration in a few representative time steps, are shown in Table 4. It can be observed that the convergence of the Newton method by use of BPTS is slightly slower than the use of CTS because CTS can guarantee quadratic convergence of the Newton method, while BPTS is obtained by the ANN model that approximates CTS but cannot achieve strict quadratic convergence. However, because the After the training is performed, the CTS is replaced by the tangent stiffness predicted by the BP neural network (i.e., BPTS) and implanted into the above wheelset to verify its responses. Figure 12 shows the comparison of the lateral displacement of the left wheel obtained by two tangent stiffness. Figure 13 shows the comparison of the acceleration of the left wheel obtained by two tangent stiffness. The wheelset response calculated by BPTS agrees well with that calculated by CTS.   To compare the convergence rate by using the two tangent stiffness, the convergence data, i.e., the norm of the unbalance force at each iteration in a few representative time steps, are shown in Table 4. It can be observed that the convergence of the Newton method by use of BPTS is slightly slower than the use of CTS because CTS can guarantee quadratic convergence of the Newton method, while BPTS is obtained by the ANN model that approximates CTS but cannot achieve strict quadratic convergence. However, because the calculation of BPTS is much simpler in each iteration step, it has great advantages in overall computational efficiency. In this example, the calculation time using CTS is 7618 s, while the calculation time using BPTS is 4959 s, which saves 35.0% of the calculation time compared with the former.

A Light Rail Vehicle Running on the Seven-Span Bridge
Using the same trained BPTS above (i.e., that obtained by a wheelset running with a constant speed), an ML95 trailer vehicle running on the bridge is studied to verify the BPTS derived herein, as shown in Figure 14. The detail of the model of the car body and bogie can be found in [22][23][24]. The vehicle consists of the car body, bogies, and wheelsets, respectively. The vehicle body and bogies are modeled using lumped mass points and rigid beams, while the wheelsets are simulated by the wheel nodes mentioned in Section 2.1 and the elastic beam (axis of the wheelsets). The primary suspension of the vehicle is simulated by spring-dampers I between the wheelsets and the bogies, which appear in black. The secondary suspension of the vehicle is simulated by spring-dampers II between the bogies and the car body, which appear in pink. The parameters of spring-dampers I and II can be found in Table 5, where k I denotes the stiffness of spring-dampers I, c I represents the damping of spring-dampers I, k II denotes the stiffness of spring-dampers II, and c II represents the damping of spring-dampers II. Both k and c are three-dimensional herein and correspond to x, y, z in Table 5. The rail profile used herein is the UIC50. The element of rail is simulated by an elastic beam, the parameters of which can be found in Table 6, where A r is the cross-sectional area of the rail, I rx , I ry , I rz represent the area moment of inertia about the x, y, and z-axes, respectively. During the static analysis, the lateral forces applied to each wheelset to make the lateral displacement of the vehicle is 1 mm, and then the forces are removed at the beginning of the dynamic analysis.

A Light Rail Vehicle Running on the Seven-Span Bridge
Using the same trained BPTS above (i.e., that obtained by a wheelset running with a constant speed), an ML95 trailer vehicle running on the bridge is studied to verify the BPTS derived herein, as shown in Figure 14. The detail of the model of the car body and bogie can be found in [22][23][24]. The vehicle consists of the car body, bogies, and wheelsets, respectively. The vehicle body and bogies are modeled using lumped mass points and rigid beams, while the wheelsets are simulated by the wheel nodes mentioned in Section 2.1 and the elastic beam (axis of the wheelsets). The primary suspension of the vehicle is simulated by spring-dampers I between the wheelsets and the bogies, which appear in black. The secondary suspension of the vehicle is simulated by spring-dampers II between the bogies and the car body, which appear in pink. The parameters of spring-dampers I and II can be found in Table 5, where kI denotes the stiffness of spring-dampers I, cI represents the damping of spring-dampers I, kII denotes the stiffness of spring-dampers II, and cII represents the damping of spring-dampers II. Both k and c are three-dimensional herein and correspond to x, y, z in Table 5. The rail profile used herein is the UIC50. The element of rail is simulated by an elastic beam, the parameters of which can be found in Table 6, where Ar is the cross-sectional area of the rail, rx I , ry I , rz I represent the area moment of inertia about the x, y, and z-axes, respectively. During the static analysis, the lateral forces applied to each wheelset to make the lateral displacement of the vehicle is 1 mm, and then the forces are removed at the beginning of the dynamic analysis. The BPTS trained previously using a wheelset is re-used in the TTB system; the lateral displacement of the left wheel of wheelset 1 is shown in Figure 15 for cases of running speed at 30 and 50 m/s, respectively. Additionally, the lateral displacement of the car body is shown in Figure 16. The responses obtained by the two tangent stiffness are identical.   The BPTS trained previously using a wheelset is re-used in the TTB system; the lateral displacement of the left wheel of wheelset 1 is shown in Figure 15 for cases of running speed at 30 and 50 m/s, respectively. Additionally, the lateral displacement of the car body is shown in Figure 16. The responses obtained by the two tangent stiffness are identical.    The convergence data are also recorded. The convergence rates are compared between the use of the two tangent stiffness, as shown in Tables 7 and 8. It can be seen that the convergence rate of BPTS is slightly less than that of CTS. However, the total calculation time is faster when using BPTS over CTS. The total calculation time using CTS is 11,508 s when the speed of the train is 30 m/s, while the calculation time using BPTS is 8264 s, which saves 28% of time. When the train speed is 50 m/s, the calculation times are 16,528 and 13,004 s for the use of CTS and BPTS, respectively, and the use of BPTS saves 21.3% of time. The efficiency improvement by BPTS reduces slightly with the incensement of running speed. The above results show that the BPTS has good applicability in the TTB system. Compared with the CTS, the BPTS can improve computational efficiency and maintain good accuracy. The convergence data are also recorded. The convergence rates are compared between the use of the two tangent stiffness, as shown in Tables 7 and 8. It can be seen that the convergence rate of BPTS is slightly less than that of CTS. However, the total calculation time is faster when using BPTS over CTS. The total calculation time using CTS is 11,508 s when the speed of the train is 30 m/s, while the calculation time using BPTS is 8264 s, which saves 28% of time. When the train speed is 50 m/s, the calculation times are 16,528 and 13,004 s for the use of CTS and BPTS, respectively, and the use of BPTS saves 21.3% of time. The efficiency improvement by BPTS reduces slightly with the incensement of running speed. The above results show that the BPTS has good applicability in the TTB system. Compared with the CTS, the BPTS can improve computational efficiency and maintain good accuracy.

Conclusions
This paper presents a tangent stiffness calculation method based on a backpropagation neural network (i.e., BPTS) for a wheel-rail interaction (WRI) element and implements it into a general finite element software, OpenSees. The BPTS obtained using a wheelset running on a bridge at a constant speed is used for the dynamic response analysis of a realistic train-track-bridge (TTB) system with different vehicle speeds. The accuracy and efficiency of BPTS are studied by comparing the results with those using consistent tangent stiffness (CTS). In the example of a single wheelset, the lateral displacement and acceleration obtained by the two methods are compared and proven to be close. The convergence rate of BPTS is also verified by comparing the norm of the unbalance force and the total computation time. The results show that BPTS needs more convergence numbers but can significantly save the computation time, e.g., when BPTS is used, 35% of time is