E3S Web of Conferences 415, 01007 (2023)
DFHM8

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

Numerical investigation of particle dynamic behaviours ingeophysical flows considering solid-fluid interaction
Jun Fang1, Yifei Cui1*, Xinyue Li2, and Hui Tang3
1State Key Laboratory of Hydroscience and Engineering, Tsinghua University, China
2Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai, China
3Section 4.7: Earth Surface Process Modelling, German Research Centre for Geosciences (GFZ), Potsdam 14473, Germany

Abstract. Solid-fluid interaction vitally influences the flow dynamics of particles in a geophysical flow. A
coupled computational fluid dynamics and discrete element method (CFD-DEM) is used in this study tomodel multiphase geophysical flow as a mixture of fluid and solid phases. The two non-Newtonian fluids(i.e., Bingham and Hershcel-Bulkley fluids) and water mixed with particles are considered in the simulation,
while dry granular flow with the same volume is simulated as a control test. Results revealed that the solidfluid interaction heavily governs the particle dynamic behaviours. Specifically, compared to dry case,particles in three multiphase cases are characterized by larger flow mobility and greater shear rate whilesmaller basal normal force. In addition, a power-law distribution with a crossover to a generalized Pareto
Distribution is recommended to fit the distribution of normalized interparticle contact force.

1 Introduction
Geophysical flows contain a wide type of debris flows,rock avalanches, and snow avalanches based on theclassification proposed by Hungr et al. [1]. They canbehave like a fluid and move far away from their initial
positions with high velocities. Over the past few
decades, there have been considerable numerical andexperimental investigations to reveal key mechanisms ingeophysical flows during their transportation stage,including size segregation, particle breakage andentrainment. Most of these studies have consideredgeophysical flows into either pure fluids or dry particles,with only a small amount in water-particle mixtures [25]. However, it is well known that the solid and fluid
phases vitally influence geophysical flow dynamics [6],and gross simplification will inhibit a comprehensive
understanding of two-phase interaction.
In this study, the focus is placed on the effect of solidfluid interaction on particle mechanical behaviour in a
mixture composed of granular and fluids, where a
Newtonian (water) and two non-Newtonian fluids(Bingham and Herschel-Bulkley fluids) are adopted tomodel the liquid phase, corresponding to Case W-P,
Case B-P, Case HB-P, respectively. In reality, fluids in
geophysical flows, such as slurry, are commonlycomplicated in terms of rheological properties andexhibit non-Newtonian behaviour. Considering nonNewtonian fluids provide better explorations of more
complex natural flows than Newtonian fluids like water
only. In addition, the dry granular flow (Case P) with thesame volume is simulated as a control test.

2 Methodology and model setup
2.1 Methodology
The solid particles are regarded as a discrete phase thatis modelled by the discrete element method (DEM),
while the fluid is a continuous system simulated by thecomputational fluid dynamics method (CFD). It is
assumed that the translational and rotational motions ofeach particle are controlled by Newton’s equations ofmotion, and the continuous fluid system can bedescribed by locally averaged Navier-stokes equation.
The key to the coupling between the CFD and DEM is
to properly consider the particle-fluid interaction forces,including drag force and buoyancy force. More details
about CFD-DEM can be found in the literature [5, 7].
2.2 Model setup
Fig. 1 shows the schematic of the CFD-DEM model
simulating fluid and particle mixture flow collapse onto
an inclined channel. A gravity-driven mixture sample (H= 0.3 m, L = 0.5 m, W = 0.3 m) composed of viscous
fluid and particles is first produced in a pre-defined area
before being released to flow along the flume base toimpact the rigid barrier. Note that for dry granular flow,
only particles are generated in the pre-defined area.
After the packing is stable, the flume will be rotated by20°. Then the door will be removed, and the particles
will flow along the flume base. In the initial state, formixture flows, it is assumed that only the fluid cells in* Corresponding author: yifeicui@mail.tsinghua.edu.cn© 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, 01007 (2023)
DFHM8

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

