Inconsistent effect of dynamic load waveform on macro- and micro-scale responses of ballast bed characterized in individual cycle: a numerical study

Cyclic load is widely adopted in laboratory to simulate the effect of train load on ballast bed. The effectiveness of such load equivalence is usually testified by having similar results of key concerns of ballast bed, such as deformation or stiffness, while the consistency of particle scale characteristics under two loading patterns is rarely examined, which is insufficient to well-understand and use the load simplification. In this study, a previous laboratory model test of ballast bed under cyclic load is rebuilt using 3D discrete element method (DEM), which is validated by dynamic responses monitored by high-resolution sensors. Then, train load having the same magnitude and amplitude as the cyclic load is applied in the numerical model to obtain the statistical characteristics of inter-particle contact force and particle movements in ballast bed. The results show that particle scale responses under two loading patterns could have quite deviation, even when macro-scale responses of ballast bed under two loading patterns are very close. This inconsistency indicates that the application scale of the DEM model should not exceed the validation scale. Moreover, it is important to examine multiscale responses to validate the effectiveness of load simplification.


Introduction
Discrete element method (DEM) has been applied in modelling ballast for decades.The application of DEM enables us to observe the ensemble of ballast particles [1,2].In DEM modelling, each particle serves like a sensor to obtain the tremendous details in ballast bed, including inter-particle contact force [3,4], particle movements [5,6], etc. Besides, various influence factors and complex conditions can be well-controlled in DEM modelling.A key concern of modelling ballast bed is to validate the effectiveness of DEM model.Generally, a numerical model of ballast bed can be validated by comparing with the results from laboratory or in situ check tests.In laboratory test, it is difficult to reproduce train load, especially high-speed train load.Hence, cyclic load is widely adopted to simulate the effect of train load in experimental study of ballast behavior [7,8,9] or ballast bed responses [10,11], especially for tests with representative elements, such as the triaxial test, direct shear test, and ballast box test.The effectiveness of such load equivalence is usually testified by having similar results of key concerns of ballast bed, such as deformation [12], stiffness [13] or strength [14,15,16], sleeper displacement [17], and acceleration [18], while the consistency of particle scale characteristics under two loading patterns is rarely examined.Since the same or similar dynamic behaviors of ballast bed from limited cases could result from very different microscopic processes, numerical modelling without examining particle scale responses tends to be insufficient in well-understanding the microscopic behaviors of ballast bed [5] and furtherly the equivalence of train load and cyclic load.
The field or laboratory measurement of ballast particles' behavior is the pre-condition to validate a DEM model of ballast bed.Initially, wired accelerator was employed to capture the acceleration characteristic in ballast bed [19,20].Recently, wireless sensing stone in particle size has been employed to characterize the movement of individual particles in ballast bed [21,22].The obtained ballast particle movements show great potential in nondestructive disease identification [23], fouling evaluation [22,24], and reinforcement design [25] of ballast bed.Ballast movements measured by sensing stone provide particle-scale indicators to validate the DEM model concerning railway ballast [5,26].By now, ballast particle movements have been obtained under train load in situ [23].In laboratory test, ballast particle movements are obtained mainly under cyclic load [6,27], rarely under train load [28], while DEM modelling of ballast bed is validated mainly by particle movements obtained under cyclic load [5,6].This owes mainly to that laboratory test is usually easily controlled and more productive in data collection.Besides, a representative segment is usually adopted in DEM modelling due to the limitation of computation capacity of computers, whose constraints are more similar to those of laboratory test [2,4,15,29].In order to well-understand and use the load simplification in exploring the micro-mechanism of ballast bed behavior, the consistency of particle scale characteristics under cyclic load and train load is worth being examined.
In this study, a full-scale ballast bed model with a half sleeper is built via PFC3D and validated by the dynamic responses monitored by high-resolution sensors in a previous model test [5,15,27].Then, two widely used waveforms, i.e., cyclic load and train load, are applied in the numerical model.According to the equivalence of loading effect, cyclic load and train load are assigned the same magnitude and amplitude for comparison [12].Then, we obtained not only macroscopic responses like dynamic stiffness and subballast stress, but also statistical characteristics of inter-particle contact force and particle movements in ballast bed.By comparing the values of each indicator under two loading patterns, the influence of loading pattern on multiscale responses is discussed.

