A CFD-DEM Based Numerical Investigation of Debris Flow on Ballasted Railway Track

Debris flow occurring in mountainous areas can cause issues to railway tracks. Debris flow may cause large track deformation and even track breakage and introduce server ballast fouling afterwards. After the flash of a debris flow, the fine particles can be retained in the ballast layer and significantly reduce track drainage, leading to lower bearing capacity and a higher risk of track lateral stability problems. Moreover, these solid particles may deposit on the railway surface and endanger the train directly. Unfortunately, those debris flow introduced track issues have not been thoroughly investigated. This study presents a numerical investigation of the impact of the debris flow on the railway track. Various factors governing the debris flow are considered, including particle size and solid fraction. Besides, those factors affecting the ballast are also discussed, such as fouling condition and initial void ratio. A coupled computational fluid dynamics and discrete element method (CFD-DEM) approach is developed to capture the interactions between particles/particles, water/air, and particles/fluid. The results from this study may help the railway to improve track resilience before the debris flow and to improve maintenance strategy after the debris flow flashing.


Introduction
The debris flow threads the safety of railway engineering in mountainous areas [1][2][3][4]. Not only the direct impact of the debris flow onto the train but also the impairing on the structure of the railway, especially the substructure. This may cause service disruption and even derailments and fatalities.
When fine particles in debris flow are detained by the ballast aggregates, they fill the voids and cause ballast fouling [5]. Fouled ballast often leads to excessive settlement, degraded track modulus, and poor drainage.
Traditionally, the finite element method (FEM) and computational fluid dynamics (CFD) are popular numerical methods to investigate the fouling effects on the ballast [6][7][8][9][10]. However, these methods use the porous media to represent granular materials and ignore the movement of fine particles within the ballast. For a better maintenance strategy, it is critical to understand how and where fines would initiate settlement and cumulate in the ballast under a certain driven force, like rainfall seepage, train-induced vibration, and debris flow flashing.
The coupled computational fluid dynamics and discrete element method (CFD-DEM) has recently been used widely in the geotechnical and geological areas [11][12][13]. With coupled CFD-DEM, the motion of particles can be captured, and the influence of particles on the fluid can also be simulated. This study employs an improved CFD-DEM model by including two different fluids, water and air. Hence the migration of fine particles in the ballast layer can be simulated without the assumption of a fully saturated track. Moreover, the quantified fine particle accumulation and water content distribution throughout the ballast layer can assist in decision-making for railway maintenance.

Computational domain
The computational domain in this study contains a typical railway cross-section in North America [9,10,14,15]. For a newly constructed ballast layer, the size of aggregates ranges from 0.019 m to 0.051 m but is mostly uniformly graded. This study simplifies the complexity by assuming all the ballast aggregates share the same cubical shape having a side of 0.025 m. All these aggregates are subtracted from the computational domain and left with their outlines as a wall boundary. Totally, 1,292 aggregates are introduced and arranged in a staggered formation, forming a porous region with a porosity of 0.4, similar to the field condition. Fig. 1 shows the sizes of the computational domain and the ballasted track. The thickness of the domain in the longitudinal direction is 0.025 m to save the computational resource and to include one full ballast. The top, bottom, and ballast outlines are set to be wall boundaries. A vertical slope with size 1.5m x 2 m filled with backfillings is set to fail after saturation, then the debris flow forms and advances towards the railway ballast layer. The right is set to be pressure outlet. During the whole simulation, the solid particles from the debris flow are tracked and used to quantify the performance of the rail track.

Governing equations of the CFD-DEM method
The CFD-DEM method involves two different types of phases: the continuous phase-fluid and the discrete phase-solid. In this study, the continuous phase is composed of two different fluids: water and air, and its governing equations are formulated as follows: where is the volume fraction of the continuous phase in every cell, is the density of the fluid mixture, and is the velocity of the continuous phase. On the right-hand side of Eq. (1), is the shear stress tensor, is the gravitational acceleration vector, is the coordinate vector, ℎ = − • is the modified pressure, is the fluid surface tension coefficient, is the curvature at the water-air interface, 1 is the volume fraction of the water in the continuous phase, and is the cell-averaged solid-fluid interaction force. Moreover, 2 is the volume fraction of the air in the continuous phase, correspondingly.
Besides, in order to increase the resolution of the volume fraction at the water-air interface, this modified volume of fluid (VOF) equation is employed: where is a coefficient controlling the interface resolution.
For the discrete phase, Newton's second law of motion is used to formulate the equations for every single particle: where , , , and are a particle's mass, moment of inertia, velocity, and angular velocity. and are the contact force and contact moment within particles.
Two forces are considered in this research to account for the solid-fluid interaction force, the drag force and the pressure gradient force ∇ . Ergun-Wen-Yu drag model is widely used and has the expression: where , , and are the mass, volume, and diameter of a particle, is the dynamic viscosity, and is the drag coefficient.
depends on the value of the particle Reynolds number = | − |/ and is determined by the equation [ The pressure gradient force is formulated as: Detailed parameters employed in modelling the debris flow are summarized in Table 1. The test matrix listed in Table 2 is performed to study the influence of the debris flow on the ballasted track parametrically and quantitatively.

Preliminary Discussion
According to the definition of fouling index [5], fouling index (FI) equals the summation of the accumulative passing of No.4 (4.75 mm) and No.200 (0.074 mm) sieves. Because it is assumed the ballast particles are cubic with the same size. Assuming that the density of these ballast stones is 2.7e3 kg/m 3 , the density of granite, the total mass of ballast is: = 1292 × 2700 × 0.0254 3 = 57.16 kg (9) Since the particle diameter in the debris flow is no larger than the No.4 sieve, the flowing equation can be used to quantify the FI in the track after the debris flow flashing: FI = + = 3 3 + 0.041 (10) where Nr is the number of solid particles retained with the track region. Besides, the distribution of the retained solid particles changes as time evolves. To capture the spatial distribution of particles, the ballast domain is discretized into a mesh of size 60 × 6, as shown in Fig. 2. In each cell, a new index-local fouling occupancy (LFO) is defined to represent the number of particles in every cell: LFO = 4.7244 × 0.381 × 0.025 30 × 6 × 100% (11) Also, with the new post-processing mesh, water content distribution can be quantified if any water is left within the ballast layer.