the initial deposition domain contain water while the restof the CFD fluid cells is air. The DEM time step isadopted based on the criteria of the Rayleigh time[8],whilst the CFD time step is determined as 100 times
larger than that of the DEM to ensure sufficient accuracy
and efficiency.
As summarized in Table 1, besides water, there are two
non-Newtonian fluids (Bingham and Herschel-Bulkelyfluids) adopted to simulate the complex fluid in
geophysical flow. The particle parameters have also
been listed in Table 1. Interested readers can find thedetails of benchmarks of particle and fluid properties inliterature (e.g., Li and Zhao [9], Fang et al. [5]).

3 Results
3.1 Flow regimes
Fig. 2a shows the computed flow dynamics of dry
granular flow (Case P) at three typical instants. It can beseen that, after releasing, the right portion of the particle
column is mobilized by gravity with a wedge front,
while the left still remains static. During the flow
process, the depositing height reduces, and moving
particles become dominant. For comparison, the flow
evolutions of the water and particle mixture flow (CaseW-P) are presented in Fig. 2b. Unlike dry granular flow,
water and particle mixture flow exhibits a larger
spreading velocity, and a longer flow distance can beobserved simultaneously. In addition, due to solid-fluid
interaction, the collisions between particles areweakened, and the phenomenon where dispersedparticles move ahead of the flow front and detach from
subsequent particles is not obvious in Case W-P.

Fig. 1. Model setup for CFD-DEM simulation(a)

Table 1. Input parameters adopted for the CFD-DEM model
Material
Parameters
ValuestypeParticle
Barrier Young’s modulus (Pa)
3.2×109
Particle Young’s modulus
6×1010(Pa)
Coefficient of restitution
0.78
Particle density (Kg/m³)
2550
Gravitational acceleration
9.81(m/s²)
Inter-element friction
0.48
coefficientInterface-element friction
0.4
coefficientParticle diameter (mm)10Particle Poisson ratio
0.25Barrier Poisson ratio
0.35
Air
Viscosity (Pa·s)
1.48×10-5
Density (Kg/m³)
1.0
Water
Viscosity (Pa·s)
1×10-3
Density (Kg/m³)
1000
BinghamDensity ρ (Kg/m3)
1400
fluid
Consistency index k (Pa·sn)
4
Flow index n
1
Yield stress τ0 (Pa)
2.1
HerschelDensity ρ (Kg/m3)
1000
Bulkley
Consistency index k (Pa·sn)
4.279
fluid
Flow index n0.479
Yield stress τ0 (Pa)30.002
Simulation
DEM time step (s)5×10-7control
CFD time step (s)5×10-5
Simulated time interval
100(b)
Fig. 2. Flow process for Case P (a) and Case W-P (b)It is apparent the solid-fluid interaction can influence theflow dynamics. We further analyze the vertical centre ofparticle position z, as shown in Fig. 3. The verticalcentres of particle position of three mixture cases arelower than that of Case P, which signifies that theparticles in the three mixture cases collapse faster. Thisis caused by the solid-fluid interaction, where the fluid
can impose a driving force to enhance the flow dynamicsof the particles and a buoyant force to reduce theinterparticle collision [10]. Furthermore, it should be
noted that there are also some differences among thethree mixture flows, where vertical centres of particle
position in Case B-P slightly deviate from the other two
mixture flow (Case W-P and Case HB-P). The higher
density of Bingham fluid may explain this phenomenon
compared to the other two fluids in the mixture cases,
which can provide greater dragging and buoyancy force
on the particles [9].

2E3S Web of Conferences 415, 01007 (2023)
DFHM8

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