DEM modelling of ballast bed
Considering the symmetry of ballast bed and the loading capacity, a half sleeper model is built in laboratory, as shown in Fig. 1a.The establishment details and test schemes of the full-scale laboratory test have been reported in previous work [27,29].The gradation of reproduced ballast conforms to the first class gradation requirements of the Chinese criterion 'Railway Ballast (TB/T 2140-2008)', as shown in Fig. 1b.Cyclic load is applied by an electrohydraulic servo actuator.The loading frequency varies from 1 to 18 Hz, with a loading magnitude of 30 kN.
The half sleeper model shown in Fig. 1 is built via PFC 3D, as shown in Fig. 2. Since single track railway line is considered, the bottom of ballast bed is about 2.5 m in width (along x direction), with a slope of 1:1.75.Ballast aggregate beneath the sleeper is approximately 0.35 m in height (along the z direction) and 0.6 m in thickness (along the y direction).The model contains more than 11,000 ballast particles, where particles are reproduced on the basis of the profiles scanned from real ballast.Here, three dominant ballast shapes are considered: quasi-cubic type, flat type and elongated type, whose aspect ratios are about 2:2:1, 1:1:1 and 2:1:1, respectively.The shape database includes shape samples from 20 real ballast particles.For each diameter range shown in Fig. 1b (for example 25-36 mm), particles are generated with shape randomly assigned from the shape database, whose size follows a normal distribution.The simulation accuracy of ballast particles is governed mainly by the ratio of the smallest to largest pebble and the distance (angular measure) between overlapping pebbles in the bubble pack algorithm [30].In this study, the ratio and distance are taken as 0.4 and 100°, respectively.
The ballast bed is filled in three layers, where enough ballast particles are generated and preloaded for each layer.The final bulk density of ballast bed is controlled no lower than 1750 kg/m 3 , which conforms to the density requirement in the Chinese criterion 'Railway Ballast (TB/T 2140-2008)'.The sleeper is composed of mono-size balls in rectangular arrangement, where the parallel contact model is adopted to allow deflection of the sleeper.A loading plate is fixed on the sleeper at the sub-rail location.Theoretically, the symmetry plane of a sleeper is deemed having no rotation under symmetric wheel load, and balls at the cut-section of the half sleeper are allowed only free vertical translation.The numerical parameters of ballast are calibrated by triaxial test [5].The whole numerical model is validated by macro-scale responses obtained in laboratory test [27], including sub-ballast stress and dynamic stiffness of the ballast bed.The model parameters are finally determined as shown in Table 1, where damping coefficients take the values suggested in Ref. [5].Four walls are set in DEM model to simulate the boundary condition shown in Fig. 1b: The normal and tangential contact stiffnesses for the two walls along the y direction are all 2.3 × 10 6 N/m, which are set as 1 × 10 7 N/m for the other walls.
Figure 3 shows the comparison of macroscopic responses obtained in laboratory test shown in Fig. 1a and  In the following, the validated DEM model is adopted to obtain multiscale responses of ballast bed under moving cyclic load and train load.The cyclic load applied in DEM modelling is identical to that in the laboratory test [27,29], given by where P d is the applied cyclic load, P amp is the loading amplitude, P 0 is to simulate the static seating force in laboratory, and f is loading frequency.Here we have P amp = 24 kN and P 0 = 6 kN.
The train load applied in DEM modelling is formulated by series concerning train speed and length [31], given by where P axle is the axle load, α is the spacing between adjacent sleepers, β is the stiffness ratio between rail and sub-rail base, x i is the horizontal distance between wheel i and the rail seat, L is the length of carriage, t is loading time shown in Fig. 4, v is train speed, and < λ > is the integer component of λ.Here, we have P 0 = 6 kN, P axle = 140 kN, α = 0.6, β = 1.18,L = 25 m, i = 1-4, x 1 -x 4 are −5.0,−2.5, 2.5 and 5.0, respectively (i.e., the target sleeper locates in the middle of two carriages), and v = 250 km/h.
The loading frequency of the cyclic load is usually estimated from the train speed by f = kv/L, where k is the loading cycles per carriage.Since a carriage usually has two bogies, we have k = 2 (with the assumption of even bogie spacing) to simulate the dynamic effect from train load [12].Then, the loading frequency f in Eq. ( 1) is calculated as 5.5 Hz.The applied loading curves are shown in Fig. 4. It is observed that the cyclic load and train load have the same magnitude.
Generally, applying measured dynamic loads or a numerical load simulated from the train-track dynamic interaction matches the real condition better [12,32].Since the train load calculated by Eqs. ( 2)-( 3) differs a lot from the cyclic load, as shown in Fig. 4, it is qualified to be used (1)  in investigating the influence of load pattern on multiscale responses of ballast bed.Before applying dynamic load, the ballast bed model (see Fig. 2) has reached the initial equilibrium state.After then, the two loading waves are applied on sleeper to start the loading procedure.In other words, the two numerical tests are exactly the same before loading.The following analysis characterizes the loading-induced variation of multiscale responses.Since the DEM model has more than 11,000 ballast particles, 22 loading cycles are applied within 4 s to obtain a reasonable computation efficiency.There are many macro-scale indicators to characterize the mechanical behaviors of ballast bed, such as the elastic vertical deformation, permanent vertical deformation, stress distribution, stiffness, and lateral confinement loss.Since only 22 loading cycles are applied in DEM modelling, indicators having insignificant variation with loading cycles are analyzed, including sub-ballast stress and ballast bed stiffness.

