An efficient method for train–track–substructure dynamic interaction analysis by implicit-explicit integration and multi-time-step solution

In this work, a method is put forward to obtain the dynamic solution efficiently and accurately for a large-scale train–track–substructure (TTS) system. It is called implicit-explicit integration and multi-time-step solution method (abbreviated as mI-nE-MTS method). The TTS system is divided into train–track subsystem and substructure subsystem. Considering that the root cause of low efficiency of obtaining TTS solution lies in solving the algebraic equation of the substructures, the high-efficient Zhai method, an explicit integration scheme, can be introduced to avoid matrix inversion process. The train–track system is solved by implicitly Park method. Moreover, it is known that the requirement of time step size differs for different sub-systems, integration methods and structural frequency response characteristics. A multi-time-step solution is proposed, in which time step size for the train–track subsystem and the substructure subsystem can be arbitrarily chosen once satisfying stability and precision demand, namely the time spent for m implicit integral steps is equal to n explicit integral steps, i.e., mI = nE as mentioned above. The numerical examples show the accuracy, efficiency, and engineering practicality of the proposed method.

and time-dependence, which must be solved by time-domain step-by-step integral methods. When a large number of time steps are required, the dynamic equations as algebraic forms must be solved many times and the parallel computer is incapable of working well because of the sequential solution of the systems.
Scholars have done pioneering work to address the high-efficient computation problem concerning massive freedoms. In early stage, Hughes et al. [10] and Liu et al. [11] paid attention on putting forward a new family of implicit-explicit finite elements partition procedures. Using this approach, the stiff and flexible subdomains in system interaction are solved independently. Moreover, Belytschko and Mullen [12] discussed the stability in energy of explicitimplicit mesh partitions where the central difference explicit integration and trapezoidal implicit integration methods are considered properly. Smolinski et al. [13] proved the stability of an explicit multi-time-step scheme of the secondorder differential equations. Initiating from [10][11][12][13], multitime-step (MTS) solution has emerged, also called mE-I partition, where m explicit time steps are contained in each implicit time step. Within the Newmark family, Gravouil and Combescure [14] developed a multi-time-step explicitimplicit method for nonlinear structural dynamics. The dual Schur formulation was applied to model the interfaces with interfacial forces represented by Lagrange multipliers. Based on nodal partition method and modified trapezoidal rule, Wu and Smolinski [15] proposed an explicit subcycling integration method. Based on the method presented in Ref. [14], Prakash and Hjelmstad [16] further proposed an MTS coupling method using finite element tearing and interconnecting decomposition [17]. To solve large-scale composite material problems, Beneš et al. [18] introduced subcycling algorithm and using three-filed decomposition method to achieve the implementation of different time steps and possible different numerical methods on various computational domain.
As to the reality of train-track-substructure dynamic interaction, Zhai [19] proposed a new simple explicit two-step method, and by using the integration scheme, the simultaneous solving of algebraic equations is avoided, which is rather suitable for dynamic systems with diagonal mass matrix [5,6].
Based on the precise integration method [20], Zhang et al. [21] proposed a method to obtain accurate and efficient solutions for vehicle-track systems, but it belongs to uniform time step solutions. Recently, Zhu et al. [22] applied the MTS method into train-track-bridge dynamic analysis by implementing fine and coarse time steps to solve the train-track subsystem and bridge subsystem. Besides, scholars developed iterative algorithms to obtain TTS dynamic solutions [9,23].
Summarily, pioneering work presented above intend to reduce the computational cost in the train-track-substructure dynamic analysis. In the work of Zhai et al. [5], the explicit-implicit hybrid integration method is proposed and proved to be accurate and efficient for the solution of train-track-bridge dynamic interaction. Following this methodology, MTS solution method is applied as a mI-nE type, m and n can be arbitrary real numbers, and the explicit Zhai method is particularly introduced to efficiently solve substructure responses that occupies the most part of the degrees of freedom (DOFs) in the TTS system, so as to maximizing the computational efficiency.
The following part of this paper is organized as follows. In Sect. 2, a brief introduction of the establishment of the TTS interaction model is presented. The solution method for TTS interaction systems is developed in Sect. 3. In Sect. 4, examples are provided to validate the proposed method. Conclusions are drawn in the last section.

