Train–track coupled dynamics analysis: system spatial variation on geometry, physics and mechanics

This paper aims to clarify the influence of system spatial variability on train–track interaction from perspectives of stochastic analysis and statistics. Considering the spatial randomness of system properties in geometry, physics and mechanics, the primary work is therefore simulating the uncertainties realistically, representatively and efficiently. With regard to the track irregularity simulation, a model is newly developed to obtain random sample sets of track irregularities by transforming its power spectral density function into the equivalent track quality index for representation based on the discrete Parseval theorem, where the correlation between various types of track irregularities is accounted for. To statistically clarify the uncertainty of track properties in physics and mechanics in space, a model combining discrete element method and finite element method is developed to obtain the spatially varied track parametric characteristics, e.g. track stiffness and density, through which the highly expensive experiments in situ can be avoided. Finally a train–track stochastic analysis model is formulated by integrating the system uncertainties into the dynamics model. Numerical examples have validated the accuracy and efficiency of this model and illustrated the effects of system spatial variability on train–track vibrations comprehensively.


Introduction
Originated from the uncertainty of manufacturing error, material fatigue and damage, complex excitations, environmental effects, etc., the stochasticity of the train-track interaction becomes an essential characteristic for this dynamics system. Generally, the random evolution of system properties in geometry, physics and mechanics is aroused by the dynamic interaction between train and track in space. The system spatial variability will inversely influence the system vibration, accelerating the system property evolution.
Unlike most of the other dynamic systems, train-track interaction takes place in a longitudinally large field with viscoelasticity, nonlinearity and high-dimensional degrees of freedom (DOF). Consequently, the stochastic characteristics of geometric, physical and mechanical parameters of the system are actually scattered in a wide range. With consideration of the system property evolution, it is anticipated that train-track interaction is both a random process and possessing abundant random information. See for instance, the rail profile irregularities, generally denoted as track irregularities, exist inevitably and randomly along the whole railway line which might be hundreds or thousands of kilometres long, not to mention the track irregularities that vary constantly with train moving loads, track settlement, structural deformation, etc.
Train-track interaction is intrinsically an assemblage of system parameters with various mechanical, physical and geometrical properties, which can be effectively coupled and characterized by dynamics methodologies. Aiming at clarifying system stochastic behaviours, the spatial variability of system geometric and parametric excitations will be considered with an emphasis in the modelling construction and the dynamics analysis.

Spatial variation of system geometric property
It is well known that the geometric deformation inevitably exists in dynamical systems to accommodate forces applied to the structures or systems. In the train-track interaction system, one of the most important and unavoidable excitation is definitively the track irregularities and also the system dynamic behaviour is indeed to be characterized on spatially varied track portions. The track irregularity has long been observed, measured and defined by worldwide researchers and institutions [1,2]; generally, power spectral density (PSD) function can be used to statistically represent the track irregularities at frequency-amplitude field by definitive formula, such as the US railway classes 1-6 and German low-and high-speed spectrums. Using time-frequency transformation method [3], time-domain track irregularities, loaded as an input of train-track interaction, can be equivalently obtained from the PSD function.
Commonly track irregularities measured from an entire railroad belong to the big data category, and consequently, it is time-consuming and less efficient if loading all the time-domain track irregularity signals into the dynamical computations. Besides it is unreliable to reveal all system stochastic behaviours if only the statistical average PSD function or its relevant time-domain track irregularity is used as the excitation. Highlighting this concern, pertinent work on sampling track irregularity random sets from massive data becomes an interesting topic. Perrin et al. [4] were devoted to develop a stochastic model for track irregularities with experimental validations. In this work, random field theories, such as expansions of Karhunen-Loève and Polynomial Chaos, are introduced to account for the statistical variability and dependency of track irregularities along track abscissa. In recent years, Xu et al. [3] strode forward to propose a practical track irregularity probabilistic model by utilizing the shape similarity of track irregularity PSD function in a specific railroad. It is proved that this approach can reach approachable accuracy as the model in [4] but with higher efficiency.