Quantification of particle-scale characteristics
Experimental studies show that the dynamic load from sleeper follows a diffusion angle of 35° in ballast bed [19,29].Since the thickness of ballast bed beneath the sleeper is 0.35 m, five typical sub-regions with side length of 35 cm in ballast bed along lateral direction are enough to cover the core region, as shown in Fig. 5.The 5 regions include three sub-rail regions (R 1 , R 2 and R 3 ), a sleeper-end region (R 4 ), and a shoulder region (R 5 ).For each region, multiscale responses are characterized.Meanwhile, in order to emphasize the influence of dynamic load on ballast behavior, the dynamic component of each response is characterized.In this section, quantifications of micro-scale indicators adopted in this work are introduced.The specific results are analyzed and interpreted in Sect. 4. Such arrangement of °Fig.
5 Five sub-regions in ballast bed along lateral direction content helps to compare the inconsistent variation of different indicators simultaneously.

Particle acceleration
The translational acceleration of the probe particle 10 cm beneath the sleeper end is taken as examples to show the determination of particle acceleration amplitude ( A i,j ).Particle acceleration and rotation movements become stable after the initial 6 loading cycles.Figure 6 shows the variation of particle acceleration with loading time during the stable phase, where the direction of rotation motion conforms to the right-hand rule.It is observed that particle acceleration has the same or similar waveform under the applied load after initial several loading cycles.It is observed from Fig. 6 that for both cyclic and train load: a x and a y are close, which are smaller than a z .This indicates that the unbalanced force is mainly in z direction, i.e., the loading direction.For a loading cycle, the amplitude of acceleration is characterized as the difference between the maximum and minimum values, as highlighted in Fig. 6.
Here, 'a loading cycle' for train load is taken as the effect from a bogie.Since a single particle's acceleration may not be representative, the regional mean amplitude of particles' translational accelerations is adopted to evaluate the motion intensity of each region.For a sub-range R i , the mean amplitude of particles' translational acceleration is calculated by Eq. ( 4): where i is the sub-range ID, j is the particle ID in R i ; N P i is the quantity of ballast particles in R i ; A i,j refers to the (4) amplitude of the translational acceleration on particle j during a loading cycle; A i describes the mean acceleration level in R i .