Introduction of train-track-substructure interaction model
Train-track-substructure dynamic interaction model is developed from the vehicle-track coupled dynamics [6,24]. It treats the train, track and various substructures (bridge, subgrade, etc.) as an entire system. Generally, the train-track interaction system, the substructure system and their coupled solution can be, respectively, depicted as below.

Modeling of the train-track interaction
In the train-track interaction model subject to the substructure, as shown in Fig. 1, multi-rigid-body dynamics and finite element theory are applied. The train is modeled as a series of vehicles. Each vehicle has one car body, two bogie frames and four wheelsets, which are assumed as rigid bodies connected by primary and secondary suspension systems. Each body has six degrees of freedom (DOFs), i.e., displacements in the longitudinal x , lateral y , vertical z directions, and angles of yaw , pitch and roll ; thus each vehicle has 42 DOFs.
The track system mainly includes two types: one is ballasted track and the other is the ballastless track. Ballasted track system consists of the rail, sleeper, track bed and spring-dashpot elements, i.e., the rail and the sleeper are all assumed as spatial Bernoulli-Euler beams characterizing the lateral and vertical motion for the rail, and longitudinal and vertical displacement for the sleeper. The longitudinal motion of the rail and the lateral motion of the sleeper are assumed as a bar motion. As to the longitudinal and lateral motion of the rail and sleeper, bar elements are introduced, namely each node of the rail and sleeper beam has 6 DOFs. The rail pads between the rail and the sleeper are modeled as spring-dashpot elements. The track beds are treated as mass elements with 3 DOFs, i.e., longitudinal-, lateral-and 22 L. Xu et al.

3
vertical-displacements, and considering the occlusal actions between ballasts represented by stiffness and damping coefficients.
While for the ballastless track system, it mainly consists of the rail, track slab, filling layer and concrete basement, and interlayer spring-dashpot elements. The rail is modeled as a Bernoulli-Euler beam, the track slab and concrete basement are modeled as an assemblage of thin-plate elements. The interaction for the rail-track slab subsystem, and the track slab-basement subsystem is, respectively, modeled as discrete and facial supported spring-dashpot elements.
To couple the train and the track, wheel-rail spatially coupled dynamics model developed in [24] has been applied to clarify the wheel-rail contact geometries and parameters, and then wheel-rail coupling matrices derived by energy variation principle are formulated (see Ref. [25]). In this way, the train-track interaction model is established. More details have been presented in [9], which can be referred to.

Modeling of the substructure
Generally, the substructures with regular geometry morphology can be conveniently modeled by self-compiled program [7][8][9]. However, when it involves complex structural geometries, element types and local defects, the commercial finite element software seems to be more practical in engineering application. In this work, ABAQUS @ is applied to model the substructures, and then the mass M ss , damping C ss and stiffness K ss matrices are exported, also with the exportation of the nodal coordinates used to localize the nodal DOFs imposed forces from the track.

Time-dependent coupling procedure between the train-track system and the substructure system
Unlike conventional models [24,25], the moving length of the train is equal to that of the track in the train-track-substructure coupled dynamic analysis. In this work, the track length is modeled as a minimum value of L T + 2l t no matter how long the train moves forward, where L T and l t denote, respectively, the train and the boundary length. The train-track system and the substructure system are intrinsically coupled by their interaction forces. When there is no contact area between the track and substructure, such as scenarios A and C shown in Fig. 2, the track-substructure interaction force vector can be obtained by where k ts and c ts denote the track-substructure interaction stiffness and damping coefficient respectively; the subscripts t and s denote the track and substructure system, respectively; the subscript denotes a representation of the forces in the longitudinal-, lateral-and vertical-directions; X t, and Ẋ t, denote the displacement and velocity of the track bed with respect to various directions; X t, and ̇X t, denote the displacement and velocity of the track bed element; N t denotes the shape function, and N t = 1 in this work.
(1) F ts, = k ts, X t, + c ts, Ẋ t, , When there is contact area between the track and substructure, such as scenario B shown in Fig. 2, the track-substructure interaction force can be obtained by where X s and Ẋ s denote the nodal displacement and velocity of the substructure with respect to various directions, and for non-contacting area in scenario B, Eq. (1) is applied.
In the time-dependent coupling procedure, another key except for the calculation of the interaction forces is accurately locating the track and substructure DOFs at the contact area on which the interaction forces are exerted. Besides, there exists DOFs mapping between the global coordinate The detailed coupling method has been presented in [26], here not presented for brevity.