Spatial variation of system properties in physics and mechanics
Except for track irregularities, parametric excitation, triggered by the spatial variability of system properties in physics and mechanics, is another objectively existing excitation type. This work has early been studied by Náprstek and Frýba [5]. In [5], the railway track and its substructure are simplified as an infinite beam resting on Winkler foundation. The stiffness of which is randomly generated by the spectral density of diffuse type with an exponential correlation. Later, Fröhling [6] used a mathematical model to discuss the effects of spatially varying track stiffness on vehicle-track dynamic behaviours. Besides, Oscarsson and Nielsen [7][8][9] made strides in investigating the influence of the scatter in track properties with a numerical simulation model for train-track vertical interaction; mostly important, full-scale measurements in the field and laboratory measurements were carried out to show the unevenness of track properties such as rail pad stiffness, ballast stiffness and sleeper spacing. Also, Andersen and Nielsen [10] dealt with a problem by perturbation analysis, where a one-DOF vehicle moves along a simple track beam with stochastically varying support stiffness. Moreover, Wu and Thompson [11] treated the sleeper spacing and ballast stiffness as random values at each support point. It is found that the point receptance and the vibration decay rate are changed by the variability of track properties.
In the last decade, researchers began to notice the spatial variability of track properties and its influence on track deterioration, settlement and train-track dynamic behaviours, etc. For example, Dahlberg [12] assumed that the track stiffness irregularities mainly originate from two aspects: one is the track superstructure, e.g. rails, rail pads and ballast; the other is the substructure, e.g. foundation and subgrade soil. In his work, it is demonstrated that the variation of wheel-rail interaction may be considerably reduced by an optimized design of the track stiffness variation such as the use of grouting or under-sleeper pads. More specifically, some researchers as shown in [13] paid special attentions on track settlement, track transition, switches, turnout and rail joints, etc. The track stiffness irregularities of these track portions are born to exist and mainly analysed by deterministic dynamic analysis methods. The relevant work is rather abundant, see the literature review by Sañudo et al. [13].

