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

Abstract

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–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, thin-plate 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 multi-rigid-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 three-dimensional (3D) space. Regarding the differences in wheel–rail contact/creepage descriptions, i.e. rigid contact [7–9, 14–17], elastic contact [3, 18–20] and elastic-plastic contact [21–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 physics–mechanics uncertainty have been studied as a hot topic but generally in an independent way in the train–track 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

1

$$\left[\begin{array}{cc}{\mathit{M}}_{\text{VV}}& {\mathit{M}}_{\text{VT}}\\ {\mathit{M}}_{\text{TV}}& {\mathit{M}}_{\text{TT}}\end{array}\right]\left[\begin{array}{c}{\ddot{\mathit{X}}}_{\text{V}}\\ {\ddot{\mathit{X}}}_{\text{T}}\end{array}\right]+\left[\begin{array}{cc}{\mathit{C}}_{\text{VV}}& {\mathit{C}}_{\text{VT}}\\ {\mathit{C}}_{\text{TV}}& {\mathit{C}}_{\text{TT}}\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{X}}}_{\text{V}}\\ {\dot{\mathit{X}}}_{\text{T}}\end{array}\right]+\left[\begin{array}{cc}{\mathit{K}}_{\text{VV}}& {\mathit{K}}_{\text{VT}}\\ {\mathit{K}}_{\text{TV}}& {\mathit{K}}_{\text{TT}}\end{array}\right]\left[\begin{array}{c}{\mathit{X}}_{\text{V}}\\ {\mathit{X}}_{\text{T}}\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_{\text{V}}\\ {\mathit{F}}_{\text{T}}\end{array}\right],$$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; $\ddot{\mathit{X}}$, $\dot{\mathit{X}}$ and X denote the acceleration, velocity and displacement vectors, respectively; F denotes the loading vector; ${\mathit{M}}_{\text{VT}}$ and ${\mathit{M}}_{\text{TV}}$ are non-zero matrices for wheel–rail rigid contacts.

Fig. 1 [Images not available. See PDF.]