Solution for train-track-substructure dynamic interaction
The explicit-implicit hybrid integration method originated from Ref. [5], in which the train-track system is solved by explicit Zhai method, and the bridge structures are solved by well-known implicit Newmark-β method. Regarding that this work put an emphasis on the structural dynamics instead of the very high frequency wheel-rail vibrations, the responses of upper structures built by MATLAB® with relatively low DOFs are obtained by implicit integral scheme, and the responses of the sub-structures built by ABAQUS® with very large DOFs are obtained by the explicit Zhai method by treating the mass matrix of the substructure as diagonal form. Furthermore, different time step sizes are applied for the upper-and substructures to improve the computational efficiency.
It should be noted that wheel-rail interaction generally shows high frequency vibrations, especially with consideration of rail short-wavelength roughness, and applying an explicit integral scheme is more efficient [6]. However, if mainly concerning the structural dynamic performance such as the track structures and substructures, the vibration frequency analyzed can be decreased properly. Moreover, the beam element with the consistent mass matrix is applied in this model and is coupled to the train model in matrix formulations. With above consideration, the wheel-rail coupled interaction is solved by implicit integration method.

Implicit-explicit hybrid-integral schemes
To achieve the solution for such a large and complex train-track-substructure dynamic system, computational efficiency becomes a concern apart from the accuracy. Hence, an explicit-implicit hybrid-integration scheme combining the Park method [27] and Zhai method [19] is developed to efficiently obtain the train-track solution and substructure solution, respectively. Except for the wheel-rail interaction, the train-track system also undergoes time-dependent reaction forces from the supporting layer under the track beds.
With acquisition of the wheel-rail forces [28], the initial force vector for the train-track system at the time step n can be obtained as where t denotes an assemblage of track layer elements contacting the substructure surface.
To obtain the train-track dynamic responses at time step n + 1, the following integral schemes are applied by introducing Park integration method: Step 1 Calculate the equivalent stiffness matrix by where Δt denotes the time step size.
Step 2 Calculate the equivalent force vector by , where the superscripts i , i − 1 and i − 2 denote the three adjacent time steps before the (n-1)th time step.  Step 3 The displacement, velocity and acceleration response vector can be obtained by Step 4 Update the train-track system responses at adjacent three steps: The implicit Park method can be represented as a function form: As to the substructural system, its force vector can be assembled as where s denotes an assemblage of substructural nodal DOFs contacting to the upper tracks.
By introducing explicit Zhai's integration method, the displacement and velocity response vectors of the substructural system can be obtained by [19] where = = 0 n = 1 0.5 n > 1 .
The acceleration responses can be further derived as The explicit Zhai method can be represented as a function form: Through Eqs. (3)-(11), it can be cognized that all subsystem responses of the train, track and substructure can be solved time-dependently by updating the wheel-rail forces and track-substructure interaction forces by the implicitexplicit hybrid-integral schemes.