Characterization of train-track interaction
Apart from random simulations of system geometric, physical and mechanical properties, it is known that another key work is to integrate these uncertainties into a dynamics system with satisfaction of mechanics principle, that is, the modelling of train-track interaction with uncertainties.
Generally track structures are modelled by beam, thinplate and mass elements, and the interaction between track components is depicted by linear/nonlinear spring-dashpot elements. In general, the vehicle is modelled as a multirigid-body system with two-stage suspensions. The methodologies to build the vehicle and the track system can be generally classified as mode superposition method and finite element method. The difficulty widely known is the characterization of the wheel-rail interaction in a threedimensional (3D) space. Regarding the differences in wheel-rail contact/creepage descriptions, i.e. rigid contact [7][8][9][14][15][16][17], elastic contact [3,[18][19][20] and elastic-plastic contact [21][22][23], the complexity, accuracy and efficiency of the model are significantly different.
To railway dynamics subject to large-scale stochastic problems, the elastic-plastic model is accurate but inefficient to train-track dynamics problems at a macro-level, where large random samples may be accounted for and the random variable vector is of high dimensionality in space. Instead the wheel-rail rigid and elastic contact model is far more efficient and generally accurate enough to clarify the external loads of rail substructures to explore the internal strain-stress variation in engineering practices.

Outline of this work
It can be observed from the state-of-the-art work presented above that the system geometric uncertainty and physicsmechanics uncertainty have been studied as a hot topic but generally in an independent way in the train-track (a)  Fig. 1 Train-track interaction model (a side view; b end view) (the symbols x, y, z, w, b and h denote the longitudinal, the lateral, the vertical, the yaw, the pitching and the rolling motion of the bodies, respectively; the subscripts 'c', 'b' and 'w' denote the car body, the bogie frame and the wheelset, respectively) Rail. Eng. Science (2020) 28(1): [36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53] interaction, and sometimes without quantitatively qualifying the uncertainty effects from statistics.
In this work, an integrated research is presented with a goal of constructing a train-track coupled stochastic analysis system, where system uncertainties in geometry, mechanics and physics as the system excitation in a realistic and representative way have been wholly accounted for. The organization of the following paper is as follows: • In Sect. 2, the modelling method for the train-ballasted track interaction is presented with brevity. • In Sect. 3, the spatial uncertainty of system excitations on geometry, physics and mechanics is illustrated and the quantification methods for these random excitations are elaborated. Besides the framework is formed by integrating the presented works on dynamics model construction and uncertain parameter simulation. • In Sect. 4, numerical examples are conducted to validate the proposed method and to survey the influence of spatial uncertainty of system geometry and physics-mechanics property on train-track interaction. • Finally in Sect. 5, some concluding remarks are presented.

Dynamics model for train-track interaction
Based on the theory of vehicle-track coupled dynamics [18,19], the coupled matrices, as representation of the dynamic equations of motion for the train-track interaction shown in Fig. 1, can be established as where M, C and K denote the mass, damping and stiffness matrices, respectively; the subscript 'V' and 'T' indicate quantities for the systems of train and track, respectively; 'VV' and 'TT' indicate matrices for the train and the track, and 'VT' and 'TV' indicate matrices for the interaction between the train and the tracks; € X, _ X and X denote the acceleration, velocity and displacement vectors, respectively; F denotes the loading vector; M VT and M TV are non-zero matrices for wheel-rail rigid contacts.

Construction of matrices for the train system
A train includes a group of vehicles, including motor cars and trailers that are connected by the coupler and draft gear system. The vehicle is modelled as a multi-rigid-body system consisting of a car body, two bogie frames and four wheelsets. The system components are mechanically connected by two-stage suspension systems.
The detail method for constructing train matrices M VV , C VV , K VV and F V has been presented in [24] for references.

Construction of matrices for the track system
The track is modelled as a commonly used ballasted track system by FEM. The rails are models as Bernoulli-Euler beam, the nodal displacements and rotations towards X-, Yand Z-axes are considered. The sleeper and the track bed are modelled as a rigid body and a mass, respectively. The linear displacements of the sleeper along Y-and Z-axes and angular displacement around X-axis are considered. The vertical displacement of the track bed is accounted for, including the shearing effect between track beds.
The detail method for the establishment of the track matrices, i.e. M TT , C TT and K TT , has been presented in [25] by finite elemental formulations.

Coupling method for train-track interaction
The following work is coupling the vehicle subsystem and the track subsystem through the wheel-rail interaction as shown in Fig. 2. The wheel-rail lateral interaction induced by the tangential creepages (lateral, longitudinal and spin), closely correlating to the wheel-rail relative velocity and creep coefficient, has been described by the fundamental work of Kalker [26], Vermeulen and Johnson [27], etc. While for the wheel-rail vertical interaction induced by the normal contact (or compression), different assumptions are made, such as the wheel-rail rigid contact [24], no compression between the wheel-rail normal contact and the wheel-rail elastic contact by assuming non-adhesive nonlinear normal contact for two spherical solids [18].
As a supplement to the vehicle-track coupled dynamics method in [18] and following the convenience of modelling of finite element system, the wheel-rail interaction is characterized by matrix coupling formulations in energy principle instead of deriving wheel-rail forces explicitly.

Wheel-rail vertical coupling matrices
Following the Hertzian contact theory, the wheel-rail vertical force F wr can be expressed by where G is the wheel-rail contact constant (m/N 2/3 ); Dz is the wheel-rail elastic compression (m). From Eq.
(2), the wheel-rail equivalent contact stiffness k wr,z can be obtained: Accordingly, the wheel-rail coupling matrix K wr can be obtained as where the symbol U j;g denotes the displacement vector corresponding to the shape function; i denotes the ith vehicle; j denotes the jth wheelset and g ¼ 1; 2 denote, respectively, the left side and right side of a wheelset; d lr denotes half distance between the left-and right-side wheel-rail contact points; the subscripts 'w' and 'r' denote the wheelset and the rail, respectively; z and u denote the vertical motion and the roll motion, respectively; and v denotes the virtual coordinate of track irregularity.

Wheel-rail lateral coupling matrices
Based on the vehicle-track coupled dynamics theory [18,19], the wheel-rail lateral coupling matrices have been derived by energy variation principle in [25], which can be introduced to this model accordingly.
Till now the methods for coupling the train and the tracks have been illustrated with brevity. The train and the tracks are effectively united as an entire system by coupled matrix formulations. Through numerical validations, it had been proved that the dynamically coupled matrices of Eq. (1) can be solved by numerical integral schemes with high stability and accuracy even at a relatively larger time step size; besides no iterative procedures are required in the numerical integration, indicating that one can simultaneously obtain the dynamic responses of the train and the track responses at each time step.
3 Spatial uncertainty and quantification of system excitation on geometries, physics and mechanics In course of a train-track interaction event, system excitations (H), including track geometries and system parameters of physics and mechanics, show randomness in the spatial axle, which can be assembled by where q denotes the sampling number: q ¼ 1; 2; . . .; g, and g is the total number; the subscript 1 denotes the track irregularity vector, and n denotes the system physics and mechanics parametric vector.

Generation of correlated pseudo-random variables following arbitrary probability distribution function
In Eq. (5), it can be known that there is a high-dimensional random vector consisting of various random variables following arbitrary probability distributions and probably possessing correlations between random variables. For generating random variables with m-dimensional correlations, a linear and nonlinear two-step transformation method [28,29], where the sampling sequences of multidimensional correlated random variables with specified edge distribution and correlation coefficient can be obtained, is applied in this present study. The detail method is presented in ''Appendix 1''.

Random simulation model for track irregularities
Xu et al. [3] have previously proposed a model for obtaining highly representable random samples from massive track irregularity data by probabilistic methodologies, but it ignores the spatial correlation of track irregularities of different types; moreover, the track quality representation of track irregularity spectral density function has not been illustrated. For solving above issues, a new model for random simulation of track irregularities is extensively presented on the foundation of spectral representation method. In a 3D space, the left-and right-side rails both have vertical and lateral irregularities, and the track irregularity vector can be therefore assembled by H q;1 ¼ ðk l ðsÞ; m l ðsÞ; k r ðsÞ; m r ðsÞÞ; s 2 ð0; where c 1 , c 2 , c 3 and c 4 denote the vertical profile, alignment, cross-level and gauge irregularities, respectively; k and m denote the lateral and vertical track irregularities, respectively; the subscript 'l' and 'r' denote the left-and right-side rail, respectively; s is the longitudinal abscissa of the track, S is the total length.
The PSD function of track irregularity vector can be expressed by where Ç ðÁÞ denotes the PSD estimator.

Spectral correlation of track irregularities
Divide the total length S into N segments, i.e. S ¼ ðn À 1Þ S=N; n S=N ð , n ¼ 1; 2; . . .; N. There will be N PSD functions for each track irregularity type.
According to discrete Parseval theorem [30], one can derive the following expression as where S w ðÁÞ denotes the discrete Fourier transform (DFT) of track irregularities; Ds and Dx denote the time and frequency interval, respectively; w denotes the discrete time-domain track irregularity signals; t denotes the total discrete number; P denotes the PSD vector. Equation (8) has exposed a close relation between time signals and its PSD. An integration indicator U w ðÁÞ used to determine the PSD representing certain track geometric status can be therefore defined by where b is the PSD sample; M denotes the total number of frequency points.
With the definition of correlation coefficient, the spectral correlation of track irregularities can be calculated by ; ði; jÞ 2 ½c l ; m l ; c r ; m r ; i 6 ¼ j; where C i;j denotes correlation coefficient; U i denotes the integration indicator, CovðÁÞ denotes the covariance function, and VarðÁÞ denotes the variance function Based on Eq. (10), the correlation coefficient matrix (CCM) for the track irregularity integration indicator U i can be obtained as where the subscript W denotes the distribution type of the track irregularity integration indicator; 'C i;j ' with subscript i; j¼ ðk l ; k r ; m l ; m r Þ denotes the correlation coefficient between the ith and the jth irregularity type. When i ¼ j, C i;j ¼C j;i ¼ 1.

Random simulation of PSD with specific correlations
From Eq. (9), one knows that track irregularity PSD can be approximately defined by an integration indicator. Thus the simulation of PSD is equivalent to the simulation of U w ðÁÞ. Generally, the statistical PDF of U w ðÁÞ differs from a normal distribution. The method presented in Sect. 3.1 can be therefore applied to choose the PSD samples U w ðÁÞ.

Random simulation of track irregularities
With the work above, one can naturally obtain the random variables satisfying specific correlations. The algorithm is listed concisely as below.
In Algorithm 1, E is an indicator for judging the deviation between the simulations and the standards; e 0 is the original value assigned to E; e is a tiny value, e\e 0 ; Z n9m is a real number matrix with order n 9 m for the variable z following normal distribution; } denotes the random simulation operator for standard normal distribution; R W!N denotes the CCMs correlation coefficient matrices transformed from arbitrary distribution (W) to normal distribution (N); CholðÁÞ denotes the Cholesky decomposition and CorðÁÞ denotes the operator of correlation coefficient.
Because the track irregularities k l;r ðsÞ and m l;r ðsÞ can be realistically measured, the PDF of the PSD indicator U w ðbÞ, b ¼ 1; 2; . . .; N, can therefore be obtained statistically and discretely, represented as f U ðbjwÞ.
With acquisition of f U ðbjwÞ, the random simulation of the PSD indicator can be calculated by where F denotes the cumulative distribution function of the PSD indicator. Obviously U w ðbÞ corresponds to the PSD function P w ðxÞ, as shown in Eq. (9). In train-track interactions, track irregularities, as the system excitation, must be equivalently transformed into time-domain processes. The algorithm is presented as below. In Algorithm 2, f y is the frequency domain of P w ; M il is the length of the time-domain track irregularities; N r is the sampling number; f y is the frequency range; f l and f u denote the lower-and upper-limit frequency; conjðÁÞ, interp1ðÁÞ and ifftðÁÞ denote the conjugate function, 1D interpolating function and inverse Fourier transformation function, respectively; U follows uniform distribution at (0, 1) with length of L p ; and u is the phase distribution.
Till now, the methods for random simulation of track irregularities have been constructed entirely. In this method, the random variables are no longer limited to normal distribution. Instead, arbitrary distributions with/ without correlations can be simulated efficiently and accurately. Set 1 km as a unit length of the track irregularity PSD. The CCM of track irregularities can be calculated out. Through transformations, the equivalent CCM of the normal distribution can be obtained consequently.

