E3S Web of Conferences 415, 07007 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341507007

A CFD-DEM Based Numerical Investigation of Debris Flow on
Ballasted Railway Track
Yufeng Gong1, and Yu Qian1*
1Department of Civil and Environmental Engineering, University of South Carolina, Columbia, South Carolina 29208, USA

Abstract 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. Afterthe flash of a debris flow, the fine particles can be retained in the ballast layer and significantly reduce trackdrainage, 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, thosedebris flow introduced track issues have not been thoroughly investigated. This study presents a numericalinvestigation 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 anddiscrete 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.

1 Introduction
The debris flow threads the safety of railway
engineering in mountainous areas [1-4]. Not only thedirect impact of the debris flow onto the train but also
the impairing on the structure of the railway, especiallythe substructure. This may cause service disruption andeven 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 toexcessive settlement, degraded track modulus, and poor
drainage.
Traditionally, the finite element method (FEM) andcomputational fluid dynamics (CFD) are popular
numerical methods to investigate the fouling effects on
the ballast [6-10]. However, these methods use theporous media to represent granular materials and ignore
the movement of fine particles within the ballast. For abetter maintenance strategy, it is critical to understandhow and where fines would initiate settlement andcumulate in the ballast under a certain driven force, likerainfall seepage, train-induced vibration, and debris
flow flashing.
The coupled computational fluid dynamics anddiscrete element method (CFD-DEM) has recently beenused widely in the geotechnical and geological areas[11-13]. With coupled CFD-DEM, the motion ofparticles 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 offine particles in the ballast layer can be simulated
without the assumption of a fully saturated track.
Moreover, the quantified fine particle accumulation andwater content distribution throughout the ballast layer
can assist in decision-making for railway maintenance.