Multi-time-step solution procedures
In the multi-time-step solution procedures, iteration procedures are not of necessity within the linear elastic dynamics. However, the time step size required to guarantee numerical solution stability possesses differences with respect to different integral schemes. With this consideration, multi-time-step solution procedures are proposed to improve the computational efficiency.
Set a time axis to label the time step vector to characterize train-track system solution T t and substructural system solution T s , denoted by T t = 0, Δt t , 2Δt t , ⋯ and T s = 0, Δt s , 2Δt s , ⋯ , respectively, Δt t and Δt s are, respectively, the time step size for the train-track implicit-integral solution and substructural explicit-integral solution.
Two steps below can be followed to achieve multi-timestep solution: Step 1 Set the step number N t = 1 and N s = 1 for the upper-structure (abbreviated as 'Us') and the sub-structure (abbreviated as 'Ss'). Select a start moment as the initial time point at which the dynamic equations of the Us and the Ss are, respectively, solved by Park method and Zhai method, as illustrated in Fig. 3. The time step vectors for the Us and Ss are extended from T t = (0) → T t = 0, Δt t and T s = (0) → T s = 0, Δt s , and integral step numbers for the Us and Ss are updated as N t = N t + 1 , and N s = N s + 1.
Step 2 When N t ≥ 2 and N s ≥ 2 , the relative position of the integral time points for the Us system and Ss system at the time axis is judged and solution is obtained as follows.
a T t N t = T s N s , the Us system and the Ss system will be solved simultaneously by implicit-explicit hybrid integration method. 1. The interaction forces between the Us system and the Ss system for obtaining the next step solution can be calculated by With acquisition of the interaction forces, the implicitexplicit hybrid-integral schemes presented in Sect. 3.1 can be used to obtain all sub-system responses.
The Ss responses are chosen from the last step solution, i.e., X N t t and Ẋ N t t . 2. The Us-Ss interaction forces can be consequently obtained by 3. Solving the dynamic equations of the Ss systems by Zhai method as 4. Update the step number and time step vector for Ss system N s = N s + 1,  , The Us system responses are chosen from the last step solution, i.e., X N s s and Ẋ N s s . 2. The Us-Ss interaction forces can be consequently obtained by The system gravitational forces should also be considered for the Us system as shown in Eq. (2).

Solving the dynamic equations of the Us systems by
Park method as 4. Update the step number and time step vector for Us In the train-track-substructure dynamic interaction, the train is constantly moving on the tracks in real scene. However, it is noted that the dynamic integral processes for the train-track system, i.e., the Us system, can be only implemented in conditions satisfying T t N t ≤ T s N s , or the train should be static without moving forward.
Summarily the implicit-explicit integration and multitime-step solution can be embedded in the TTS dynamic interaction as shown in Fig. 4.

Numerical examples
Three examples are presented to show the effectiveness, applicability and practical application of this work in evaluating the large-scale train-track-substructure interaction system, even with complex soil foundations.

Example 1: effectiveness of the proposed implicit-explicit integration method
The train and track parameters are listed in Appendix A, respectively. The substructure, as shown in Fig. 5, is treated as an isotropic soil sub-system with parameters listed in Appendix B. The soil substructure is modeled as a lumped mass system. The total number of DOFs for the train, the track and the soil substructure are, respectively, 126, 12,438 and 72,414. The train consists of three identical vehicles with a constant running velocity of 200 km/h. Measured track irregularities on a ballasted line are considered as the excitation. The geometry configuration for the model is shown in Fig. 6.
To validate the accuracy and efficiency of the proposed implicit-explicit integration method, and to show the robustness of the method in arbitrarily choosing the time step sizes once satisfying integration stability for Us (train-track system) and Ss (substructure system), several cases are implemented, as illustrated in Table 1. From the time steps set in Table 1, it is known that the time steps of the Us and Ss can be non-integer multiple relations, namely an mI-nE integration mode. Besides, the multi-time-step solution procedures can be also applied to the unified implicit solutions, namely the Us and Ss are solved by Park method only, namely an mI-nI integration mode, and the results are treated as comparisons. Figure 7 shows the comparisons on the displacement and acceleration of the soil foundation with respect to different solution methods and time step sizes. It can be obviously seen from Fig. 7 that stable solution for the displacements and accelerations can be obtained by different solution methods though time step sizes are varied greatly. The results obtained by the mI-nE solution method and the mI-nI solution method coincide significantly well. Moreover, it is noted that the responses show slight deviations due to the difference of excitation frequencies caused by the variation of time step sizes. Generally when the time step size is varied from 1.5 × 10 -4 s to 5.5 × 10 -4 s, the response difference is smaller than 0.25%.
Apart from the comparisons on soil vibrations, the comparisons on car body accelerations and wheel-rail forces can be also performed both from the lateral vibration and vertical vibration. As illustrated in Fig. 8, the car body lateral and vertical acceleration show approachable results with respect to different solution methods. The relative deviation is generally smaller than 0.1%. Besides, Fig. 9 further presents comparisons on wheel-rail forces, from which it can be clearly seen that the time step sizes show remarkable influence on wheel-rail forces, especially with participation of short-wavelength irregularity excitations at high frequencies. Obviously local wheel-rail forces at time step size of 1.5 × 10 -4 s represented by the red solid and dashed lines are larger than those by time step sizes 3.5 × 10 -4 s and 5.5 × 10 -4 s, as indicated by the power spectral density (PSD) of wheel-rail forces in Fig. 10.
To show the computational efficiency of the proposed mI-nE-MTS method, the time spent by it and the mI-nI-MTS method is listed in Table 2. It is shown in Table 2 that the computational time is significantly reduced by the usage of explicit method, moreover, the computational efficiency is also improved by adopting the MTS strategy.