Fig. 3. Time evolution of the vertical centre of particle
𝟏
position z, where 𝐳 = ∑𝒏𝒊=𝟏 𝒛𝒊 , where n and zi denotes𝒏particle number and vertical centre height of ith particle,
respectively. The normalizer of time is calculated as 𝐓 =√𝒈𝑯, where 𝑯 is the flow thickness and 𝒈 is gravitational
acceleration(a)During the transportation process, shearing is developed
along with the flowing thickness and the traveling
direction. Fig. 4 shows the computed velocity
distribution along the normal and traveling direction for
four geophysical flows at t=0.6s, and the corresponding
flow profile is also presented. Different symbols areadopted to represent the particle velocities at thedifferent positions along flow thickness, where thedistance of particles is scaled away from the rear (L* =0) to the front (L* = 1). The reciprocal values of the fitted
curve slopes represent shear rates. The results show thatthe shear rate in the particles increases from the rear to
the front along the flow direction for all cases. An
apparent increase in the shear rate near the bottom isobserved, which indicates that there is a boundary layerdeveloping. The boundary effect of the flume base may
be responsible for this.
In this study, it can be found the velocity profiles aremainly in power-law distribution:
H ~ (v)k(b)(1)

The nature seems to be independent on the properties of
the incoming flow. Compared with the three mixture
cases, the dry granular flow (Case P) has a smaller kvalue along the flow direction, which implies aweakened shearing of the granular materials. Especiallyat the location of L* = 0.1, the fitting curve of Case Pshows upper convex with k smaller than 1.0, while thelower convex tendency is observed in other fitting
curves with k larger than 1.0. Because of different fluid
properties, including density and viscosity, the minor
difference can also be found among the three mixture,
where Case B-P has a maximum shear rate.

(c)

(d)
Fig. 4. Computed velocity distribution along the height andtraveling direction for four geophysical flows at t = 0.6s: (a)Case W-P; (b) Case HB-P; (c) Case B-P; (d) Case P3

E3S Web of Conferences 415, 01007 (2023)
DFHM8

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

3.2 Basal normal force
Basal normal force Fn is caused by the impact ofparticles, which is composed of single free particle
impact 𝐹𝑛𝑠 and multiple particle impact 𝐹𝑛𝑐 . The single
free particle impact 𝐹𝑛𝑠 is caused by particles free fromsurrounding particles. And the multiple particle impact
𝐹𝑛𝑐 results from the impact of particles that also have acollision with other particles with a force chain effect[11]. Fig. 5 shows the development of the basal normal
force of four geophysical flows, normalized by the total
weight of particles. And the ratio of 𝐹𝑛𝑠 and Fn is also
added in the figure.
As shown, at the initial stage, the basal normal force Fn
has a small fluctuation, then reaches a static state. Afterstabilization, the magnitude of Fn of four cases evolves
in the following order: Case P > Case W-P > Case HBP > Case B-P. The differences are mainly attributable tothe collective effect of the density and viscosity of thefluid, which further influences the solid-fluid
interaction. In addition, it can be found that at thebeginning, the basal normal force generated by multiple
particle impact 𝐹𝑛𝑐 is the main component, thengradually decreases during the flowing process. The
corresponding impact load of a single free particle 𝐹𝑛𝑠
plays a more important role, with 𝐹𝑛𝑠 /Fn approaching
1.0. This phenomenon indicates that more particles tend
to be dispersed and free from surrounding particles in
the flow evolution. But for dry granular flow (Case P),
it would take longer for particles to become dispersed.

(d)Fig. 5. Computed velocity distribution along the height andtraveling direction for four geophysical flows: (a) Case W-P;(b) Case HB-P; (c) Case B-P; (d) Case P