2 Methodology
2.1 Computational domain
The computational domain in this study contains atypical railway cross-section in North America[9,10,14,15]. For a newly constructed ballast layer, thesize of aggregates ranges from 0.019 m to 0.051 m but
is mostly uniformly graded. This study simplifies thecomplexity by assuming all the ballast aggregates sharethe 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 ina staggered formation, forming a porous region with aporosity 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* Corresponding author: yuqian@sc.edu© The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0(https://creativecommons.org/licenses/by/4.0/).

E3S Web of Conferences 415, 07007 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341507007

Fig. 1. Computational domain of fluid zone populated with the ballast

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:

computational resource and to include one full ballast.
The top, bottom, and ballast outlines are set to be wallboundaries. A vertical slope with size 1.5m x 2 m filledwith backfillings is set to fail after saturation, then thedebris 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 thedebris flow are tracked and used to quantify theperformance of the rail track.

𝑚𝑖𝐼𝑖

2.2 Governing equations of the CFD-DEMmethod𝜕𝛼𝑐+ ∇ ∙ (𝛼𝑐 𝐔𝑐 ) = 0𝜕𝑡d𝛚𝑖𝑐= Σ𝐌𝑖𝑗𝑗d𝑡

(4)
(5)

where 𝑚𝑖 , 𝐼𝑖 , 𝐮𝑖 , and 𝛚𝑖 are a particle's mass, moment𝑐of inertia, velocity, and angular velocity. 𝐅𝑖𝑗𝑐 and 𝐌𝑖𝑗
arethe contact force and contact moment within particles.
Two forces are considered in this research to accountfor 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:

The CFD-DEM method involves two different types ofphases: the continuous phase-fluid and the discrete
phase-solid. In this study, the continuous phase iscomposed of two different fluids: water and air, and its
governing equations are formulated as follows:𝜕(𝛼𝑐 𝜌𝑐 𝐔c )+ ∇ ∙ (𝛼𝑐 𝜌𝑐 𝐔c 𝐔c ) =𝜕𝑡∇ ∙ (𝛼𝑐 𝛕) − 𝐠 ∙ 𝐡∇𝜌𝑐 − ∇𝑝𝑟𝑔ℎ + 𝜎𝜅∇𝛼1 − 𝐅𝑠d𝐮𝑖∇𝑝= Σ𝐅𝑖𝑗𝑐 + 𝑚𝑖 𝐠 + 𝐟𝑖𝑑 + 𝐟𝑖
𝑗d𝑡

(1 − 𝛼𝑐 )𝜇𝑐𝑚𝑖𝜌𝑐 |𝐮𝑖 − 𝐔𝑐 |[150+ 1.75] , 𝛼𝑐 < 0.82 𝑑2𝜌𝑑𝑖𝛼𝑖𝑐 𝑖(6)𝐟𝑖𝑑 =3 𝑚𝑖 𝜌𝑐 |𝑢𝑖 − 𝑈𝑐 |𝐶𝑑,𝛼
≥0.8𝑐
4 𝜌𝑖
𝛼𝑐2.65 𝑑𝑖{(1)(2)

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 stresstensor, 𝐠 is the gravitational acceleration vector, 𝐡 is
the coordinate vector, 𝑝𝑟𝑔ℎ = 𝑝 − 𝜌𝐠 ∙ 𝐡 is themodified 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 continuousphase, 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 thevolume fraction at the water-air interface, this modified
volume of fluid (VOF) equation is employed:

where 𝑚𝑖 , 𝜌𝑖 , and 𝑑𝑖 are the mass, volume, and diameterof a particle, 𝜇𝑐 is the dynamic viscosity, and 𝐶𝑑 is thedrag coefficient. 𝐶𝑑 depends on the value of the particle
Reynolds number 𝑅𝑒 = 𝛼𝑐 𝜌𝑐 𝑑𝑖 |𝐮𝑖 − 𝐔𝑐 |/𝜇𝑐 and isdetermined by the equation [16]:0.44, 𝑅𝑒 > 1000
𝐶𝑑 = { 24(1 + 0.15𝑅𝑒 0.687 ), 𝑅𝑒 ≤ 1000𝑅𝑒(7)The pressure gradient force is formulated as:∇𝑝𝐟𝑖=−𝑚𝑖∇𝑝𝜌𝑖(8)Detailed parameters employed in modelling thedebris flow are summarized in Table 1.

𝜕(𝛼𝑐 𝛼1 )∇𝛼1+ ∇ ∙ (𝛼𝑐 𝛼1 𝐔c ) + ∇ ∙ (𝛼𝑐 𝛼1 (1 − 𝛼1 )𝑐|𝐔𝑐 |) = 0 (3)|∇𝛼1 |
𝜕𝑡2Table 1. Parameters used in the physical modelling

Items

Properties

Value

E3S Web of Conferences 415, 07007 (2023)
DFHM8Particle

WaterAirhttps://doi.org/10.1051/e3sconf/202341507007

Density, kg/m3
Young’s modulus, GPa
Poisson’s ratio
Restitution coefficient
Rolling friction coefficient
Density, kg/m3
Viscosity, Pa∙s
Density, kg/m3
Viscosity, Pa∙s

2.65e3
70
0.3
0.7
0.15
1e3
1e-06
1
1.48e-05

LFO =References
1.

M. Pastor, J. A. F. Merodo, M. Quecedo, M. I.
Herreros, E. González, P. Mira, Revue française degénie civil 6 (6), 1213 (2002).
2. V. Vilímek, J. Klimeš, J. Vlcko, R. Carreño,Environmental Geology 50 (7), 1041 (2006).
3. A. C. Tillery, M. J. Darr, S. H. Cannon, J. A.
Michael, US Geological Survey Open-File Report
1257 (9) (2011).
4. L. Fan, P. Lehmann, B. McArdell, D. Or,Geomorphology 280, 1 (2017).
5. E. T. Selig, J. M. Waters, Track geotechnology
and substructure management. (1994).
6. E. T. Selig, presented at the Proceedings 6thInternational Heavy Haul Railway Conference,Cape Town, South Africa, 1997 (unpublished).
7. K. Rushton, G. Ghataora, presented at theProceedings of the Institution of Civil EngineersTransport, 2009 (unpublished).
8. N. Tennakoon, B. Indraratna, C. Rujikiatkamjorn,
S. Nimbalkar, T. Neville, Geotechnical Testing
Journal 35 (4), 629 (2012).
9. J. I. Shu, Y. Wang, Y. Qian, J. A. Khan, S. R.
Schmidt, Transportation Geotechnics 28, 100532(2021).
10. Y. Gong, Y. Qian, X. Chen, Transportation
Research Record 0 (0), 03611981221098659(2022).
11. J. Zhao, T. Shan, Powder Technology 239, 248(2013).
12. L. Jing, C. Kwok, Y. Leung, Z. Zhang, L. Dai,
Journal of Geophysical Research: Earth Surface
123 (8), 2004 (2018).
13. X. He, H. Xu, W. Li, D. Sheng, Computers andGeotechnics 128, 103818 (2020).
14. Y. Qian, S. Zhu, Y. Wang, D. C. Rizos, presented
at the ASME/IEEE Joint Rail Conference, 2018
15. G. Cui, M. Liu, W. Dai, Y. Gan, International
Journal of Multiphase Flow 114, 19 (2019).
16. C. T. Crowe, J. D. Schwarzkopf, M. Sommerfeld,
Y. Tsuji, Multiphase flows with droplets andparticles. (CRC press, 2011).

Table 2. Tests matrix

123

Debris Flow
Diameter, m Solidfraction
4.75e-3
0.5
4.75e-3
0.3
4.75e-30.5

44.75e-3

0.3

Ballasted TrackFouling
Void
Profile
Ratio
Clean
0.67
Clean
0.67
Moderately 0.67
FouledModerately 0.67
Fouled3 Preliminary Discussion
According to the definition of fouling index [5], foulingindex (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 ofthese ballast stones is 2.7e3 kg/m3, the density of granite,
the total mass of ballast is:𝑀𝑏 = 1292 × 2700 × 0.02543 = 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 toquantify the FI in the track after the debris flow flashing:
FI =𝑀𝑟𝑁𝑟 𝑑𝑖3=3
𝑀𝑟 + 𝑀𝑏 𝑁𝑟 𝑑𝑖 + 0.041(11)

Also, with the new post-processing mesh, watercontent distribution can be quantified if any water is left
within the ballast layer.

The test matrix listed in Table 2 is performed tostudy the influence of the debris flow on the ballasted
track parametrically and quantitatively.

No.

𝑉𝑟𝐶× 100%
4.7244 × 0.381 × 0.025
30 × 6(10)where Nr is the number of solid particles retained withthe 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 eachcell, a new index-local fouling occupancy (LFO) isdefined to represent the number of particles in every cell:Fig. 2. Illustration for the local fouling occupancy (LFO)

3