Model validations
The detailed values of CCM are presented as follows: Based on the methods presented above, the random samples of PSD indicators can be selected by Algorithm 1, as shown in Fig. 4, where / v l , / v r and / c l denote the vertical irregularities at the left and right sides and the lateral irregularities at the right side, respectively. It can be clearly observed from Fig. 4 that the vertical irregularities, respectively, at the left and right sides of the rails are highly correlated; however, the rail irregularities at different directions (lateral and vertical) are correlated weakly. Besides, the correlation coefficients C m l ;m r and C m l ;c l are, respectively, with values of 0.914 and 0.222, which are close to the reference values 0.918 and 0.224, respectively.
The PSD functions can certainly be determined using the integration indicators demonstrated in Fig. 4. Then Algorithm 2 is applied to obtain time-domain track irregularities that are randomly transformed from PSD functions. To validate the accuracy and efficiency of the From the comparisons, it can be concluded that the proposed methods have high accurateness and efficiency. Moreover, instead of data-based methods, the present model is purely developed from random theory and statistics; therefore, it stands for the more general track geometry status, providing a foundation for random excitation simulation of train-track interactions.

DEM-FEM model
Instead of the expensive full-scale experiments, an alternative, a combination of discrete element method (DEM) and finite element method (FEM) to model the ballasted track system, is therefore developed to obtain sufficiently enough track properties to achieve statistical evaluation of uneven track parameters.
The DEM, developed by Cundall and Strack [31], provides an effective and reliable approach to simulate the granular characteristic of ballast assemblies. For the railway ballasted track, the random property of ballast layer is much larger than that of the subgrade. It is a fact that the  coarse granular aggregates are difficult to compact uniformly. Besides, research from the EUROBOLT project recommended that the variation of the subgrade stiffness should be limited to less than 10% of the mean value [32]. Therefore, regardless of the random nature of the subgrade system, a railway ballast track DEM model consisting of rail, sleepers and ballasted bed is therefore developed to reveal the non-uniform property of the supporting stiffness of the ballast layer. Figure 7 shows the DEM model of railway ballasted track. The rail is simulated as a beam using the parallel bond, which can transmit both force and moment between particles, and the sleeper and fastener are modelled by bonding small particles together as a clump (super rigid body) [33]. Especially, the ballast layer is modelled by the compacted well-graded ballast particles with irregular shape, which is simulated by integrating the DEM with the image processing techniques. For more details of the DEM model of railway ballast track, one can refer to Ref. [34].
Applying DEM to establish a large-scale track model is yet impractical due to the large computational cost, thus only 57 sleepers are simulated in the DEM model. The length of the model is 35 m (including 57 sleepers) with the sleeper spacing interval of 0.6 m, and the thickness of ballast bed is 0.35 m. To obtain sufficiently enough track data of sleeper support stiffness, 19 identical DEM models of the same scale have been established simultaneously. The differences between these track models are the random characteristics of ballast particles in compactness, irregular shape, size and location, etc. Therefore, a total length of 1995-m ballast track can be obtained by assembling these segment tracks, which contacts 1083 sleepers.

