A new time-delayed periodic boundary condition for discrete element modelling of railway track under moving wheel loads

A new time-delayed periodic boundary condition (PBC) has been proposed for discrete element modelling (DEM) of periodic structures subject to moving loads such as railway track based on a box test which is normally used as an element testing model. The new proposed time-delayed PBC is approached by predicting forces acting on ghost particles with the consideration of different loading phases for adjacent sleepers whereas a normal PBC simply gives the ghost particles the same contact forces as the original particles. By comparing the sleeper in a single sleeper test with a fixed boundary, a normal periodic boundary and the newly proposed time-delayed PBC (TDPBC), the new TDPBC was found to produce the closest settlement to that of the middle sleeper in a three-sleeper test which was assumed to be free of boundary effects. It appears that the new TDPBC can eliminate the boundary effect more effectively than either a fixed boundary or a normal periodic cell.


Introduction
One major challenge for discrete element modelling (DEM) of practical particulate systems is the low computational efficiency as the real field system normally involves millions of particles. Element testing models are the most commonly used method to reduce the number of modelled particles and thus increase the computational efficiency. For example, the railway track system has been investigated by simulating 'box tests' [1][2][3][4][5][6], cone penetration tests (CPT) and triaxial tests [7][8][9][10][11] to understand the ballast mechanics. However, the artificial boundaries in such simulations can have significant effects on granular movement [12,13]. To eliminate the boundary effects, periodic boundary conditions (PBCs) have been recognized as a useful tool [11,12,[14][15][16][17] for some models. PBCs are normally used when the whole particle system can be represented by repeated and identical representative unit cells along the whole space [18]. The PBC was initially applied to rectangular samples [10,11]. Figure 1a illustrates the principle of the rectangular PBC which simply allows the particles located at one of the periodic boundaries to interact with the particles located at the other, opposite periodic boundary, which is achieved by duplicating the particles in the simulation cell itself.. Cui et al. (14) proposed a circumferential PBC adapted to the axisymmetric samples where the periodic domain was represented as a sector model. The principle is similar and illustrated as Fig. 1b: the particles next to one boundary are duplicated symmetrically about the symmetry axis and are imposed against the other boundary.
However, the above PBCs are not appropriate for nonsymmetrical loads, such as the problem of loads travelling along a railway track. In this case, as Fig. 2 illustrates, the sleepers transfer load from the rail to the ballast, and the instantaneous loads transmitted from the sleepers vary based on their distances to the moving wheel. Typically, in both experimental tests and in DEM simulations, this loading regime is implemented by applying out-of-phase sinusoidal loads to the sleepers [19][20][21], which is demonstrated in \* MERGEFORMAT Fig. 3. Sub-figure (a) shows (as an example) that when the wheel is directly above a sleeper, this sleeper has the maximum load, whilst the adjacent sleepers carry approximately ½ load [22][23][24]. \* MERGEFORMAT Fig. 3b shows the non-uniform applied loads for when the wheel is between sleepers.
A DEM model including only one sleeper and the surrounding ballast (labelled simulation domain in Fig. 2) was used as the simulation cell in this study and is referred to as a 'box test'. The box test has been widely used to study railway track loading both experimentally and numerically [2,4,[25][26][27][28][29][30][31][32][33], due to time/cost efficiency. In this case, the modelled representative element is not completely repeatable, as adjacent sleepers are subjected to different phases of the cyclic moving loads ( \* MERGEFORMAT Fig. 3b), and thus the particles located at the boundary (e.g. particle A in Fig. 1a) and its duplicated particle at the other side (particle A') would be subjected to different stress conditions by sleepers at the different loading phases. Although the dynamics of structures under periodic (cyclic) moving loads have been studied by quite a number of researchers, most of them were employing mathematical methods such as a Fourier-series technique [20,21,34], the Floquet theorem [35,36] and Floquet decomposition [37]. To the authors' knowledge, there is no published DEM study investigating the application of PBC under moving loads. This study aims to propose a novel 'time-delayed' PBC for DEM modelling of a non-symmetrical periodic structure under moving loads.
Two box test configurations were conducted in this study, one with a single sleeper and the other with three sleepers (in which loads can be applied out of phase to the three individual sleepers). The mechanical behaviour of the middle sleeper in the three-sleeper box test was regarded as being sufficiently far from the boundaries for the boundary effect to be negligible; the adjacent sleepers create the necessary conditions for the central sleeper to be subject to the correct loading conditions. The behaviour of this test, specifically that of the central sleeper, is compared to that of the single sleeper box test, using various boundary conditions. Three alternative boundary conditions were used: fixed boundary, normal PBC and the new time-delayed PBC. By comparing the behaviour of the central sleeper in the three-sleeper test to the single sleeper

Box test conditions
The box test simulates the ballast-sleeper interaction occurring under the rail seat of a track (Fig. 4). Previous studies [2,[38][39][40] have shown the apparatus used in this study to be a useful method of simulating realistic track loads on railway ballast in order to assess the effects of specific track variables. It simply consists of loading cyclically a section of sleeper which is embedded into ballast confined in a 0.3 × 0.7x0.45 m box (Fig. 5a). The sleeper section is loaded vertically with a 3 Hz cyclic load oscillating between 3 and 40kN. In the present numerical approach a three-sleeper case has also been simulated using a box three times as wide (Fig. 5b), similar to the Nottingham Railway Test Facility [19] as shown in Fig. 6. In this case the sleepers are loaded sequentially with 90 degrees phase angle delay along the axis of the propagation of the load which are similar to the conditions for the Railway Test Facility, as shown in \* MERGEFOR-MAT Fig. 7. To measure the effect of the boundary conditions on the settlement of the sleeper, results of the single sleeper box test are compared with the results obtained with the middle sleeper of the three-sleeper test.

DEM model
The commercial DEM code PFC3D [41] was used in this study. DEM considers granular materials like ballast as an assembly of objects interacting through a contact law. A Hertz-Mindlin [41] contact law characterised by a Youngs modulus E and a Poisson ratio v has been used to calculate, from the overlap between two objects, normal and shear elastic forces. The ratio between the shear and normal contact force components is limited using a Mohr-Coulomb sliding criterion characterised by a friction coefficient f. This is the contact model embedded in  the DEM code PFC 3D used in this analysis. The acceleration of each object is deduced from the contact forces acting on it using Newton's second law, and then its velocity and displacement by integration using a timestepping scheme. The displacement of the objects then leads to new overlaps and new contact forces creating a complete DEM calculation cycle. The displacement of the objects is then determined as a function of time (number of timesteps x time increment). A non-viscous damping force reducing the acceleration of the particles by 70% is introduced into the model to avoid non-physical oscillation within the assembly of objects and allow further dissipation of energy in addition to friction at contacts [41].
The particles are modelled as 'clumps' with a single, realistic shape as scanned by Li et al. [42,43] as \* MERGEFORMAT Fig. 8 shows. For a given scanned particle surface, PFC3D is able to create a 'clump' of the same shape by using the algorithm of Taghavi [44] (see also Li and McDowell [2]). The clumps used in the simulations have a particle size distribution typical of ballast used by Network Rail UK [45]. The mechanical parameters of the simulated ballast particles were determined in previous studies [2] and are summarized in Table 1. The same contact model operates for the walls. Figure 9 shows the generated DEM samples for the two test configurations. The top sample was generated with a simple deposition method; details are given in Li and McDowell [2]. The bottom sample is generated by deleting the two entire outer "boxes" (i.e. outer thirds of the sample) including the boundaries, particles and sleepers of top sample and applying alternative boundary conditions (fixed boundary, normal PBC; new time-delayed PBC) to the left and right hand sides of the sample.

Boundary effect for fixed boundary and normal periodic boundary
Firstly, fixed boundaries and standard (no time-delay) PBCs were applied to single sleeper tests. The fixed boundary here means rigid walls with the same parameters as listed in Table 1. \* MERGEFORMAT Fig. 10 shows the sample with standard PBCs-which has periodic boundaries on either side. The light grey particles are the standard particles and constitute the bulk of the sample. The red particles are the 'mirrored' or duplicated particles, which will be called 'ghost particles' in this paper. Only the 2 lateral boundaries shown in the figure were used as periodic boundaries, as this is the direction in which traffic loads move along (not into/out of the page). The standard PBCs were applied according to the following procedure. The position and surfaces of all particles were continuously monitored every timestep. When any particle surface comes in to contact with (touches or overlaps) the lateral boundaries/planes (e.g. Particles A and B in \* MERGEFORMAT Fig. 10), a corresponding particle is created at the opposite side of the box (e.g. Particle A' and B'). This duplicated particle has an identical shape and is created with the same orientation. This ghost particle therefore provides a physical boundary at the opposite side of the box, preventing the internal normal particles from escaping the sample.
The net force acting on the ghost particle has the same magnitude and direction as those acting on the original particle. The ghost particle also has an identical velocity to the original particle. For position relocation, it was treated as a shift along the periodic direction. For the sample used in this study, with a coordination system shown in \* MERGEFOR-MAT Fig. 10, the coordinates of ghost particle (x g , y g , z g ) for the original particle (x o , y o, z o ) could be represented as: where L is the length of the sample in the periodic direction.
\* MERGEFORMAT Fig. 11 shows the deflections of the sleepers as a function of time for the cases of the central sleeper in three-sleeper test, the single sleeper box test with the fixed boundary and also with normal PBCs. The central sleeper in the three-sleeper test, which is considered as being without boundary effects in this study, clearly presents a much larger settlement than the single sleeper box test with the fixed boundaries. That is to say, the application of simple fixed boundaries does cause a boundary effect-clearly reducing the sleeper settlement. The settlement curve for standard PBC lies in the middle, which means although this boundary effect has been attenuated by changing the fixed boundary to a normal PBC, the PBC is not able to eliminate the boundary effect. (1)

The time-delayed periodic boundary (TDPBC)
The reason why the normal periodic boundary cannot eliminate the boundary effect is that it assumes the particles moving outside the boundary and the corresponding reintroduced particles have the same mechanical behaviour (velocity and net contact force), which means it assumes the three sleepers in this study are always being loaded in phase, as \* MERGEFORMAT Fig. 12a illustrates. However, with a moving load, the actual loads on the three sleepers are out of phase as \* MERGEFORMAT Fig. 12b shows. Therefore, to get a more compatible PBC, the mechanical condition applied on the ghost particles should consider the change of loading phase for adjacent periodic elements. The investigations started from the analysis of the contact force chains on the single sleeper test with fixed boundary. \* MERGEFORMAT Fig. 13 shows the contact force chains at four different loading phases of the 10 th  (19) loading cycle at the same scale (where thickness represents force magnitude). It can be seen for all the four time points in the cycle, the shape of contact force chains are quite similar; the major difference is the magnitude of contact force which changes with the sleeper load: the contact forces increase with increasing sleeper loads and should be proportional to and in phase with the sleeper load. This fact was then used to explore the possibility of predicting the contact forces based on the known sleeper loads. \* MERGEFORMAT Fig. 14 plots both the average magnitude of the contact forces for the particles in contact with both the left and right boundaries combined and the sleeper load as a function of number of loading cycles. \* MERGEFORMAT Fig. 14 clearly shows that the average magnitude of the boundary contact forces is approximately proportional to and in phase with the sleeper loads.
If the average magnitude of contact forces F Ave on each boundary is assumed to be a sinusoidal function of loading time during any single loading cycle as observed from \* MERGEFORMAT Fig. 14, the average magnitude of contact forces for both boundaries combined F Ave_current at any time point would be expressed as: where Φ is the phase angle for this time point, F Ave_max and F Ave_min are the average magnitudes of the boundary contact forces at the maximum and minimum sleeper loads during the sleeper loading cycle.
In DEM simulations, for the value of current sleeper load during the current cycle, Φ is a known value, and the F Ave_current is easily calculated directly from the DEM model. F Ave_min is also a known value as it is the first data point for each sinusoidal cycle, so F Ave_max is unknown for the first half of the loading cycle where the sleeper load hasn't reached the maximum point. In this case, F Ave_max could be estimated by solving Eq. 2: In fact, to allow for the fact that the average boundary contact force and maximum average contact force will not be precisely proportional to current and maximum sleeper loads respectively, the current value of F Ave_current for both boundaries combined is updated each timestep and then substituted into Eq. (3) to get the corresponding F Ave_max . Hence, for the known minimum Average contact force for the current cycle and phase angle together with the current average, the average magnitude of boundary contact force as a sinusoidal function of time for the whole loading cycle can be predicted, following Eq. (2) for different phase angles.
The traffic loading is travelling from left to right. So given the calculated boundary contact forces for the fixed wall case above, for the case of a periodic cell which accounts for the travelling load and the out of phase loading of what would    (4) For the rightmost ghost particles the average magnitude of contact forces is therefore: Knowing the average boundary contact forces, it is possible therefore to calculate a current value of contact force for each ghost particle individually: where f detected represents the current total detected contact force acting on the corresponding boundary particle (e.g. particles A and B) and f left_ghost and f right_ghost are the calculated and applied forces on the ghost particles B' and A' at the left and right hand sides of the sample respectively.
By applying f left_ghost and f right_ghost to the ghost particles, the modified periodic boundary was called a 'time-delayed periodic boundary'. It is also noted here the velocities and positions are assumed to be those calculated as if there was a normal PBC. This time-delayed PBC was then applied to three different samples (different initial particle locations, as would be the case in a laboratory sample) to verify its effectiveness. \* MERGEFORMAT Fig. 15   It clearly shows the time-delayed periodic boundary presents the closest settlement curve when compared to the threesleeper test. Although there is still a slight gap between the three-sleeper test and the time-delayed boundary, the new time-delayed PBC predicts the sleeper settlement much better than the normal PBC and the result seems to be repeatable for different initial samples. Therefore, the time-delayed PBC is more useful when considering an element test with a moving wheel load.

Conclusions
A novel time-delayed periodic boundary has been proposed for DEM modelling cyclic loading in an element test when the particulate system is subject to travelling loads such as in a railway track. The investigations were based on a box test which is normally used as an element testing model for railway track in DEM studies and also in some laboratory tests. The proposed time-delayed PBC was conducted by predicting the forces acting on ghost particles with the consideration of different loading phases for adjacent sleepers while the normal PBC simply gives the ghost particles the same contact forces as the original particles at Single sleeper with rigid boundary the boundaries. The settlement of the middle sleeper in a three-sleeper test was assumed to give the result of settlement without boundary effects (for that sleeper) and it was compared with the sleeper settlements in a single sleeper test with a fixed boundary, a normal periodic boundary and the newly proposed time-delayed PBC. The comparison shows the time-delayed PBC increased the accuracy of predicting sleeper settlement compared to the normal PBC or the fixed boundary, which proves that this time-delayed PBC could effectively reduce the boundary effect in DEM modelling of box tests. This provides a new option for researchers to model other similar periodic structures with travelling loads using DEM.