Contact network
Contact network is generally defined to describe the internal structure of granular material.There are many quantifications to describe contact network, such as the directional distribution of contact normal and contact force [33,34,35], statistical analysis of branch vector [36], coordinate number, and fabric tensor [37].In this study, the coordinate number, contact normal and contact force are adopted to characterize contact network, where contact normal (force) refers to the normal contact vector (force) at each contact.The mean coordinate number of each region is defined in Eq. ( 5), which describes the contact efficiency in R i ; the regional mean contact force is defined in Eq. ( 6), which evaluates the contact force level in R i .
where c i is the mean coordinate number of R i , N C i,j is the quantity of active contacts on particle j in R i , F i is the mean normal contact force of R i , N C i is the contact quantity on particles in R i , and F i,j is the normal contact force on contact j subjected to particles in R i .Figure 7a and b shows the variation of F i under cyclic loading and train load, respectively, which show similar waveforms as that of the applied loads.Then, the difference between the maximum and minimum ( 5) .Despite the mean contact force or coordinate number, directional distribution of contact normal is also an important aspect to understand the performance of contact network.For 3D particle assembly, the directional distribution of contact force is usually described by probability density shown in a spherical histogram.However, due to the small quantity of ballast particles in each sub-region, the 3D directional distribution of contact force shows great difference among directions even a large interval of steradian is adopted [4].Considering that the ballast bed is a periodic structure in the longitudinal direction, the contact forces are projected on the x-z plane to mainly analyze the directional distribution of contact force in vertical and lateral directions.The probability of contact normal oriented in the angle interval [θ 1 , θ 2 ] is expressed by P c (θ), as defined in Eq. ( 7) [33].Meanwhile, the directional probability distribution of contact normal is weighted by force intensity to better consider load propagation, see P f (θ) defined in Eq. ( 8).
where N( − Δ ∕ 2, + Δ ∕ 2) is the number of contact normal oriented in the interval − Δ ∕ 2, + Δ ∕ 2 , and F( − Δ ∕ 2, + Δ ∕ 2) is the arithmetic sum of forces ori- ented in the interval − Δ ∕ 2, + Δ ∕ 2 .In the following analysis, Δ is selected as 10º.( 7) With P c (θ) and P f (θ), the probability density function (PDF) describing the directional distribution of contact normal or contact force is calculated by derivation [33], which is generally given by where i could be 'c' or 'f' to describe the distribution of contact normal or contact force, respectively, θ i is the major direction of the probability density distribution, and δ i is the deviator of the probability density distribution.Larger δ i indicates a more anisotropic distribution of p i (θ). Figure 8a and b shows the PDF of contact force of R 3 at a moment during cyclic loading and train load, respectively.It is observed that the directional probability density distributions of contact force for two loading patterns have the same dominant direction (loading direction), which indicates that the contact force among ballast particles mainly locate along the loading direction.
The contact evolution is characterized by the major direction (θ i ) and deviator (δ i ).For all five regions, the deflection of θ i induced by dynamic loading is no higher than 10°, which is even smaller than the statistical angle interval ( Δ ) of the probability density distribution.Hence, θ i is not discussed in the following analysis.In other words, the major direction of contact network is deemed the same throughout the loading cycle.Figure 9 shows the variation of δ f under cyclic load and train load, respectively, which show similar waveforms as that of the applied loads.Then, the difference between the maximum and minimum values of δ f during a loading cycle is characterized as the amplitude of δ f , denoted as amp f (also highlighted in Fig. 9).The same process is used to obtain   With the above quantifications of particle-scale characteristics and conventional evaluation indicators (sub-ballast stress and dynamic stiffness of ballast bed), the influence of dynamic load waveform on macro-and micro-scale characteristics of ballast bed is quantitatively discussed in the following section.For comparison purpose, all indicators are collected at the same loading cycle.