Comparison with experimental results
The sleeper support stiffness is a comprehensive and basic index for evaluating track performance, and it can be defined as the capacity of a track to bear a load. Herein, the gradual increase load P (from 0 to 40 kN) is applied to a certain sleeper unfastened with the rail (removing the fastener particle), and the corresponding sleeper displacement y is recorded. Thus, the sleeper support stiffness under this sleeper can be expressed as 2P/y approximately (only half of the track is considered in the model). Examples of loaddeflection characteristics for four sleepers are shown in Fig. 8, the ''Sleeper 16'' means the No. 16 sleepers counted from left to right facing X-axis. Note that the deflection of sleepers varies significantly from 0.56 to 0.67 mm for the same load level and the sleeper stiffness values are in the range from 119 to 142 kN mm -1 , which is consistent with previous filed test results conducted by Oscarsson and Dahlberg [35,36] and Ma [37]. Using this DEM-FEM model, the statistical distribution of the support stiffness for the track can be obtained, as illustrated in Fig. 9, where the mean and the standard deviation of the support stiffness are 112.95 and 16.46 MN m -1 , respectively.

Random field of system parameters
Based on random field theory, the system parameter space is denoted by S. If H q,n is a set of dimension N, where N is the number of physical and mechanical quantities, then each quantity is a stochastic process of f(s) with a vector of dimension d. Random field of system parameters means a collection of random variables: According to the statistics of DEM-FEM model solutions, the PDF of sleeper support stiffness and ballast mass can be obtained.  The random sets of system parameters have already been given in the random vector H n 2 H.
For each random sample, the parameters related to the physical properties of systems will be considered in the dynamic matrices, where the parameters related to the wheel-rail contacts/creepages, e.g. the wheel-rail friction coefficient, all can be characterized in the derivation of the dynamic loads.