3.3 Interparticle contact force
Fig. 6 shows the semilogarithmic plots of probability
density function P(f) of the interparticle contact forces N
in the four cases at three typical moments. The contact
force N is normalized by the average contact force <N>.
According to Radjai et al. (1996) [12], the data for forces
lower than the mean have a power-law distribution(Power), while the rest data can be well fitted by an
exponential decay (Exp). However, except at the initialinstants, the force distributions above the mean deviate
from this behaviour, which decays more slowly than anexponential. A generalized Pareto Distribution (GPD) isintroduced to fit the pdfs of contact forces larger thanthe mean [13], which shows a better performance thanexponential behaviour. We conclude that thedistribution of normalized contact forces is independentof our incoming flow properties and can beapproximated by a power-law decay with a crossover toa GPD distribution:
𝑃(𝑓) = 𝑎 ∗ (𝑁/< 𝑁 >)𝑏{1𝑃(𝑓) = ( ) (1 + 𝐵
𝐴,𝑁 < < 𝑁 >1(𝑁/<𝑁>)−𝐶 (−1−𝐵)
𝐴),𝑁 > < 𝑁 >(2)where a, b, A, B and C are constant parameters. The
fitting lines have been reported in Fig. 6 and the valuesof parameters a, b, A, B and C are also indicated.

(a)

(b)

(a)

(c)

4

E3S Web of Conferences 415, 01007 (2023)
DFHM8

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

References
1.

O. Hungr, S. Leroueil, L Picarelli, Landslides 11(2),167–194 (2014)
2. N. Balmforth, R. Kerswell, J. Fluid Mech. 538,399–428 (2005)
3. Z. Lai, L. E. Vallejo, W. Zhou, G. Ma, J. M. Espitia,
B. Caicedo, X. Chang, Geophys. Res. Lett., 44(24),12–181 (2017)
4. A. Leonardi, G. R. Goodwin, M. Pirulli, Acta
Geotech. 14(6), 1949–1963 (2019)
5. J. Fang, L. Z. Wang, Y. Hong, J. D. Zhao,Géotechnique 72(5), 1–62 (2022)
6. R. M. Iverson, Reviews of Geophysics 35(3), 245296 (1997)
7. J. Zhao, T. Shan, Powder Technol., 239(17), 248258 (2013)
8. E. M. Smuts, D. A. Deglon, C. J. Meyer,
Methodology for CFD-DEM modelling ofparticulate suspension rheology, in Proceedings ofthe Ninth International Conference on CFD in theMinerals and Process Industries CSIRO,Melbourne, Australia, 37–53 (2012).
9. X. Li, J. Zhao, Powder Technol. 338, 493–505(2018)
10. J. Fang, Y. F. Cui, X. Y. Li, J. Y. Nie,. Computersand Geotechnics 148, 104790, (2022)
11. J. Estep, J. Dufek. Journal of Geophysical Research
117, F01028 (2012)
12. F. Radjai, M. Jean, J.-J. Moreau, S. Roux, Phys.
Rev. Lett. 77 (2), 274–277 (1996)
13. S.W. McCoy, G. Tucker, J. Kean, J. Coe. J.
Geophys. Res. Earth Surf. 118 (2), 589–602 (2013).

(b)

(c)
Fig. 6. Computed probability density functions of interparticle normalised contact force: (a) t=0.0s, a=3.31, b=-0.67,
A=1.12, B=-0.05, C=2.53; (b) t=0.6s, a=3.17, b=-0.85,
A=1.76, B=0.57, C=2.35; (c) t=1.0s, a=3.27, b=-0.79,
A=1.90, B=0.50, C=2.58

4 Conclusions
The effect of solid-fluid interaction on particle
behaviours in geophysical flows has been examinedusing a coupled CFD-DEM model. It was found that thesolid-fluid interaction has a remarkable influence on theflow mobility and flow regimes of particles ingeophysical flows. Compared with the dry granular
flow, particles in three mixture cases have a larger flow
dynamic and an enhanced shear rate along flow
direction and flow thickness. In addition, thegeophysical fluid can accelerate the dispersion ofparticles, and the impact of single dispersed particles
would become the main source of basal normal force. A
combination of power-law decay and GPD is suggested
to fit the probability density function of the normalizedinterparticle contact force for four cases.
This study is funded by the National Natural Science
Foundation of China (42120104002, 41941019), the Second
Tibetan Plateau Scientific Expedition and Research Program(2019QZKK0903), and the Research Fund Program of theState Key Laboratory of Hydroscience and Engineering (2021KY-04).

5