Comparison of multiscale characteristics of ballast bed under cyclic load and train load
The mean particle acceleration under two loading patterns is shown in Fig. 10.It is observed that the general variation and distribution of particle acceleration are independent from loading pattern.For example, no matter for cyclic or train load, the z-acceleration is larger than x-and y-accelerations; particle accelerations at sub-rail regions (R 2 and R 3 ) are much larger than those at other regions, while, quantitatively, the particle acceleration along a direction under cyclic load is lower than that under train load, except the z-acceleration at R 2 and R 3 (which are close).
Figure 11 shows the comparison of contact state under cyclic and train loads.The coordinate number under cyclic load is larger than that under train load, as shown in Fig. 11a.This indicates that more inter-particle contacts form to resist the increase in cyclic load.Owing to the same magnitude of two load patterns, the amplitude of contact force under cyclic load should be smaller.In fact, the amplitude of contact force under cyclic load is indeed smaller than that under train load, as shown in Fig. 11b.
Figure 12 shows the lateral distribution of the amplitudes of δ c and δ f under two loading patterns, which are close for R 1 , R 2 , R 4 and R 5 , while, for the sub-rail region (R 3 ), amp c under train load is much larger than that under cyclic load.It is also observed that amp f is much larger than amp c , indicating that the directional distribution of contact force is more anisotropic than contact normal.In other words, not only more inter-particle contacts form along vertical direction under vertical load, but also the mean normal force on these contacts is considerably larger than the mean normal force on other contacts.This reflects how granular materials self-organize to resist the vertically applied load.
For a region R i , the sub-ballast stress is calculated as the ratio of the overall contact force on bottom wall within the range of R i to that on the bottom area of R i , where the bottom area refers to the product of the width of R i (0.35 m) and the thickness of model box (0.6 m). Figure 13 shows the comparison of the amplitudes of sub-ballast stresses under cyclic and train loads, which are very close in both variation tendency and magnitude.Besides, the dynamic stiffness (K) of ballast bed is also calculated following the method in [15].Then, we have K = 81 MN/m and K = 87 MN/m for cyclic and train load, respectively.This indicates the equivalent cyclic load performs well in simulating train load from the perspective of macroscopic responses.
Since the above multiscale responses of ballast bed under cyclic load are different in magnitude from those under train load, the deviation is defined as (10) i ( ) = where i ( ) is the deviation of at R i ; refers to any indi- cator describing the dynamic response of ballast bed, here including particle acceleration, mean contact force, coordinate number, anisotropy of contact network, stress amplitude and track stiffness; refers to the obtained at R i under cyclic loading, and tra i refers to the obtained at R i under train load.The positive value of ( ) indicates that using equivalent cyclic load tends to overestimate the real intensity of , while the negative value of ( ) indicates that using equivalent cyclic load tends to underestimate the real intensity of .
Then for each indicator, the deviation is calculated and shown in Fig. 14.It is observed that the deviation of each indicator varies by region.Meanwhile, the maximum deviation of each indicator occurs at different regions.For macroscopic responses (σ z and K), the deviations are below 6%.This indicates that using cyclic load to simulate train load is appropriate in obtaining ballast bed deformation and regional stress, while, for particle-scale responses, the absolute values of deviations are much larger.The maximum deviation of F are about 19.5% (R 3 ), 44.7% (R 2 ), 47.1% ( A x at R 1 ), 39.7% (R 5 ), and 24.1% (R 5 ), respectively.This evidences that the similar macroscopic responses (σ z and K) result from quite different micro-mechanical process.Particle motion and inter-particle contact force has significant influence on abrasion and breakage of ballast particles.Contact force between particles at shoulder region (R 5 ) are important to lateral confinement and stability of ballast bed.Therefore, the deviations of particle-scale responses under cyclic loading and train load are worth paying attention.
It is concluded from Fig. 14 that compared with macroscopic responses, particle-scale responses are more sensitive to the variation of loading type.The inconsistency of multiscale responses to the variation of loading pattern means that the application scale of DEM model should not exceed the validation scale when involving load simplification.Otherwise, it is important to examine multiscale responses to validate the effectiveness of load simplification.
It should be noted that the deviations shown in Fig. 14 are obtained from the numerical test model shown in Fig. 1.For DEM simulations with different individual features, such as model setup, ballast shape, particle parameters, and loading conditions, the deviation of multiscale responses induced by load waveform variation could be different from Fig. 14.
At least two possible reasons account for this inconsistency of multiscale responses to the variation of loading pattern: (1) The two loading waveforms have different impulse effects, which lead to inconsistent effect on multiscale  responses.In this case, the above results obtained from DEM modelling just show the inconsistency.This indicates that the equivalent cyclic load of train load performs well in a specific scale (macro-or micro-scale).( 2) This inconsistency may lie only in DEM modelling; e.g., DEM modelling under two loading patterns has accumulated difference during force-displacement iteration.Experimental study of particle-scale characteristics under two loading patterns is expected to illustrate whether the inconsistency exists only in numerical simulation.Before targeting reasons for the inconsistency of multiscale responses induced by load simplifications, it is necessary to pay some attentions to the influence of such inconsistency on characterization of multiscale dynamic responses of a granular system.