Framework for system stochastic analysis
Abovementioned work has fundamentally addressed the issues in the modelling of train-track interaction and the simulation of uncertain geometric, physical and mechanical excitations quantitatively and with statistical significance.
The modelling framework for train-track coupled stochastic analysis can be formed as shown in Fig. 10.

Numerical examples
In the numerical examples, a train consisting of three identical vehicles is set to run with a constant speed of 160 km/h. The main parameters of the ballast tracks and the vehicle are listed in ''Appendix 2''.

Uncertainty quantification of track geometric excitation
The most significantly important parts in constructing a random analysis model for train-track interactions have been illustrated. However, a validation should be further investigated: the practicability of the random simulation method for track irregularities in statistically evaluating the train-track dynamic behaviours. In this example, random analysis is mainly conducted from statistics and frequency prospects, in which the random simulation method of track irregularities and the train-track interaction model are integrated effectively.

PDF comparison
For comparisons, two forms of irregularity excitation are considered: • C 1 : Excitation of the measurement: a long length of track irregularities that are experimentally measured, as shown in Fig. 3, is used as the system excitation. • C 2 : Excitation of the simulation: the track irregularities simulated by the random simulation method, in which the irregularity correlation and amplitude and wavelength characteristics are extracted from the measured irregularities of C 1 , is used as the system excitation. Figure 11 plots the comparisons of the vertical acceleration of the bogie frame, the lateral acceleration of the wheelset and the wheel-rail forces between C 1 and C 2 . It can be observed from Fig. 11 that the PDFs of dynamic indices excited by the measured track irregularities coincide well with those excited by the simulated track irregularities. However, the sampling number of the track geometry sets used in C 2 is significantly smaller than that in C 1 (50 vs 230). Thus, the statistical distribution of the dynamic responses of train-track systems, where the amplitude and probability can be fully depicted, is efficiently revealed by the proposed models.   Fig. 13 System responses with respect to spatial variation of track parametric excitations (a rail displacement; b wheel-rail force)