Train–track interaction model (a side view; b end view) (the symbols $x$, $y$, $z$, $\psi$, $\beta$ and $\theta$ 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)

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 ${\mathit{M}}_{\text{VV}}$, ${\mathit{C}}_{\text{VV}}$, ${\mathit{K}}_{\text{VV}}$ and ${\mathit{F}}_{\text{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-, Y- and 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. ${\mathit{M}}_{\text{TT}}$, ${\mathit{C}}_{\text{TT}}$ and ${\mathit{K}}_{\text{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].

Fig. 2 [Images not available. See PDF.]

Wheel–rail interaction model (the symbols ${\delta }_{\text{l}}$ and ${\delta }_{\text{r}}$ denote the wheel–rail left and right contact angle respectively; ${\theta }_{\text{ct}}$ denotes the rail cant, and $X_{\text{r}}$, $Y_{\text{r}}$, $Z_{\text{r}}$, ${\theta }_{\text{r}x}$, ${\theta }_{\text{r}y}$ and ${\theta }_{\text{r}z}$ denote the longitudinal displacement, the lateral displacement, the vertical displacement, the angle around X-axis, the angle around Y-axis and the angle around Z-axis of the rail centroid, respectively)

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

2

$$F_{\text{wr}}(\mathrm{\Delta }z)={\left(\frac{1}{G}\mathrm{\Delta }z\right)}^{3/2},$$where $G$ is the wheel–rail contact constant (m/N^2/3); $\mathrm{\Delta }z$ is the wheel–rail elastic compression (m).

From Eq. (2), the wheel–rail equivalent contact stiffness k_wr,z can be obtained:

3

$$k_{\text{wr},z}=\frac{F_{\text{wr}}(\mathrm{\Delta }z)}{\mathrm{\Delta }z},$$

Accordingly, the wheel–rail coupling matrix ${\mathit{K}}_{\text{wr}}$ can be obtained as

4

$${\mathit{K}}_{\text{wr}}=\sum _i\sum _{j=1}^4\sum _{g=1}^2{k}_{\text{wr,}z,j,g}{\mathit{N}}_{j,g}^{\text{T}}{\mathit{N}}_{j,g}$$with ${\mathit{N}}_{j,g}=\left\{\left[\begin{array}{cccc}1& {(-1)}^{g+1}{d}_{\text{lr}}& -1& -1\end{array}\right]|{\mathit{\Phi }}_{j,g}\right\}$, ${\mathit{\Phi }}_{j,g}=\left[\begin{array}{cccc}{z}_{\text{w},j}& {\phi }_{\text{w},j}& v_{\text{r},j,g}& z_{\text{r},j,g}\end{array}\right]$, $j=1,2,3,4$, $g=1,2$, where the symbol ${\mathit{\Phi }}_{j,g}$ denotes the displacement vector corresponding to the shape function; $i$ denotes the $i$th vehicle; $j$ denotes the $j$th wheelset and $g=1,2$ denote, respectively, the left side and right side of a wheelset; $d_{\text{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 $\phi$ 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.

Spatial uncertainty and quantification of system excitation on geometries, physics and mechanics

In course of a train–track interaction event, system excitations ($\mathit{\Theta }$), including track geometries and system parameters of physics and mechanics, show randomness in the spatial axle, which can be assembled by

5

$${\mathit{\Theta }}_q=({\mathit{\Theta }}_{q,\mathit{\varsigma }},{\mathit{\Theta }}_{q,\mathit{\xi }}),$$where $q$ denotes the sampling number: $q=1,2,\dots ,\overline{g}$, and $\overline{g}$ is the total number; the subscript $\mathit{\varsigma }$ denotes the track irregularity vector, and $\mathit{\xi }$ 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 multi-dimensional 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

6

$${\mathit{\Theta }}_{q,\mathit{\varsigma }}=({\lambda }_{\text{l}}(s),{\nu }_{\text{l}}(s),{\lambda }_{\text{r}}(s),{\nu }_{\text{r}}(s)),s\in (0,\overline{S}]$$with ${\lambda }_{\text{l,r}}={\gamma }_1\pm \frac{1}{2}{\gamma }_3$, ${\nu }_{\text{l},\text{r}}={\gamma }_2\pm \frac{1}{2}{\gamma }_4$,where ${\gamma }_1$, ${\gamma }_2$, ${\gamma }_3$ and ${\gamma }_4$ denote the vertical profile, alignment, cross-level and gauge irregularities, respectively; $\lambda$ and $\nu$ 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, $\overline{S}$ is the total length.

The PSD function of track irregularity vector can be expressed by

7

$${\mathit{P}}_{q,\mathit{\varsigma }}={\rm Y}({\mathit{\Theta }}_{q,\mathit{\varsigma }}),$$where ${\rm Y}(·)$ denotes the PSD estimator.

Spectral correlation of track irregularities

Divide the total length $\overline{S}$ into $N$ segments, i.e. $\overline{S}=\left((n-1)\overline{S}/N,n\overline{S}/N\right]$, $n=1,2,\dots ,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

8

$$\sum _{k=0}^{\upsilon -1}{\left|\mathit{\psi }(k\mathrm{\Delta }s)\right|}^2=\frac{\mathrm{\Delta }\omega }{\upsilon }\sum _{q=1}^{\upsilon }{\left|{\mathit{S}}_{\mathit{\psi }}(q\mathrm{\Delta }\omega )\right|}^2\approx \mathrm{\Delta }\omega \sum _{q=1}^{\upsilon }\mathit{P}(q\mathrm{\Delta }\omega ),$$where ${\mathit{S}}_{\psi }(·)$ denotes the discrete Fourier transform (DFT) of track irregularities; $\mathrm{\Delta }s$ and $\mathrm{\Delta }\omega$ denote the time and frequency interval, respectively; $\mathit{\psi }$ denotes the discrete time-domain track irregularity signals; $\upsilon$ denotes the total discrete number; $\mathit{P}$ denotes the PSD vector.

Equation (8) has exposed a close relation between time signals and its PSD. An integration indicator ${\mathit{\Phi }}_{\mathit{\psi }}(·)$ used to determine the PSD representing certain track geometric status can be therefore defined by

9

$${\mathit{\Phi }}_{\mathit{\psi }}\left(b\right)=\mathrm{\Delta }\omega \sum _{q=0}^{M-1}{\mathit{P}}_{\mathit{\psi }}(q\mathrm{\Delta }\omega |b),$$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

10

$$C_{i,j}=\frac{\text{Cov}({\mathit{\Phi }}_i,{\mathit{\Phi }}_j)}{\sqrt{\text{Var}[{\mathit{\Phi }}_i]\text{Var}[{\mathit{\Phi }}_j]}},(i,j)\in [{\gamma }_{\text{l}},{\nu }_{\text{l}},{\gamma }_{\text{r}},{\nu }_{\text{r}}],i\ne j,$$where $C_{i,j}$ denotes correlation coefficient; ${\mathit{\Phi }}_i$ denotes the integration indicator, $\text{Cov}(·)$ denotes the covariance function, and $\text{Var}(·)$ denotes the variance function

Based on Eq. (10), the correlation coefficient matrix (CCM) for the track irregularity integration indicator Φ_i can be obtained as

11

$${\mathit{R}}_{\Psi }=\left[\begin{array}{cccc}1& C_{{\lambda }_{\text{l}},{\lambda }_{\text{r}}}& C_{{\lambda }_{\text{l}},{\nu }_{\text{l}}}& C_{{\lambda }_{\text{l}},{\nu }_{\text{r}}}\\ C_{{\lambda }_{\text{r}},{\lambda }_{\text{l}}}& 1& C_{{\lambda }_{\text{r}},{\nu }_{\text{l}}}& C_{{\lambda }_{\text{r}},{\nu }_{\text{r}}}\\ C_{{\nu }_{\text{l}},{\lambda }_{\text{l}}}& C_{{\nu }_{\text{l}},{\lambda }_{\text{r}}}& 1& C_{{\nu }_{\text{l}},{\nu }_{\text{r}}}\\ C_{{\nu }_{\text{r}},{\lambda }_{\text{l}}}& C_{{\nu }_{\text{r}},{\lambda }_{\text{r}}}& C_{{\nu }_{\text{r}},{\nu }_{\text{l}}}& 1\end{array}\right],C_{i,j}=C_{j,i},$$where the subscript Ψ denotes the distribution type of the track irregularity integration indicator; ‘$C_{i,j}$’ with subscript $i,j=({\lambda }_{\text{l}},{\lambda }_{\text{r}},{\nu }_{\text{l}},{\nu }_{\text{r}})$ denotes the correlation coefficient between the $i$th and the $j$th 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 ${\mathit{\Phi }}_{\mathit{\psi }}(·)$. Generally, the statistical PDF of ${\mathit{\Phi }}_{\mathit{\psi }}(·)$ differs from a normal distribution. The method presented in Sect. 3.1 can be therefore applied to choose the PSD samples ${\mathit{\Phi }}_{\mathit{\psi }}(·)$.

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; ${\epsilon }_0$ is the original value assigned to $E$; $\epsilon$ is a tiny value, $\epsilon <{\epsilon }_0$; Z_n×m is a real number matrix with order n × m for the variable z following normal distribution; $\wp$ denotes the random simulation operator for standard normal distribution; ${\mathit{R}}_{\Psi \to \text{N}}$ denotes the CCMs correlation coefficient matrices transformed from arbitrary distribution ($\Psi$) to normal distribution (N); $\text{Chol}(·)$ denotes the Cholesky decomposition and $\text{Cor}(·)$ denotes the operator of correlation coefficient.

Because the track irregularities ${\lambda }_{\text{l,r}}(s)$ and ${\nu }_{\text{l,r}}(s)$ can be realistically measured, the PDF of the PSD indicator ${\mathit{\Phi }}_{\mathit{\psi }}(b)$, $b=1,2,\dots ,N$, can therefore be obtained statistically and discretely, represented as $f_{\mathit{\Phi }}(b|\mathit{\psi })$.

With acquisition of $f_{\mathit{\Phi }}(b|\mathit{\psi })$, the random simulation of the PSD indicator can be calculated by

12

$${\mathit{\Phi }}_{\mathit{\psi }}(b)=F_{\mathit{\Phi }}^{-1}(F_{{\overline{Z}}_{n\times m}}({\overline{z}}_i)),$$ where F denotes the cumulative distribution function of the PSD indicator. Obviously ${\mathit{\Phi }}_{\mathit{\psi }}(b)$ corresponds to the PSD function ${\mathit{P}}_{\mathit{\psi }}(\omega )$, 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 ${\mathit{P}}_{\mathit{\psi }}$; $M_{\mathit{il}}$ is the length of the time-domain track irregularities; $N_r$ is the sampling number; ${\overline{f}}_y$ is the frequency range; $f_{\text{l}}$ and $f_{\text{u}}$ denote the lower- and upper-limit frequency; $\text{conj}(·)$, $\text{interp}1(·)$ and $\text{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 $\phi$ 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.

Model validations

Figure 3 plots a section of track irregularities with a length of 230 km, which are measured by a track inspection car from a ballasted railway in China.

Fig. 3 [Images not available. See PDF.]

The measured track irregularities

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.

The detailed values of CCM are presented as follows:

13

$${\mathit{R}}_{\Psi }=\left[\begin{array}{cccc}1& 0.899& 0.213& 0.153\\ 0.899& 1& 0.173& 0.114\\ 0.213& 0.173& 1& 0.793\\ 0.153& 0.114& 0.793& 1\end{array}\right]\Rightarrow {\mathit{R}}_{\text{N}}=\left[\begin{array}{cccc}1& 0.918& 0.224& 0.167\\ 0.918& 1& 0.182& 0.125\\ 0.224& 0.182& 1& 0.835\\ 0.167& 0.125& 0.835& 1\end{array}\right].$$

Based on the methods presented above, the random samples of PSD indicators can be selected by Algorithm 1, as shown in Fig. 4, where ${\varphi }_{v_{\text{l}}}$, ${\varphi }_{v_{\text{r}}}$ and ${\varphi }_{{\gamma }_{\text{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_{{\nu }_{\text{l}},{\nu }_{\text{r}}}$ and $C_{{\nu }_{\text{l}},{\gamma }_{\text{l}}}$ are, respectively, with values of 0.914 and 0.222, which are close to the reference values 0.918 and 0.224, respectively.

Fig. 4 [Images not available. See PDF.]

Random selection of PSD indicator satisfying the correlation coefficient matrix (a correlation between the vertical irregularities at the left and right sides; b correlation between the vertical irregularities at left side and the lateral irregularities at the right side)

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 developed model, Figs. 5 and 6 show the comparisons of probility distribution function (PDF) and joint PDF of track irregularities that are, respectively, derived from simulation and measurement. From these comparisons, one can see obviously that the statistical distributions of simulation and measurement coincide well with each other, though with slight deviations inevitably. However, there are only 50 samples used in the random model, which is significantly lower than that of measurement (230 samples).

Fig. 5 [Images not available. See PDF.]

PDF comparison between simulation and measurement (a vertical irregularity ${\nu }_{\text{r}}$; b lateral irregularity ${\gamma }_{\text{r}}$)

Fig. 6 [Images not available. See PDF.]

Joint PDF comparison between vertical irregularities at left and right sides (a simulation; b measurement)

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.

Random simulation for system parameters

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].

Fig. 7 [Images not available. See PDF.]

The DEM for railway ballasted track

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 load-deflection 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⁻¹, 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⁻¹, respectively.

Fig. 8 [Images not available. See PDF.]

Vertical displacement–force curves of typical sleepers

Fig. 9 [Images not available. See PDF.]

Support stiffness distribution

Random field of system parameters

Based on random field theory, the system parameter space is denoted by S. If Θ_q,ξ 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:

14

$$\left\{{\mathit{\Theta }}_{\mathit{\xi }}\equiv f(s):s\in S\right\}.$$

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 Θ_ξ$\in$Θ.

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.

Fig. 10 [Images not available. See PDF.]

Modelling framework for stochastic dynamic analysis of train–track systems

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₁: 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₂: 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₁, 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₁ and C₂. 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₂ is significantly smaller than that in C₁ (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. 11 [Images not available. See PDF.]

PDF comparison between excitation types C₁ and C₂ (a vertical acceleration of the bogie frame; b lateral acceleration of the wheelset; c wheel–rail vertical force; d wheel–rail lateral force; e rail lateral displacement; f rail vertical displacement)

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₁ and C₂ is presented within the effective frequency range ($\le 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.

Fig. 12 [Images not available. See PDF.]

Average PSD comparison between C₁ and C₂ (a wheel–rail vertical force; b wheel–rail lateral force; c vertical displacement of the rail; d lateral displacement of the rail)

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.

Table 1. Parametric characteristics of the tracks

Track parameters	Mean parametric value
Support stiffness (MN/m)
Rail pad	58
Track bed	112.95
Density (kg/m3)
Track bed	1800

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 wheel–rail minimum and maximum force has also been increased from 1.93 to 6.86 kN, i.e. a 2.55 times increase.

Fig. 13 [Images not available. See PDF.]

System responses with respect to spatial variation of track parametric excitations (a rail displacement; b wheel–rail force)

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.

Fig. 14 [Images not available. See PDF.]

Interaction forces between system components (a rail–sleeper interaction force; b track bed–subgrade interaction force)

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.

Fig. 15 [Images not available. See PDF.]

Car body vertical acceleration

Contribution of track geometric and parametric excitations on system performance

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 $\gamma$, 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 $\zeta$.

Set $\gamma$ to be varied from 0.01 to 0.99, and $\zeta$ 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.

Fig. 16 [Images not available. See PDF.]

Car body vertical acceleration against the changes of track irregularity probability level $\gamma$ and parameter COV $\zeta$ (a 3D view; b side view; c end view)

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 wheel–rail 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.

Fig. 17 [Images not available. See PDF.]

Wheel–rail vertical acceleration against the changes of track irregularity probability level $\gamma$ and parameter COV $\zeta$ (a 3D view; b side view; c end view)

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.

Fig. 18 [Images not available. See PDF.]

Rail vertical displacement against the changes of track irregularity probability level $\gamma$ and parameter COV $\zeta$ (a 3D view; b side view; c end view)

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.

Acknowledgements

This work is supported by National Natural Science Foundation of China (NSFC) under Grant Nos. 51735012 and 11790283. High appreciation is given to Can Shi, Train and Track Research Institute at Southwest Jiaotong University, for helping to establish the DEM-FEM model to obtain the spatial variation distribution of track support stiffness in the numerical simulation.

References

1. Garg, VK; Dukkipati, RV. Dynamics of railway vehicle systems; 1984; Ontario, Academic Press:

2. EN14363:2005 (2005) Railway applications-testing for the acceptance of running characteristics of railway vehicles-testing of running behavior and stationary tests. CEN, Brussels, Beligium

3. Xu, L; Zhai, W; Gao, J. A probabilistic model for track random irregularities in vehicle/track coupled dynamics. Appl Math Model; 2017; 51, pp. 145-158.MathSciNet ID: 3694528zbMath ID: 07166255[DOI: https://dx.doi.org/10.1016/j.apm.2017.06.027]

4. Perrin, G; Soize, C; Duhamel, D; Fünfschilling, C. Track irregularities stochastic modelling. Probab Eng Mech; 2013; 34, pp. 123-130. [DOI: https://dx.doi.org/10.1016/j.probengmech.2013.08.006]

5. Náprstek, J; Frýba, L. Stochastic modelling of track and its substructure. Veh Sys Dy; 1995; 24, sup 1 pp. 297-310. [DOI: https://dx.doi.org/10.1080/00423119508969632]

6. Fröhling, RD; Scheffel, H; Ebersöhn, W. The vertical dynamic response of a rail vehicle caused by track stiffness variations along the track. Veh Syst Dyn; 1996; 25, Suppl pp. 175-187. [DOI: https://dx.doi.org/10.1080/00423119608969194]

7. Oscarsson, J; Nielsen, JCO. Dynamic train/track interaction with state-dependent track properties, technical report F215; 1999; Göteborg, Sweden, Department of Solid Mechanics, Chalmers University of Technology:

8. Oscarsson, J. Dynamic train–track-ballast interaction with unevenly distributed track properties. Veh Syst Dyn; 2002; 37, 1 pp. 385-396. [DOI: https://dx.doi.org/10.1080/00423114.2002.11666248]

9. Nielsen, JCO; Oscarsson, J. Simulation of dynamic train–track interaction with state-dependent track properties. J Sound Vib; 2004; 275, 3–5 pp. 515-532. [DOI: https://dx.doi.org/10.1016/j.jsv.2003.06.033]

10. Andersen, L; Øren, S; Nielsen, RK. Vibrations of a track caused by variation of the foundation stiffness. Probab Eng Mech; 2003; 18, pp. 171-184. [DOI: https://dx.doi.org/10.1016/S0266-8920(03)00012-2]

11. Wu, TX; Thompson, DJ. Influence of random sleeper spacing and ballast stiffness on the vibration behaviour of railway track. Acta Acust United Acust; 2000; 86, 2 pp. 313-321.

12. Dahlberg, T. Railway track stiffness variations: consequences and countermeasures. Int J Civil Eng; 2010; 8, 1 pp. 1-12.MathSciNet ID: 3404460

13. Sañudo, R; dell’Olio, L; Casado, JA; Carrascal, IA; Diego, S. Track transitions in railways: a review. Constr Build Mater; 2016; 112, pp. 140-157. [DOI: https://dx.doi.org/10.1016/j.conbuildmat.2016.02.084]

14. Nielsen, JCO; Abrahamsson, TJS. Coupling of physical and modal components for analysis of moving non-linear dynamic systems on general beam structures. Int J Numer Meth Eng; 1992; 33, 9 pp. 1843-1859. [DOI: https://dx.doi.org/10.1002/nme.1620330906]

15. Yang, YB; Chang, CH; Yau, JD. An element for analysing vehicle-bridge systems considering vehicle’s pitching effect. Int J Numer Method Eng; 1999; 46, pp. 1031-1047.zbMath ID: 0964.70004[DOI: https://dx.doi.org/10.1002/(SICI)1097-0207(19991110)46:7<1031::AID-NME738>3.0.CO;2-V]

16. Zhai, W; Han, Z; Chen, Z; Ling, L; Zhu, S. train–track-bridge dynamic interaction: a state-of-the-art review. Veh Syst Dyn; 2019; 57, 7 pp. 984-1027. [DOI: https://dx.doi.org/10.1080/00423114.2019.1605085]

17. Chen, XM; Deng, X; Xu, L. A three-dimensional dynamic model for railway vehicle-track interactions. J Comput Nonlinear Dyn; 2018; 13, 071006. [DOI: https://dx.doi.org/10.1115/1.4040254]

18. Zhai, W. Vehicle-track coupled dynamics theory and applications; 2020; Singapore, Springer: [DOI: https://dx.doi.org/10.1007/978-981-32-9283-3]

19. Zhai, W; Wang, K; Cai, C. Fundamentals of vehicle-track coupled dynamics. Veh Syst Dyn; 2009; 47, 11 pp. 1349-1376. [DOI: https://dx.doi.org/10.1080/00423110802621561]

20. Arvidsson, T; Andersson, A; Karoumi, R. Train running safety on non-ballasted bridges. Int J Rail Transp; 2019; 7, 1 pp. 1-22. [DOI: https://dx.doi.org/10.1080/23248378.2018.1503975]

21. Zhao, X; Li, Z. The solution of frictional wheel–rail rolling contact with a 3D transient finite element model: Validation and error analysis. Wear; 2011; 271, pp. 444-452. [DOI: https://dx.doi.org/10.1016/j.wear.2010.10.007]

22. Vermeulen, PJ; Johnson, KL. Contact of non-spherical bodies transmitting tangential forces. J Appl Mech; 1964; 31, pp. 338-340.zbMath ID: 0133.18105[DOI: https://dx.doi.org/10.1115/1.3629610]

23. Sladkowski, A; Sitarz, M. Analysis of wheel–rail interaction using FE Software. Wear; 2005; 258, pp. 1217-1223. [DOI: https://dx.doi.org/10.1016/j.wear.2004.03.032]

24. Xu, L; Zhai, W. A three-dimensional dynamic model for train–track interactions. Appl Math Model; 2019; 76, pp. 443-465.MathSciNet ID: 3979362[DOI: https://dx.doi.org/10.1016/j.apm.2019.04.037]

25. Xu, L; Yu, ZW; Shi, C. A matrix coupled model for vehicle-slab track-subgrade interactions at 3-D space. Soil Dyn Earthq Eng; 2020; 128, 105894. [DOI: https://dx.doi.org/10.1016/j.soildyn.2019.105894]

26. Kalker, JJ. Survey of wheel–rail rolling contact theory. Veh Syst Dyn; 1979; 5, pp. 317-358. [DOI: https://dx.doi.org/10.1080/00423117908968610]

27. Vermeulen, PJ; Johnson, KL. Contact of nonspherical elastic bodies transmitting tangential forces. J Appl Mech; 1964; 86, pp. 338-340.zbMath ID: 0133.18105[DOI: https://dx.doi.org/10.1115/1.3629610]

28. Wen, DZ; Zhuo, RH; Ding, DJ; Zheng, H; Cheng, J; Li, ZH. Generation of correlated pseudorandom variables in Monte Carlo simulation. Acta Phys Sin; 2012; 61, 22 220204.(in Chinese)

29. Golub, GH; Van Loan, CF. Matrix computations; 1996; 3 Baltimore, Johns Hopkins:zbMath ID: 0865.65009

30. Press, W; Flannery, B; Teukolsky, S et al. Numerical recipes in C: the art of scientific computing; 1992; 2 New York, Cambridge University Press:zbMath ID: 0845.65001

31. Cundall, PA; Strack, ODL. A discrete numerical model for granular assemblies. Geotechnique; 1979; 29, 1 pp. 47-65. [DOI: https://dx.doi.org/10.1680/geot.1979.29.1.47]

32. Hunt GA (2000) EUROBALT optimises ballasted track, Railway Gazette International: 813–816

33. Itasca Consulting Group, Inc. (2004) Particle flow code in 2 dimensions, version 3.1, user’s mannual

34. Zhang, X; Zhao, C; Zhai, W. Dynamic behavior analysis of high-speed railway ballast under moving vehicle loads using discrete element method. Int J Geomech; 2016; 17, 7 04016157. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0000871]

35. Oscarsson, J; Dahlberg, T. Dynamic train track ballast interaction computer models and full scale experiments. Veh Syst Dyn; 1998; 29, pp. 73-84. [DOI: https://dx.doi.org/10.1080/00423119808969553]

36. Oscarsson, J. Dynamic train–track interaction: variability attributable to scatter in the track properties. Veh Syst Dyn; 2002; 37, 1 pp. 59-79. [DOI: https://dx.doi.org/10.1076/vesd.37.1.59.3538]

37. Ma, C. Research on evaluation index optimization method of high-speed railway ballasted bed. Railw Stand Des; 2016; 60, 5 pp. 20-24.