Conclusions
In this numerical study, we have examined the consistency of multiscale responses under train load and the corresponding equivalent cyclic load.For cyclic load and train load with the same magnitude and amplitude, the induced multiscale responses of ballast bed are different not only in variation pattern but also in magnitude.The maximum difference of macroscopic responses under cyclic load and train load is within 6%, while the maximum difference of particle movements under the two loading patterns reaches up to about 47.1%.This indicates an inconsistent effect of dynamic load waveform on macro-and micro-scale responses: compared with macroscopic responses, particle-scale responses are more sensitive to the variation of loading type.However, the reason for this inconsistency is not fully identified yet, which could be more than the difference of load pattern.Since ballast bed is a discrete system with relative small critical length, it is important to keep the conditions adopted in DEM model the same as the practical conditions of a ballast bed.When simplification (not only loading waveform) is involved in DEM modelling, it is suggested to examine the multiscale responses in various conditions to recognize the potential scope where such simplification performs well; or at least the application scale of DEM model should not exceed the validation scale.

Fig. 1
Fig. 1 Model setup in previous laboratory test [5, 27]: a laboratory model test with semi-sleeper; b ballast gradation the DEM model shown in Fig. 2. It is observed from Fig. 3 that the macroscopic responses (sub-ballast stress distribution and dynamic stiffness of ballast bed) in DEM modelling are close to those in laboratory test.This validates the DEM model shown in Fig. 2.

Fig. 2
Fig. 2 The process of reproducing ballast particles and DEM model of ballast bed

Fig. 3 Fig. 4
Fig. 3 Comparison of the dynamic responses of ballast bed in numerical and laboratory tests: a stress distribution beneath ballast bed; b dynamic stiffness K of ballast bed

Fig. 6
Fig.6 Translational acceleration of the probe particle 10 cm beneath the sleeper end: a under cyclic load; b under train load.a x , a y , and a z are the translational accelerations in the x-, y-, and z-directions, respectively

Fig. 7
Fig. 7 Time history of contact force: a under cyclic load; b under train load

Fig. 8
Fig. 8 Directional probability density distribution of contact force of R 3 : a under cyclic load; b under train load

Fig. 9 Fig. 10 Fig. 11 Fig. 12 Fig. 13 Fig. 14
Fig. 9 Time history of the anisotropy deviators of contact force: a under cyclic load; b under train load