PSD comparison
The frequency characteristics of system responses can also be given by the time-domain solutions through spectral or frequency analysis. Figure 12 shows the statistical average PSD of dynamic indices, in which the spectral comparison between C 1 and C 2 is presented within the effective frequency range ( V=l lim , where V is the vehicle speed, and l lim is the lower-limit wavelength of track irregularities). As observed from Fig. 12, both the characteristic wavelength and spectral value coincide well for different dynamic indices. Besides, it can be noticed that the spectral distributions of the rail displacement show shape similarity to those of wheel-rail forces. Therefore, the proposed track irregularity simulation method can also be used to characterize system responses in frequency domain.

Influence of spatially uneven track parametric excitation
In the example, the track geometric excitation is temporarily ignored for clearly displaying the effects of track parametric excitation. Besides, it is remarked that the track transverse vibration is generally rather smaller than that of the vertical vibration caused by track parametric excitation. The vertical support stiffness of rail pads and ballast layer and the mass track bed are regarded as the independent uncertain parameters. The distribution characteristics of track random parameters have been illustrated in Table 1. For parameters following the normal distribution, the coefficient of variation (COV) is set to be 0.1, 0.2 and 0.3, respectively. Figure 13 shows the time-varying rail vertical displacement and the wheel-rail vertical force with respect to the spatially uneven track parametric excitations on the first wheelset of the first vehicle. It can be observed from Fig. 13 that response amplitudes in both the rail vertical displacement and the wheel-rail vertical force are enlarged by the increase in COV of the random track parameters. For example, the maximum value of rail vertical displacement is 0.605, 0.652 and 0.734 mm with respect to COV of 0.1, 0.2 and 0.3, respectively. Moreover, the maximum wheel-rail vertical force has also been increased from 51.3 to 54.5 kN; the difference between the wheelrail minimum and maximum force has also been increased from 1.93 to 6.86 kN, i.e. a 2.55 times increase.
Moreover, Fig. 14 illustrates the time-varying interaction forces between system components regarding the moving of the first wheelset. One can see that the interaction forces between the rail and sleeper and between the track bed and subgrade surface have also been significantly increased due to the longitudinal uneven distribution of the track properties in mechanics and physics.
Obviously, it can be observed from Figs. 13 and 14 that the spatial uneven track parametric excitation possesses great influence on wheel-rail interaction and track vibrations. Figure 15 further shows the performance of car body vertical acceleration against different COV of track parameters variation, from which it can be seen that the changes of car body acceleration are noticeable, but with slight response amplitudes. See for instance, the maximum car body acceleration at COV of 0.3 is about 3.8 times that at COV of 0.1. In the above analysis, it can be cognized that the system behaviours are greatly influenced by the track geometric and parametric excitations. However, these two excitation sources show different degrees of influence to different system components. Generally, track geometric excitation, such as the rail profile irregularities, exerts direct effects on whole systems, but the track parametric excitation holds more dynamic impacts on track components than those on vehicle systems. In this example, an emphasis is therefore put on clarifying the contribution of different excitation sources on system performance, but mainly on vertical vibrations. The state-of-track irregularities is quantitatively represented by the probability level, denoted by c, of their equivalent spectral density functions, as elaborated in [3], and the variation of system parameters such as rail pad stiffness, sleeper support stiffness and track bed mass is quantitatively represented by the COV denoted by f.
Set c to be varied from 0.01 to 0.99, and f to be varied from 0.01 to 0.30. Figure 16 presents the results of car body acceleration against different combining excitations of track irregularity and system longitudinal uneven parameters. It can be observed that the car body acceleration is gradually increased by the enlargement of track irregularity probability level, since the larger of the probability level, the worse of the track irregularity state. However, the spatial variation of system parameters shows little influence on car body acceleration (Fig. 16b), because the car body acceleration changes little by the COV of system parameters. Figure 17 plots the distribution of wheel-rail vertical force with respect to the variation of track geometric and parametric excitation. With the same manner, the wheelrail vertical force is also increased by the enlargement of track irregularity excitation in general. When the probability level of track irregularity is bigger than 0.65, the wheel-rail force fluctuates violently, and scenarios of wheel-rail separation occur. Unlike the characteristics of car body acceleration responding to the COV of spatial parameters, the wheel-rail force is also increased by the spatial variation of track parameters: the larger of the COV of parameters, the severer of the wheel-rail interaction. When the track irregularity probability level is below 0.65, the wheel-rail interaction force is increased by about 9.57%. If the track irregularity probability level is above 0.65, the wheel-rail force might be increased by a larger extent.
Besides, responses of track vibration indices such as rail vertical displacement can be also derived by this work, as shown in Fig. 18. Obviously, the rail vibrations are more violently influenced by the spatial variation of track parameters, compared to the car body acceleration and the wheel-rail interaction force. For example, the rail displacement can be increased by about 17.16% when the COV of parameters grows from 0.01 to 0.30.