Example 2: applicability of the explicit method in solving substructural vibration
In above studies, the mass matrix of the soil system is represented as lumped type to be solved by the explicit Zhai method. It is known that the consistent mass matrix is also widely used in modeling railway structures [5], and solved by implicit algorithms such as Park method, Newmark-β method, etc. It is therefore of necessity to perform comparison analysis between the dynamic responses of the lumped mass soil system and the consistent mass one. In this example, the soil system is modeled by MAT-LAB® program with 26,568 nodes, i.e., 79,704 DOFs, the related parameters are shown in Appendix B, and the lumped and consistent mass matrices of a solid element are shown in Appendix C. Setting time step size vector as Δt = [0.1, 0.2,     It is shown in Figs. 13 and 14 that the soil acceleration converges faster by applying implicit Park method at time step size 2 × 10 -4 s, and from viewpoint of soil responses at frequency-domain, it is known that the divergences between the implicit Park solution and explicit Zhai solution mainly appear at high frequencies larger than 260 and 343 Hz with respect to the lateral and vertical accelerations. The high frequency response is not dominant to the soil system, and when the time step size reaches 0.2 × 10 -4 s, the high frequency responses are significantly dissipated in the solution by Zhai method. Moreover, the calculation time required for Zhai method is, respectively, 2.57 and 3.46 h for time step sizes 2 × 10 -4 and 0.2 × 10 -4 s, but for Park method, the time spent is significantly increased to 100.49 h with time step  size of 2 × 10 -4 s. That is to say, the calculation efficiency is increased by 50 times by applying Zhai method comparing to Park method. Though it can be observed that there possess obvious differences for the soil acceleration between Zhai method and Park method at frequencies higher than 260 Hz for lateral acceleration and 343 Hz for vertical acceleration. However, the soil vibration at such high frequencies is generally small, and its influence on the dynamic evaluation is negligible in practice.

Example 3: train-track-bridge interaction with complex tunnel-soil-pile foundations
As a real scenario shown in Fig. 15, a shield tunneling construction is conducted to underpass a high-speed railway viaduct. Given the influence of shielding tunnels on foundation stiffness and a series of problems such as subgrade settlement, it is of necessity to evaluate the running safety of a train on the track-bridge system with complex tunnel-soil-pile foundations. To achieve this goal, the complex foundations including soil, tunnel, pile, and pile caps are established by the ABAQUS® software, as illustrated in Fig. 16. There are totally 70, 3680 DOFs. By applying the implicit-explicit integration and multi-time-step method, and coupling the train-track-bridge interaction system and the complex foundation system in the MATLAB platform by procedures presented in Sect. 2.3, where the stiffness, damping and mass matrices of the tunnel-soil-pile foundation are exported from the ABAQUS®.
To validate the accuracy of the solution by implicit-explicit multi-time-step solution, comparisons are made between this model and the ABAQUS model. In the calculation of the system responses by ABAQUS, the external excitations, recorded from the interaction forces between the bridge girder and the pier, are time-dependently loaded on the surficial nodes of the Pier elements. Figure 17 shows the comparisons of pier displacement and acceleration, from which it can be observed that the maximum displacements of this model and ABAQUS model are, respectively, 8.487 and 8.317 μm, and the maximum accelerations are, respectively, 0.02937 and 0.02784 m/s 2 . Obviously, this model can obtain approachable results comparing to those by ABAQUS model.
The response deviations lie in the difference of the time step sizes and solution methods, which directly influence the dynamic solutions. The time step size used in this model is 2 × 10 -5 s by explicit Zhai method and the time step size used in the ABAQUS model is 10 -3 s by ABAQUS/Explicit Quasi-static solution. However, CPU time consumed by this model and ABAQUS model is, respectively, 32,294 and 57,184 s; namely this model increases the computational efficiency by 43.53%, not to mention the fact that the ABAQUS model is relatively difficult to model and analyze the train-track-bridge dynamic interaction analysis considering complex wheel-rail contacts.
To show the influence of the shielding tunnel on the train running performance on the track-bridge system, conditions with and without underground shielding tunnel are considered as C 1 and C 2, respectively. Figure 18 shows the comparisons on bridge pier and girder under various foundation conditions. It is shown in Fig. 18 that  the dynamic displacement of the bridge pier is obviously enlarged in C 1 condition, which is triggered by the softening of the sub-soil foundation by shielding the tunnel. The maximum displacements for C 1 and C 2 conditions are, respectively, 49.5 and 11.2 μm, and the maximum accelerations for C 1 and C 2 conditions are, respectively, 0.12 and 0.102 m/s 2 .
As to the responses of the bridge girder illustrated in Fig. 18c and d, it is seen that the maximum displacement of C 1 is also larger than that at C 2 condition, but with slight deviations smaller than 8 × 10 -3 mm. In addition, the train running comfort and safety can be also evaluated by the proposed method. As illustrated in Fig. 19, the car body vertical acceleration and wheel unloading rate with respect to different foundation conditions are presented. The maximum differences against the car body acceleration and wheel unloading rate are, respectively, 1.4 × 10 -4 and 1.35 × 10 -3 s, indicating that the shielding tunnel exerts a slight influence on the train running performance on comfort and safety in this case study, where the effects of the periodicity of the tunnel segment are not considered.

Conclusions
In this work, an implicit-explicit integration and multi-timestep method was proposed. The implicit Park method and the explicit Zhai method were introduced into multi-time-step solution procedures, in which the time step size for the implicit integration and explicit integration can be arbitrarily chosen once satisfying the integration stability and maximum frequency interested. The examples demonstrate the practicality of the explicit Zhai method in solving large-scale substructure dynamics by satisfying the lumped mass matrix assumption. In specific conditions, it increases the computational efficiency by 50 times than the Park method. Certainly, the maximum natural frequency of vibration should also be concerned in applying the explicit integration method.
Obviously, the implicit or explicit method for solving the train-track dynamic equations can be also replaced by methods of Wilson-θ, Newmark-β, etc., because they are applicable to consistent mass system instead of lumped mass. In other words, the integral schemes applied to the train-track subsystem and substructure subsystem are alternative.
Acknowledgements Appreciation is given to Hubing Liu and Zixu Zhu, the postgraduate students under substantive guidance of Lei Xu, for their work in building and exporting the ABAQUS ® model to be capable of integrating with the self-compiled program in MATLAB ® . This work was supported by the National Natural Science Foundation of China (Grant Nos. 52008404, U1934217 and 11790283); Science and Technology Research and Development Program Project of China Railway Group Limited (Major Special Project, No. 2020-Special-02); the National Natural Science Foundation of Hunan Province (Grant No. 2021JJ30850).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.  x = x 1 N 1 + x 2 N 2 + x 3 N 3 + x 4 N 4 + x 5 N 5 + x 6 N 6 + x 7 N 7 + x 8 N 8 y = y 1 N 1 + y 2 N 2 + y 3 N 3 + y 4 N 4 + y 5 N 5 + y 6 N 6 + y 7 N 7 + y 8 N 8 z = z 1 N 1 + z 2 N 2 + z 3 N 3 + z 4 N 4 + z 5 N 5 + z 6 N 6 + z 7 N 7 + z 8 N 8 ,