Conclusions
In this work, a systematic modelling framework is developed to achieve the stochastic analysis of train-track interaction with consideration of the spatial variation of system properties. With novelty, new simulation models are developed to quantitatively characterize the randomness of system properties in track irregularities and statistical properties of railway ballast in physics and mechanics.
Apart from model validations, some conclusions can be obtained from the numerical examples: (1) The spatial variation of system physical-mechanical parameters will affect the system vibration substantially, the larger of the coefficients of variation of system random parameters, the larger of the dispersion of system responses.
(2) The extent of influence caused by the system parametric variation is absolutely different for various system components. Generally, the track parametric variation exerts slight influence on vehicle responses compared to its influence on track vibrations. (3) System geometric variation, such as the track irregularities, plays an important role in figuring out the random behaviour of the train and track. In the presented cases, it can be confirmed that when the track irregularity probability level reaches 0.65 or above, its dynamic influence is remarkably larger than those brought out by track parametric variation.
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://creativecommons. Set z mÂn ¼ ðz 1 ; z 2 ; . . .; z n Þ as the random variable sequences with order of m Â n following a normal distribution, i.e. there are n variables and m samples, with the mean value vector of u ¼ ðu 1 ; u 2 ; . . .; u n Þ and covariance matrix C z , that is, . . .
The covariance matrix C z can be decomposed into the form of It can be derived that where Y mÂn ¼ Y 1 ; Y 2 ; . . .; Y n ð Þis the n dimensional independent standard normal distribution; Ç ½Á is an operator of deriving covariance matrix.
To variables of non-normal distribution, a nonlinear transformation between normal distribution and arbitrary distribution can be obtained by where f Z i ðz i Þ is the PDF of normal distribution; f V i ðv i Þ is the PDF of an arbitrary distribution; z i and v i denote, respectively, variable value following normal distribution represented by subscript 'Z i ' and arbitrary distribution represented by subscript 'V i '.
The correlation coefficient matrix in Equation (15) will be modified by the following steps: (1) The correlation coefficient between two variables, V i and V j , can be equivalently defined by where E Á ½ is the mathematical expectation operator; u denotes the mean value of a variable; and r denotes the standard deviation of a variable.
(2) Through the nonlinear transformations, Equation (18) can be then deduced as Following Equation (19), the mathematical relation between C V i ;V j and C Z i ;Z j is established, and a mapping relation between C z of normal distribution and C v of arbitrary distribution is obtained accordingly.
By substituting C z as the modified C v into Equation (16), and by Equation (17) the correlated pseudorandom vector can be obtained.