# A multiple-time-step integration algorithm for particle-resolved simulation
with physical collision time

Zhengping Zhu  Ruifeng Hu  Xiaojing Zheng  Research Center for Applied
Mechanics, School of Mechano-Electronic Engineering, Xidian University, Xi’an
710071, China  Center for Particle-Laden Turbulence, Key Laboratory of
Mechanics on Disaster and Environment in Western China, Ministry of Education,
and Department of Mechanics, College of Civil Engineering and Mechanics,
Lanzhou University, Lanzhou 731000, China

###### Abstract

In this paper, we present a multiple-time-step integration algorithm (MTSA)
for particle collisions in particle-resolved simulations. Since the time step
required for resolving a collision process is much smaller than that for a
fluid flow, the computational cost of the traditional soft-sphere model by
reducing the time step is quite high in particle-resolved simulations. In one
state-of-the-art methodology, collision time is stretched to several times the
flow solver time step for the fluid to adapt to the sudden change in particle
motion. However, the stretched collision time is not physical, the
hydrodynamic force may be severely underestimated during a stretched
collision, and the simulation of sediment transport may be sensitive to the
stretched collision time. The proposed MTSA adopts different time steps to
resolve fluid flow, fluid-particle interaction, and particle collision. We
assessed the MTSA for particle-wall collisions as well as particle-particle
collisions, determined the optimal iteration number in the algorithm, and
obtained excellent agreements with experimental measurements and reference
simulations. The computational cost of the MTSA can be reduced to about one
order of magnitude less than that using the traditional soft-sphere model with
almost the same accuracy. The MTSA was then implemented in a particle-resolved
simulation of sediment transport with thousands of particles. By comparing the
results obtained using the MTSA and a version of the stretching collision time
algorithm similar to Costa et al. (2015), we found that stretching the
collision time reduced particle stiffness, weakened particle entrainment, and
affected some turbulence and particle statistics.

##  1 Introduction

Particle-laden two-phase flows are widely encountered in industrial
engineering and natural environments Balachandar and Eaton (2010); Brandt and
Coletti (2022), such as aeolian sand and dust movement, sediment transport,
and fluidized-bed processes. When simulating particle-laden flows, if the flow
field surrounding particles is fully resolved, it is called a particle-
resolved simulation. A number of different particle-resolved simulation
approaches have been developed, like the immersed boundary (IB) method Mittal
and Iaccarino (2005); Uhlmann (2005); Luo et al. (2007); Breugem (2012); Kempe
and Fröhlich (2012b); Zhou and Fan (2014); Tao et al. (2018); Huang and Tian
(2019); Wang et al. (2019); Verzicco (2023) and the distributed Lagrange
multiplier (DLM)/fictitious domain approach Glowinski et al. (1999, 2001); Yu
and Shao (2007, 2010); Shao et al. (2012). Particle-resolved simulations have
been widely used to obtain high-fidelity simulation data, which is greatly
helpful for developing reduced-order models Bagchi and Balachandar (2003);
Eaton (2009); Homann et al. (2013); Tenneti and Subramaniam (2014); Zhou and
Fan (2015); Akiki et al. (2017); Seyed-Ahmadi and Wachs (2020); Zhang et al.
(2020); Wang et al. (2022); Xia et al. (2022a, b) or studying the mechanism of
the fluid-particle interactions Uhlmann (2008); Kidanemariam and Uhlmann
(2014b); Vreman (2016); Vowinckel et al. (2016); Luo et al. (2017); Peng et
al. (2020); Costa et al. (2021); Yu et al. (2021).

The focus of the present paper is on the collision algorithm for particles in
particle-resolved simulations. Particle collision modeling is essential in
particle-resolved simulations, whether a particle-laden flow is in the dense
regime Kidanemariam and Uhlmann (2014b); Vowinckel et al. (2014); Picano et
al. (2015) or the dilute regime Costa et al. (2020). Soft-sphere model Hertz
(1882); Mindlin and Deresiewicz (1953) is one of the most popular collision
models in particle-resolved simulations when multi-particle contact frequently
occurs. However, the traditional soft-sphere model (_i.e._ , the soft-sphere
model with physical collision time) can be inefficient when applied to
particle-resolved simulations, as the time scale associated with particle
collision is typically much smaller than that of fluid flow, particularly for
rigid particles. This mismatch can lead to expensive computation if using the
same time step to resolve particle collisions and advance fluid flow
simultaneously. On the other hand, the stretching collision time algorithm
(SCTA) Feng et al. (2010); Papista et al. (2011); Kempe and Fröhlich (2012a);
Costa et al. (2015); Biegert et al. (2017); Rettinger and Rüde (2022) was
proposed as a remedy and has been shown to be successful in reproducing the
particle trajectories of particle-wall Joseph et al. (2001); Gondret et al.
(2002); Joseph and Hunt (2004) and particle-particle Yang and Hunt (2006)
collisions. In the SCTA, the collision time $T_{c}$ was stretched to $N$ times
the flow solver time step $\Delta t_{f}$ (_i.e._ , $T_{c}=N\Delta t_{f}$)
since the fluid requires time to adapt to the sudden change in particle
velocity during collisions. By stretching the collision time, the
computational efficiency can be improved compared with the traditional soft-
sphere model. The particle stiffness and damping coefficient are determined by
the desired restitution coefficient $e_{n,d}$ and collision time $T_{c}$.
However, the SCTA can only retain the particle kinematic properties rather
than the mechanical properties as it greatly reduces the physical stiffness.
Most of the SCTA is based on a linear spring-dashpot system, and the stiffness
and damping coefficients were determined by an algebraic expression Izard et
al. (2014); Kidanemariam and Uhlmann (2014b); Costa et al. (2015); Ardekani et
al. (2016). Other implementations are based on the Hertzian contact theory,
which is nonlinear, and the stiffness and damping coefficients were determined
by an optimization procedure Kempe and Fröhlich (2012a); Ray et al. (2015);
Biegert et al. (2017). Rettinger and Rüde (2022) developed an efficient four-
way coupled lattice Boltzmann-discrete element method using the stretching
collision time algorithm. They found that the collision time $T_{c}$ is a
critical parameter of the method, which has a significant effect on the
collision dynamics. The collision time should be stretched as
$T_{c}=2.31D_{p}/c_{s}$ to allow the fluid to adapt to the sudden change in
particle velocity, where $D_{p}$ is particle diameter and $c_{s}$ is the
lattice speed of sound. Recently, different from the SCTA developed based on
the soft-sphere model, Jain et al. (2019) proposed an algorithm based on the
hard-sphere model, and the hydrodynamic forces are incorporated during the
collision using a semi-implicit IB method Tschisgale et al. (2018).

Although the SCTA can significantly reduce the computational cost compared
with the traditional soft-sphere model and reproduce the experimental results
of particle-particle and particle-wall collisions successfully Joseph et al.
(2001); Gondret et al. (2002); Joseph and Hunt (2004); Yang and Hunt (2006),
it has several significant drawbacks. First, particle stiffness in the SCTA is
greatly reduced by stretching the collision time, which is not physical.
Second, the artificially stretched collision time $T_{c}$ should be large
enough to get converged trajectory Rettinger and Rüde (2022). However, the
hydrodynamic force during such an artificially long collision process is
usually neglected and significantly underestimated, which may lead to
erroneous predictions of turbulence statistics Xia et al. (2020) and can be
expected to also affect particle statistics in massive-particle-laden flows,
like sediment transport. To avoid introducing a stretched collision time, Li
et al. (2020) were the first to use the physical collision time $T_{c}$ in
particle-resolved simulations. They adopted single non-uniform grid points
with $O\left(10^{5}\right)$ grid elements to resolve the particle and
lubrication force, yielding excellent agreement with experimental results.
However, such an extra-high grid resolution is almost impossible in cases with
a large number of particles.

The present study proposed a multiple-time-step integration algorithm (MTSA)
to reproduce the experimental results Joseph et al. (2001); Gondret et al.
(2002); Joseph and Hunt (2004); Yang and Hunt (2006) and simulations Costa et
al. (2015) using the soft-sphere model with physical collision time on a
uniformly coarse grid. The grid resolution is about $D_{p}/\Delta x=O(10)$,
where $\Delta x$ is the grid spacing. To overcome the time scale mismatch
problem, three different time steps are employed in the MTSA to resolve fluid
flow, fluid-particle interaction and particle motion during particle
collisions. The pressure Poisson equation is only solved at the fluid flow
time step to reduce computational cost. When time is advanced to the same
instant, the fluid and particle information will be updated and synchronized.
The MTSA can reduce the computational cost of the traditional soft-sphere
model by nearly one order of magnitude and avoids artificially reducing
particle stiffness like the SCTA. Meanwhile, we validated the accuracy and
convergence of the MTSA with experimental data Joseph et al. (2001); Gondret
et al. (2002); Joseph and Hunt (2004); Yang and Hunt (2006), numerical data
Costa et al. (2015), and a reference simulation using the traditional soft-
sphere model with fine time steps. We also presented guidelines for
determining the optimal values of the parameters in the MTSA. The MTSA was
then applied to a sediment transport case with thousands of finite-size
particles. The results of the MTSA were compared for the first time with those
obtained by the SCTA to investigate the effects of particle stiffness in a
particle-resolved simulation of sediment transport.

The structure of the paper is arranged as follows. The governing equations and
physical models for fluid flow and particle motion are given in §2. In §3, the
MTSA is introduced. In §4, the influence of the time substep size is
systematically analyzed and discussed in particle-wall and particle-particle
collisions to yield a universally optimal time substep size. The MTSA is
further applied to a sediment transport case with thousands of finite-size
particles, and the effect of particle stiffness is investigated. Finally, we
present the main conclusions of the paper in §5.

##  2 Models

###  2.1 Governing equations

The particle-laden flow considered here is governed by the Navier-Stokes
equations for the carrier phase and the Newton-Euler equations for the
dispersed particulate phase. The motion of an incompressible, Newtonian fluid
flow is governed by the following continuity and momentum equations:

$\nabla\cdot\bm{u}=0,$ (1) $\frac{\partial\bm{u}}{\partial
t}=-\nabla\cdot\left(\bm{uu}\right)-\frac{1}{\rho_{f}}\nabla
p+\nu_{f}\nabla^{2}\bm{u}+\bm{f},$ (2)

where $\bm{u}$ is the fluid velocity, $p$ is the pressure, $\bm{f}$ is the
volume force imposed to take into account the fluid-particle interaction or
the IB force, $\rho_{f}$ is the fluid density, and $\nu_{f}$ is the fluid
kinematic viscosity.

The translational and angular velocities of a particle are advanced by solving
the Newton-Euler equations, which for a spherical particle are reduced to

$\rho_{p}V_{p}\frac{{\rm d}\bm{u}_{p}}{{\rm d}t}={\oint_{\partial
V}{\bm{\tau}\cdot\bm{n}_{p}dA}}+\left(\rho_{p}-\rho_{f}\right)V_{p}\bm{g}+{\bm{F}_{p,lub}+\bm{F}_{p,col}},$
(3) $I_{p}\frac{{\rm d}\bm{\omega}_{p}}{{\rm d}t}={\oint_{\partial
V}{\bm{r}\times(\bm{\tau}\cdot\bm{n}_{p})dA}}+{\bm{T}_{p,col}},$ (4)

where $\rho_{p}$ is the particle density, $V_{p}$ is the volume of the
particle and equal to $(4/3)\pi R_{p}^{3}$ for a sphere with radius $R_{p}$;
$I_{p}$ is the moment of inertia of the particle and equal to
$(2/5)\rho_{p}V_{p}R_{p}^{2}$ for a spherical particle; $\bm{u}_{p}$ and
$\bm{\omega}_{p}$ are the translational and angular particle velocities,
respectively; $\bm{\tau}=-p\bm{I}+\mu_{f}(\nabla\bm{u}+\nabla\bm{u}^{T})$ is
the hydrodynamic stress tensor (the superscript $T$ indicates the
transposition of a tensor); $\mu_{f}$ is the fluid dynamic viscosity;
$\bm{n}_{p}$ is the outward-pointing unit normal vector at the surface
$\partial V$ of the particle; $\bm{r}$ is the relative position vector of the
particle surface to the particle center; $\bm{g}$ is the gravitational
acceleration; $\bm{F}_{p,lub}$, $\bm{F}_{p,col}$, and $\bm{T}_{p,col}$ are the
total lubrication force, total collision force, and total collision torque
acting on the particle $p$, respectively. The subscript $p$ refers to the
quantities of particle $p$. The carrier and dispersed phases are coupled by
the volume force $\bm{f}$, which is calculated by the multidirect forcing IB
method Luo et al. (2007); Breugem (2012). The governing equations of the
carrier phase are discretized on a uniform staggered Cartesian grid, which is
referred to as the Eulerian grid. A particle is discretized by a fixed grid
attached to the particle surface, which is referred to as the Lagrangian grid.
The quantities between the Eulerian and Lagrangian grid points are transferred
through a regularized Dirac delta function $\delta_{d}$ Roma et al. (1999).

To avoid the difficulty of directly calculating the surface integrals
$\oint_{\partial V}{\bm{\tau}\cdot\bm{n}_{p}dA}$ and $\oint_{\partial
V}{\bm{r}\times(\bm{\tau}\cdot\bm{n}_{p})dA}$, these surface integrals are
converted into volume integrals by means of a momentum balance over the volume
occupied by a particle. The volume force $\bm{f}$ is transferred from the
Eulerian grid to the Lagrangian grid. Based on the above operations, the
equations of particle motion (3) and (4) can be rewritten as Uhlmann (2005);
Kempe and Fröhlich (2012b); Breugem (2012)

$\rho_{p}V_{p}\frac{{\rm d}\bm{u}_{p}}{{\rm d}t}=\rho_{f}\frac{{\rm d}}{{\rm
d}t}\left(\int_{V_{p}}{\bm{u}dV}\right)-\rho_{f}\sum_{l=1}^{N_{l}}\bm{F}_{p,l}\Delta
V_{l}+\left(\rho_{p}-\rho_{f}\right)V_{p}\bm{g}+{\bm{F}_{p,lub}+\bm{F}_{p,col}},$
(5) $\bm{I}_{p}\frac{{\rm d}\bm{\omega}_{p}}{{\rm d}t}=\rho_{f}\frac{{\rm
d}}{{\rm
d}t}\left(\int_{V_{p}}{\bm{r}\times\bm{u}dV}\right)-\rho_{f}\sum_{l=1}^{N_{l}}\left(\bm{r}_{l}\times\bm{F}_{p,l}\right)\Delta
V_{l}+{\bm{T}_{p,col}},$ (6)

where $N_{l}$ is the number of Lagrangian grid points, $F_{p,l}$ is the volume
force on the Lagrangian grid point calculated by equation (B.3) of the IB
method, and $\Delta V_{l}=\pi\Delta x\left(12R_{p}^{2}+\Delta
x^{2}\right)/(3N_{l})$ is the volume of a Lagrangian grid cell Uhlmann (2005).
The discretization schemes for these governing equations are given in Appendix
B with details. According to the approximation of the hydrodynamic force Kempe
and Fröhlich (2012a); Costa et al. (2015); Biegert et al. (2017), the
hydrodynamic force is excluded during contact for particles with the Stokes
number $St=\rho_{p}u_{in}D_{p}/(9\rho_{f}\nu_{f})\geqslant 5$, where $u_{in}$
is the particle impact velocity.

###  2.2 Lubrication force

When particles are close to colliding in particle-resolved simulations, the
evaluation of the forces acting on a particle can be performed as follows
Breugem (2010); Kempe and Fröhlich (2012a); Costa et al. (2015); Biegert et
al. (2017): (a) when $\varepsilon\geqslant\varepsilon_{\Delta x}$, the IB
method or a similar approach can be directly used to calculate the
hydrodynamic force on the particle, where $\varepsilon=\delta_{n}/R_{p}$ is
the non-dimensional gap width, $\delta_{n}$ is the distance between two
surfaces in a collision, $R_{p}$ is particle radius, and $\varepsilon_{\Delta
x}$ is a predefined constant. (b) When $0\leq\varepsilon<\varepsilon_{\Delta
x}$, the IB method or a similar approach may underpredict the hydrodynamic
force due to a lack of spatial grid resolution. As a remedy, a lubrication
model based on asymptotic expansions of analytical solutions for the
lubrication force in the Stokes regime can be employed Brenner (1961); Cox and
Brenner (1967); Cooley and O’Neill (1969). (c) When $\varepsilon<0$, the
collision model takes control of the particle, and the hydrodynamic force is
neglected. It should be noted that neglecting the hydrodynamic force during a
collision is a model assumption.

When a particle approaches another particle or a wall with a finite relative
velocity, fluid is squeezed out of the gap and the lubrication model can be
employed to calculate the hydrodynamic force with the Stokes flow assumption
Brenner (1961); Cox and Brenner (1967); Cooley and O’Neill (1969). The total
lubrication force $\bm{F}_{p,lub}$ acting on the particle $p$ can be computed
as follows:

${\bm{F}_{p,lub}=\sum_{p,q\neq
p}^{N_{p}}\bm{F}_{n,pq}^{lub}+\bm{F}_{n,pw}^{lub},}$ (7)

where $\bm{F}_{n}^{lub}$ is the normal lubrication force for a collision pair,
$N_{p}$ is the total particle number, the subscript $n$ indicates the normal
direction, and the subscripts $pq$ and $pw$ indicate the contact of particle
$p$ with particle $q$ and the wall, respectively. Here, $\bm{F}_{n}^{lub}$ is
calculated by the following the two-parameter lubrication model proposed by
Breugem (2010) and Costa et al. (2015):

$\displaystyle\bm{F}_{n}^{lub}=\begin{cases}-6\pi\mu_{f}R_{p}\bm{u}_{cp,n}\left(\lambda(\varepsilon)-\lambda(\varepsilon_{\Delta
x})\right),\quad\varepsilon_{\sigma}\leqslant\varepsilon<\varepsilon_{\Delta
x}\\\
-6\pi\mu_{f}R_{p}\bm{u}_{cp,n}\left(\lambda(\varepsilon_{\sigma})-\lambda(\varepsilon_{\Delta
x})\right),\quad 0\leqslant\varepsilon<\varepsilon_{\sigma}\\\
0,\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\ \
\text{otherwise}\end{cases},$ (8)
$\displaystyle\lambda_{pq}(\varepsilon)=\frac{1}{2\varepsilon}-\frac{9}{20}\ln{\varepsilon}-\frac{3}{56}\varepsilon\ln{\varepsilon}+O(1),$
(9)
$\displaystyle\lambda_{pw}(\varepsilon)=\frac{1}{\varepsilon}-\frac{1}{5}\ln{\varepsilon}-\frac{1}{21}\varepsilon\ln{\varepsilon}+O(1),$
(10)

where $\bm{u}_{cp,n}$ is the normal component of the relative velocity, as
defined in Appendix A; $\lambda_{pq}$ and $\lambda_{pw}$ are the Stokes
amplification factors for the lubrication force in the interactions between
particle $p$ and particle $q$ and between particle $p$ and the wall,
respectively; $\varepsilon_{\Delta x}$ is the non-dimensional gap spacing
related to the grid size (_i.e._ , $D_{p}/\Delta x$); $\varepsilon_{\sigma}$
is a fixed minimum gap spacing related to asperities. Following Costa et al.
(2015), we adopted $\varepsilon_{\Delta x,pw}=0.075$, $\varepsilon_{\Delta
x,pq}=0.025$ for $D_{p}/\Delta x=20$ and $\varepsilon_{\Delta x,pw}=0.05$,
$\varepsilon_{\Delta x,pq}=0.025$ for $D_{p}/\Delta x=30$. In addition,
$\varepsilon_{\sigma,pw}=0.0008$ and $\varepsilon_{\sigma,pp}=0.0001$ were
calibrated by matching the trajectories of the particle–wall and particle-
particle collisions with the experimental data of Gondret et al. (2002) in the
present study. These trajectories are shown in Section 4.1.

###  2.3 Collision force

The total collision force $\bm{F}_{p,col}$ and total torque $\bm{T}_{p,col}$
acting on particle $p$ are computed as follows:

${\bm{F}_{p,col}=\sum_{p,q\neq
p}^{N_{p}}\left(\bm{F}_{n,pq}^{col}+\bm{F}_{t,pq}^{col}\right)+\bm{F}_{n,pw}^{col}+\bm{F}_{t,pw}^{col},}$
(11) ${\bm{T}_{p,col}=\sum_{p,q\neq
p}^{N_{p}}R_{p}\bm{n}_{pq}\times\bm{F}_{t,pq}^{col}+R_{p}\bm{n}_{pw}\times\bm{F}_{t,pw}^{col},}$
(12)

where $\bm{F}^{col}$ is the collision force for a collision pair, $\bm{n}$ is
the normal unit vector of the contact, and the subscripts $n$ and $t$ indicate
the normal and tangential directions, respectively. The collision force
$\bm{F}^{col}$ is calculated by the soft-sphere model. The model is composed
of a linear spring–dashpot system in the normal direction, and a linear
spring–dashpot system with a Coulomb friction slider in the tangential
direction Costa et al. (2015). The normal collision force $\bm{F}_{n}^{col}$
and the tangential collision force $\bm{F}_{t}^{col}$ are determined by

$\displaystyle\bm{F}_{n}^{col}=-k_{n}\left|\delta_{n}\right|\bm{n}-d_{n}\bm{u}_{cp,n},$
(13) $\displaystyle\bm{F}_{t}^{col}=-k_{t}\bm{\delta}_{t}-d_{t}\bm{u}_{cp,t},$
(14)

where $k_{n}$ and $k_{t}$ are the stiffnesses in the normal and tangential
directions, respectively; $d_{n}$ and $d_{t}$ are the damping coefficients in
the normal and tangential directions, respectively; $\delta_{n}$ is the
distance between two surfaces; $\bm{\delta}_{t}$ and $\bm{u}_{cp,t}$ are the
tangential displacement of the spring and the tangential component of the
relative velocity, respectively, as defined in Appendix A.

The stiffness and damping coefficients can be derived by the collision time
and the restitution coefficient Costa et al. (2015), as

$k_{n}=\frac{m_{e,n}\left(\pi^{2}+\ln^{2}e_{n,d}\right)}{T_{c}},\quad
d_{n}=-\frac{2m_{e,n}\ln e_{n,d}}{T_{c}},$ (15)
$k_{t}=\frac{m_{e,t}\left(\pi^{2}+\ln^{2}e_{t,d}\right)}{T_{c}},\quad
d_{t}=-\frac{2m_{e,t}\ln e_{t,d}}{T_{c}},$ (16)

where $m_{e,n}=(m_{p}^{-1}+m_{q}^{-1})^{-1}$ and
$m_{e,t}=(1+1/K^{2})^{-1}m_{e,n}$ are the reduced masses of the linear
spring–dashpot system; $m_{p}$ and $m_{q}$ are the masses of particles $p$ and
$q$, respectively; $K$ is the normalized particle radius of gyration with
$K^{2}=2/5$ for a homogeneous spherical particle; $e_{n,d}$ and $e_{t,d}$ are
the dry coefficients of restitution in the normal and tangential directions,
respectively.

The stretching collision time model determines the stiffness and damping
coefficients through a nonphysical stretched collision time $T_{c}$,
significantly reducing the stiffness and damping coefficients. To avoid
introducing the stretched collision time, the physical collision time should
be used Li et al. (2020). The physical collision time $T_{c}$ can be
calculated by the physical properties of particle as Zenit et al. (1997)

$T_{c}=7.894\left(\frac{m_{p}^{2}}{E^{*2}u_{in}D_{p}}\right)^{1/5},$ (17)

where $E^{*}$ is the equivalent elasticity of the collision system depending
on Young’s modulus $E$ and Poisson’s ratio $\nu$, which is expressed as
follows:

$\frac{1}{E^{*}}=\frac{1}{\pi}\left(\frac{1+\nu_{p}^{2}}{E_{p}}+\frac{1+\nu_{q}^{2}}{E_{q}}\right).$
(18)

In the present study, we use $E=200$ GPa and $\nu=0.3$ for steel particles and
$E=55$ GPa and $\nu=0.25$ for glass particles. The ratio between $T_{c}$ and
the particle response time $\tau_{p}=\rho_{p}D_{p}^{2}/18\mu_{f}$ is
$O(10^{-4}\sim 10^{-3})$ for $10$-mm steel/glass particles with $St=100$
impacting in air or water, respectively.

When
$\left|\bm{F}_{t}^{col}\right|\geqslant\mu_{c}\left|\bm{F}_{n}^{col}\right|$
($\mu_{c}$ is the friction coefficient), a particle starts to slide. The
tangential force is controlled by the slider in the tangential direction based
on the Coulomb friction law instead of the tangential displacement, as

$\bm{F}_{t}^{col}=\mu_{c}|\bm{F}_{n}^{col}|\bm{t}.$ (19)

where $\bm{t}$ is the unit vector in the tangential direction. The tangential
displacement needs to be reset in order to comply with Coulomb’s friction law,
which is

$\bm{\delta}_{t}=-\frac{\mu_{c}|\bm{F}_{n}^{col}|\bm{t}+d_{t}\bm{u}_{cp,t}}{k_{t}}.$
(20)

##  3 The multiple-time-step integration algorithm (MTSA)

###  3.1 The motivation

The mismatch problem of the characteristic time scales between particle
collision and fluid flow is crucial for particle-resolved simulations. In
fact, the characteristic time scales for fluid flow, fluid-particle
interaction and particle collision are quite different in magnitude. We have
estimated that the particle collision time can be 3 or 4 orders of magnitude
smaller than the particle response time, and thus 1 or 2 orders of magnitude
smaller than the flow time scale with $St\sim O(100)$. Using a very fine fluid
flow time step $\Delta t_{f}$ to resolve particle collision can lead to a
highly expensive computation. When employing a single large time step in a
particle-resolved simulation, the error of the soft-sphere model is also large
and mainly comes from two aspects. First, the particle motion cannot be well
resolved by a large flow time step during a collision. It will result in an
excessive overlap between colliding surfaces and overpredicting the rebounding
trajectory, as shown in figure 1 of Kempe and Fröhlich (2012a). Second,
particle velocity rapidly changes during a collision. The surrounding fluid
can hardly follow this rapid change with a large flow time step. It may cause
an overestimation of the hydrodynamic force by the IB method and underpredicts
the rebounding trajectory, as shown in figure 8 of Costa et al. (2015). These
two facts prevent the usage of the traditional soft-sphere model with a single
large time step. On the other hand, as we have discussed in the introduction,
artificially reducing particle stiffness like the SCTA is nonphysical and may
also bring in errors in simulations Xia et al. (2020).

Due to the aforementioned drawbacks of the existing methods, we propose the
MTSA employing a normal fluid flow time step while inserting additional two-
level time substeps to resolve fluid-particle interaction and particle motion
during a collision, so that the prediction errors by increasing particle
stiffness can be significantly reduced and the computation efficiency is not
degraded too much at the same time.

Figure 1: Sketch of the multiple-time-step integration algorithm (MTSA).

###  3.2 The computation procedure

The MTSA was activated once particle-particle/wall collision happens. A single
time step $\Delta t$ is replaced by three different time steps, namely a fluid
flow time step $\Delta t_{f}$ and additional two-level time substeps (_i.e._ ,
fluid-particle interaction time substep $\Delta t_{i}$ and particle motion
time substep $\Delta t_{p}$) hierarchically inserted into it. As shown in
figure 1, $\Delta t_{f}=T_{c}/R_{f}$ is employed to advance fluid flow,
$\Delta t_{i}=T_{c}/R_{i}$ to advance fluid-particle interaction, and $\Delta
t_{p}=T_{c}/R_{m}$ to advance particle motion. Here, $R_{f}$, $R_{i}$ and
$R_{m}$ are the iteration numbers of fluid flow, fluid-particle interaction
and particle motion substeps in one collision time $T_{c}$, respectively.
Based on the above three parameters, $N_{i}=R_{i}/R_{f}$ and
$N_{m}=R_{m}/R_{f}$ can be defined as the iteration numbers of fluid-particle
interaction and particle motion substeps in one fluid flow time step. It
should be noted that $N_{i}$ and $N_{m}$ should be integers, which restrict
the choices for $R_{i}$ and $R_{m}$.

The pressure Poisson equation is only solved at the fluid flow step $t=t_{f}$,
which greatly improves the computational efficiency of the MTSA. Local
pressure is assumed to be a constant at the fluid-particle interaction and
particle motion substeps in a fluid flow step. At a fluid-particle interaction
substep $t=t_{i}$, we solve the flow momentum equation, calculate the IB
force, and modify the fluid velocity by the IB force accordingly. At a
particle motion substep $t=t_{p}$, we calculate particle collision and motion.
When the time is advanced to the same instant, the information of the fluid
flow and particle will be updated and synchronized. The updating and
synchronizing process are demonstrated by the arrow lines in figure 1.
Particularly, when $t=t_{f}=t_{i}=t_{p}$, the particle position and velocity
at $t=t_{p}$ are updated to the particle-fluid interaction substep before
calculating the IB force. Second, the IB force is calculated and used to
modify the fluid velocity at $t=t_{i}$. Third, the calculated IB force at
$t=t_{i}$ is synchronized to the fluid flow step as the input for particle
motion, and the modified fluid velocity at $t=t_{i}$ is also updated as the
input for the Poisson equation. Fourth, the Poisson equation and flow field at
$t=t_{f}$ are solved and updated. Finally, the updated fluid velocity at
$t=t_{f}$ is synchronized to $t=t_{i}$ and $t=t_{p}$ as the input for the flow
momentum equation and particle motion. When $t\neq t_{f}$, we only update the
particle position and velocity.

In the MTSA, an explicit second-order Runge–Kutta (RK2) method is used for the
temporal advancement of the fluid flow Yang et al. (2017, 2018); Cui et al.
(2018); He et al. (2022); Zhu et al. (2022). At each substep of the RK2
method, the fractional-step method of Kim and Moin (1985) is applied to ensure
that the flow field is divergence-free. In order to insert enough substeps to
ensure that the flow can adapt to the change of particle motion during
collisions and reduce the computational cost, the explicit Euler method is
adopted in the particle-fluid interaction and particle-motion substeps for
time advancement. The discretized equations are presented in Appendix B.

###  3.3 The advantages

In summary, the proposed MTSA has the following advantages:

(1) Compared with the traditional soft-sphere model, the MTSA significantly
improves computational efficiency as it employs a large fluid flow time step
$\Delta t_{f}$. The pressure Poisson equation is only solved at the fluid flow
time step $t=t_{f}$, which greatly saves the computation cost. The
computational time of the MTSA can be reduced to about one order of magnitude
less than that of the traditional soft-sphere model with almost the same
accuracy, which will be shown in §4.

(2) Compared with the SCTA, the MTSA employs physical particle stiffness
without reducing particle stiffness. It is expected that the problem of
hydrodynamic force underestimation will be alleviated or resolved by the MTSA,
and the simulation of sediment transport will not be affected by nonphysically
reduced particle stiffness. In addition, the parameters $R_{i}$ in the MTSA
has very good universality for different fluid time steps, which will also be
shown in §4.

In addition to the aforementioned algorithms proposed in particle-resolved
simulations, some numerical algorithms have also been proposed in point
particle simulations to address the mismatch of the time scales between
particle collision and the fluid flow Deen et al. (2007); Capecelatro and
Desjardins (2013); Finn et al. (2016). They generally employed two-time
scales, where a large time step was adopted to advance fluid flow and the
inter-phase coupling, while a small time step was employed for particle
motion. Some researchers even used a smaller time scale to calculate particle
collisions Darmana et al. (2006); Marshall (2009). However, as will be shown
in §4, these algorithms are ineffective in particle-resolved simulations,
since fluid-particle interaction is not resolved during a collision.

##  4 Results

In this section, the accuracy and convergence of the traditional soft-sphere
model and the MTSA are assessed through simulations of particle-wall and
particle-particle collisions. Furthermore, we discuss the impact of $R_{f}$
and $R_{i}$ on the prediction results using the MTSA. Finally, a sediment
transport case with thousands of finite-size particles is simulated as a
practical application of the MTSA and a demonstration of the difference
between the MTSA and the SCTA. It should be noted that a version of SCTA
similar to Costa et al. (2015) is adopted here to be compared with the MTSA,
which incorporates a two-parameter lubrication model, a linear spring-dashpot
system, and substeps for particle collisions and motions.

Because the hydrodynamic force is neglected in a particle collision, which is
totally governed by the spring-dashpot system, $R_{m}$ is not influenced by
the grid size and $St$. We directly determine $R_{m}$ by normal dry particle-
wall and particle-particle collision cases. The simulation results of the
restitution coefficient $e_{n}$ were compared with a given $e_{n,d}=0.97$ for
steel and glass particles. For both $R_{m}=40$ and $80$, it is found that the
relative error between $e_{n}$ and $e_{n,d}$ is less than $0.1\%$, which is
sufficiently small. Therefore, the impact of $R_{m}$ is not studied in this
section. And $R_{m}=40$ is used in the following cases, and $R_{m}=80$ is only
used when $R_{i}=16$.

###  4.1 Accuracy and convergence of the traditional soft-sphere model

Cases REF, S1 and S2 in table 1 are designed to verify the accuracy and
convergence of the traditional soft-sphere model as the benchmark result.
These three cases can be regarded as using the MTSA without the fluid-particle
interaction substep (_i.e._ , $N_{i}=R_{i}/R_{f}=1$), which is similar to the
algorithm proposed for the point-particle simulations Deen et al. (2007).
Therefore, there are only two-time steps in cases REF, S1 and S2.
Particularly, the Navier-Stokes and fluid-particle interaction equations are
solved with the flow time step $\Delta t_{f}$, while the Newton-Euler
equations are advanced with the particle motion time step $\Delta t_{p}$. Case
REF is regarded as the benchmark case. In case S1, $\Delta t_{f}$ was smaller
than that in case REF, which is used to verify the accuracy and convergence of
case REF. In case S2, $\Delta t_{f}$ is larger than that in case REF, which is
used to compare with the results of the MTSA with the same $\Delta t_{f}$ to
demonstrate the improvement by inserting additional fluid-particle interaction
substeps.

Table 1: The MTSA parameters $R_{f}$, $R_{i}$, and $R_{m}$ used in different cases of particle-wall and particle-particle collisions. $N_{i}=R_{i}/R_{f}=1$ indicates that the traditional soft-sphere model is used without the fluid-particle interaction (FPI) substep. Case | $R_{f}$ | $R_{i}$ | $R_{m}$ | Notes
---|---|---|---|---
REF | 8 | 8 | 40 | Benchmark case
S1 | 16 | 16 | 80 | Reduce $\Delta t_{f}$, without FPI substep
S2 | 2 | 2 | 40 | Increase $\Delta t_{f}$, without FPI substep
MTSA | 0.5,1,2 | 2,4,8,16 | 40,80 | With FPI substep

First, the bouncing motion of a single particle in a viscous fluid with
different $St$ and $D_{p}/\Delta x$ values were simulated to assess the
accuracy and convergence of the traditional soft-sphere model in particle-wall
collisions. The physical parameters in the simulation were comparable to the
experiment of Gondret et al. (2002): $\rho_{p}=7800$ kg/m3 and $e_{n,d}=0.97$,
and other parameters are listed in table 2. The computational domain was
$12.8D_{p}\times 25.6D_{p}\times 12.8D_{p}$. Two grid resolutions,
$D_{p}/\Delta x=20$ and $D_{p}/\Delta x=30$, were implemented in this
simulation. A periodic boundary condition was imposed in the horizontal
directions, and a no-slip boundary condition was imposed on both the top and
bottom surfaces. A particle was initially placed at $x=L_{x}/2$,
$y=L_{y}-0.75D_{p}$, and $z=L_{z}/2$. The falling velocity was prescribed
following Biegert et al. (2017), where it accelerated smoothly and $u_{in}$
matched the Stokes number in the previously reported experiment Gondret et al.
(2002) before the collision, as

$u_{p}(t)=u_{in}\left(e^{-40t}-1\right),\quad\text{if}\quad\delta_{n}>R_{p}.$
(21)

Once the particle reached a distance of $\delta_{n}=R_{p}$, we turned off the
prescribed velocity. Then, the particle moved under hydrodynamic,
gravitational, buoyant, and collision forces. The normal coefficient of
restitution is defined as $e_{n}=u_{out}^{*}/u_{in}^{*}$. Following Costa et
al. (2015), $u_{in}^{*}$ and $u_{out}^{*}$ are defined at the instants
$t-t_{c}=\mp f_{s}^{-1}$, where $t_{c}$ is the instant of collision and
$f_{s}$ is the sampling frequency in measurements, which is 500 Hz in the
experiment of particle-wall collision Joseph et al. (2001) and 100 Hz in the
experiment of particle-particle collision Yang and Hunt (2006).

Table 2: Parameters used in the simulation of particle-wall collisions. $St$ |  $D_{p}$ (mm) |  $\rho_{f}$ (kg/m3) |  $\mu_{f}$ (cP) |  $u_{in}$ (m/s) |  $T_{c}$ (s)
---|---|---|---|---|---
6 | 3 | 965 | 100 | 0.231 | $2.3\times 10^{-5}$
15 | 6 | 965 | 100 | 0.288 | $4.4\times 10^{-5}$
27 | 6 | 965 | 100 | 0.519 | $3.9\times 10^{-5}$
60 | 3 | 953 | 20 | 0.462 | $2.0\times 10^{-5}$
100 | 4 | 953 | 20 | 0.577 | $2.6\times 10^{-5}$
152 | 3 | 935 | 10 | 0.585 | $1.9\times 10^{-5}$
193 | 6 | 953 | 20 | 0.742 | $3.6\times 10^{-5}$
742 | 5 | 920 | 5 | 0.856 | $3.0\times 10^{-5}$

The resulting normal coefficients of restitution $e_{n}$ are compared with the
experimental data of Joseph et al. (2001) and Gondret et al. (2002) here. The
simulation results of the coefficient of restitution $e_{n}$ for normal
particle-wall collisions by the traditional soft-sphere model are shown in
figure 2. The resulting $e_{n}$ of case REF agrees well with the experimental
results over the entire range of the particle Stokes number. This verifies the
accuracy of the traditional soft-sphere model. The resulting $e_{n}$ values
are almost the same when refining the grid (_i.e._ , $D_{p}/\Delta x=30$) at
$St=15,27,152,742$ or refining the flow time step $\Delta t_{f}$ (_i.e._ ,
case S1) at $St=27,152$. This verifies the convergence of the traditional
soft-sphere model. However, $e_{n}$ is evidently underestimated with a much
larger $\Delta t_{f}$ (case S2). This is because $\Delta t_{f}$ is too large
to resolve the fluid-particle interaction so that the surrounding fluid cannot
follow the rapid change of the particle velocity during the collision, and the
velocity difference between the particle and fluid after the collision is
increased. This leads to an overestimation of the hydrodynamic force by the IB
method. The underestimated $e_{n}$ in case S2 can cause a remarkable
discrepancy in the particle bouncing trajectories, as shown in figure 3. It is
seen that the trajectory of case S1 is almost the same as that of case REF,
because the values of $e_{n}$ are close in the two cases. The trajectory of
the refined grid (_i.e._ , $D_{p}/\Delta x=30$) is not shown here, as the
corresponding $e_{n}$ and trajectory are very close to those of cases REF and
S1.

Figure 2: Wet coefficient of restitution for normal particle-wall collisions
by the traditional soft-sphere model.

Figure 3: Comparison of the particle trajectories between the experiment and
the simulation results by the traditional soft-sphere model with different
Stokes numbers and time steps: (a) $St=27$ and (b) $St=152$.
$t_{c}=t\mid_{\delta_{n}=0}$ is the instant of collision and
$t_{ref}=\sqrt{D_{p}/\left|\bm{g}\right|}$ is the reference time scale.

Next, a moving particle colliding with a steady particle in a viscous fluid
with different $St$ and $D_{p}/\Delta x$ values was simulated to assess the
accuracy and convergence of the traditional soft-sphere model in particle-
particle collisions. The physical parameters in the simulations were
comparable to the experiments of Yang and Hunt (2006): $D_{p}=12.7$ mm,
$\rho_{p}=7780$ kg/m3, $\rho_{f}=1125$ kg/m3, $\mu_{f}=45$ cP and
$e_{n,d}=0.97$, and other parameters are listed in table 3. The computational
domain was $6D_{p}\times 24D_{p}\times 6D_{p}$. Two grid resolutions of
$D_{p}/\Delta x=20$ and $30$ were implemented in the simulations. The periodic
boundary conditions were applied in all three directions. The steady and
moving particles were initially placed at $y=19D_{p}$ and $2D_{p}$,
respectively, centered at $x=L_{x}/2$ and $z=L_{z}/2$. The velocity of the
moving particle was prescribed by the negative form of equation (21) since it
was moving upward. The gravitational and buoyancy forces were not considered.

Table 3: Parameters used in the simulation of particle-particle collisions. $St$ |  $u_{in}$ (m/s-1) |  $T_{c}$ (s)
---|---|---
12.7 |  0.059 |  $1.3\times 10^{-4}$
21.5 |  0.099 |  $1.2\times 10^{-4}$
34.3 |  0.158 |  $1.1\times 10^{-4}$
52.7 |  0.243 |  $9.6\times 10^{-5}$
135.2 |  0.623 |  $8.0\times 10^{-5}$
345 |  1.591 |  $6.6\times 10^{-5}$

In figure 4 (a), the resulting normal coefficients of restitution $e_{n}$ in
the particle-particle collisions are compared with the experimental data of
Yang and Hunt (2006) and the numerical results of Costa et al. (2015). The
conclusions for the particle-particle collisions are similar to those for
particle-wall collisions. The resulting $e_{n}$ of case REF agrees well with
the experimental and numerical results over the entire range of the Stokes
number $St$. The resulting $e_{n}$ is almost unchanged by refining the grid
($D_{p}/\Delta x=30$) at $St=12,34,135,345$ or reducing the flow time step
$\Delta t_{f}$ (case S1) at $St=34$, which verifies the accuracy and
convergence of the traditional soft-sphere model. Although the resulting
$e_{n}$ by increasing $\Delta t_{f}$ (case S2) is within the range of the
widely scattered experimental results, there still was a visible deviation
compared with the results of case REF due to the insufficient temporal
resolution for fluid-particle interactions. Figure 4 (b) further shows the
contact point trajectories of particles at $St=34$. The trajectories in cases
REF and S1 have good agreement with those of Costa et al. (2015). It is seen
that the underestimated $e_{n}$ of case S2 causes a visible discrepancy in the
trajectories of the contact point compared with case REF.

The results of the particle-wall and particle-particle collisions indicate
that the traditional soft sphere model can be used to reproduce experimental
results only if the flow time step $\Delta t_{f}$ is sufficiently small, and
increasing $\Delta t_{f}$ would result in an underestimation of $e_{n}$ and
non-negligible errors in the particle trajectories. However, the computational
cost of using such a small $\Delta t_{f}$ is quite high in particle-resolved
simulations with a large number of particles. Hence, the MTSA can be a good
method to reduce the computational cost of the traditional soft-sphere model
while maintaining its accuracy. In the following, the results of case REF will
be employed as a benchmark to test the accuracy of the MTSA.

Figure 4: Simulation results of the particle-particle collisions by the
traditional soft-sphere model: (a) wet coefficient of restitution in normal
particle-particle collisions and (b) trajectories of particle contact points
at $St=34$. $t_{c}=t\mid_{\delta_{n}=0}$ is the instant of collision, and
$t_{ref}=\sqrt{D_{p}/\left|\bm{g}\right|}$ is the reference time scale.

###  4.2 Assessment of the MTSA

Based on the above results that a large flow time step $\Delta t_{f}$ cannot
be used to resolve fluid-particle interaction during a collision, we employ a
large $\Delta t_{f}$ for fluid flow advancement while inserting additional
time substeps for fluid-particle interaction in the MTSA. In this subsection,
we investigate the influence of the substep iteration number on the simulation
results with different particle Stokes number $St$ and grid resolution
$D_{p}/\Delta x$.

####  4.2.1 Particle-wall collision

For particle-wall collision, we use the same simulation setup and parameters
as in §4.1. The resulting normal coefficient of restitution $e_{n}$ using the
MTSA is compared with that of case REF. Figure 5 shows the error of $e_{n}$
between the MTSA and case REF with different $St$ and $D_{p}/\Delta x$ values.
The left and right columns show the errors with grid resolutions of
$D_{p}/\Delta x=20$ and $D_{p}/\Delta x=30$, respectively. The dashed line
corresponds to zero error.

It can be seen that the general trends of the error curves with $R_{i}$ in
each subfigure are similar in figures 5 (a)-(f). As discussed before, $e_{n}$
is underestimated ($e_{n,MTSA}-e_{n,REF}<0$) without including additional
fluid-particle interaction substeps, since the fluid flow surrounding the
particle cannot adapt to the rapid change of the particle velocity. Therefore,
the error curves are mostly under the zero-error line at small $R_{i}$ values.
If more fluid-particle interaction substeps are inserted (increasing $R_{i}$),
the absolute value of $e_{n,MTSA}-e_{n,REF}$ decreases first. However, as the
pressure Poisson equation is not solved in the particle-fluid interaction
substeps, the intermediate flow velocities are nonphysical, which leads to an
overestimation of $e_{n}$ at larger $R_{i}$ values. Therefore, the overall
trend of the absolute error of $e_{n,MTSA}-e_{n,REF}$ decreases first then
increases with $R_{i}$, as shown in figures 5 (a)-(f).

$D_{p}/\Delta x=20$

$D_{p}/\Delta x=30$

Figure 5: Prediction error of the normal coefficient of restitution $e_{n}$
between the MTSA and case REF in particle–wall collisions: (a,b) $R_{f}=0.5$,
(c,d) $R_{f}=1.0$, (e,f) $R_{f}=2.0$, and (g,h) the contours of the maximum
error $\left|e_{n,MTSA}-e_{n,REF}\right|$ with $St=15,27,152,742$.

(a,c,d,g) $D_{p}/\Delta x=20$ and (b,d,f,h) $D_{p}/\Delta x=30$. The dashed
line corresponds to zero error.

We can elucidate the increase in the error of $e_{n,MTSA}-e_{n,REF}$ with
$R_{i}$ as follows. Based on the Navier-Stokes equations discretized by the
explicit Euler method, the intermediate flow velocities at the $i$th fluid-
particle interaction substep during a collision can be expressed as

$\displaystyle\bm{u}^{i}$ $\displaystyle=\bm{u}^{i-1}+\Delta
t_{i}\left(\bm{H}^{i-1}+\bm{f}^{i-1}-\frac{1}{\rho_{f}}\nabla p^{i-1}\right),$
(22) $\displaystyle\widetilde{\bm{u}}^{i}$
$\displaystyle=\widetilde{\bm{u}}^{i-1}+\Delta
t_{i}\left(\widetilde{\bm{H}}^{i-1}+\widetilde{\bm{f}}^{i-1}\right),$ (23)

where $\bm{u}$ and $\widetilde{\bm{u}}$ are the physical (by solving the
pressure Poisson equation) and nonphysical (without solving the pressure
Poisson equation) intermediate flow velocities, respectively; $\bm{H}$ and
$\widetilde{\bm{H}}$ are the sums of the convective and viscous terms with
physical and nonphysical intermediate flow velocities, respectively; $\bm{f}$
and $\widetilde{\bm{f}}$ are the volume forces with physical and nonphysical
intermediate flow velocities, respectively. Assuming that the initial flow
velocity before a collision is $\bm{u}^{0}$, the flow velocity after the
collision ($i=R_{i}$) can be expressed as

$\displaystyle\bm{u}^{R_{i}}$
$\displaystyle=\bm{u}^{0}+\sum_{i=1}^{R_{i}}\Delta
t_{i}\left(\bm{H}^{i-1}+\bm{f}^{i-1}-\frac{1}{\rho_{f}}\nabla p^{i-1}\right),$
(24) $\displaystyle\widetilde{\bm{u}}^{R_{i}}$
$\displaystyle=\bm{u}^{0}+\sum_{i=1}^{R_{i}}\Delta
t_{i}\left(\widetilde{\bm{H}}^{i-1}+\widetilde{\bm{f}}^{i-1}\right).$ (25)

The quantities at Eulerian grid points can be transferred to Lagrangian grid
points through a regularized Dirac delta function $\delta_{d}$ Roma et al.
(1999). The flow velocity on the particle surface after the collision is

$\displaystyle\bm{U}_{l}^{R_{i}}$
$\displaystyle=\bm{U}_{l}^{0}+\sum_{i=1}^{R_{i}}\Delta
t_{i}\left(\bm{H}_{l}^{i-1}+\bm{F}_{p,l}^{i-1}-\frac{1}{\rho_{f}}\nabla
P_{l}^{i-1}\right),$ (26) $\displaystyle\widetilde{\bm{U}}_{l}^{R_{i}}$
$\displaystyle=\bm{U}_{l}^{0}+\sum_{i=1}^{R_{i}}\Delta
t_{i}\left(\widetilde{\bm{H}}_{l}^{i-1}+\widetilde{\bm{F}}_{p,l}^{i-1}\right),$
(27)

where $\bm{U}_{l}$, $\bm{F}_{p,l}$, $\bm{H}_{l}$, and $P_{l}$ are the
quantities of $\bm{u}$, $\bm{f}$, $\bm{H}$, and $p$ at Eulerian grid points
interpolated at Lagrangian grid points, respectively; and
$\widetilde{\bm{U}}_{l}$, $\widetilde{\bm{F}}_{p,l}$, and
$\widetilde{\bm{H}}_{l}$ are the quantities of $\widetilde{\bm{u}}$,
$\widetilde{\bm{f}}$, and $\widetilde{\bm{H}}$ on Eulerian grid points
interpolated on Lagrangian grid points, respectively. During a collision, the
particle motion is only determined by the collision model and not affected by
the fluid flow Kempe and Fröhlich (2012a); Costa et al. (2015); Biegert et al.
(2017), thus the particle velocities with physical and nonphysical
intermediate flow velocities are the same
($\bm{U}_{p}^{R_{i}}=\widetilde{\bm{U}}_{p}^{R_{i}}$). Based on the IB method
as in equation (B.3), the differences in the hydrodynamic force with physical
and nonphysical intermediate flow velocities after a collision can be
expressed as

$\displaystyle\Delta\bm{F}_{p,l}^{R_{i}}$
$\displaystyle=\widetilde{\bm{F}}_{p,l}^{R_{i}}-\bm{F}_{p,l}^{R_{i}}$
$\displaystyle=\frac{\widetilde{\bm{U}}_{p}^{R_{i}}-\bm{U}_{p}^{R_{i}}+\bm{U}_{l}^{R_{i}}-\widetilde{\bm{U}}_{l}^{R_{i}}}{\Delta
t_{i}}$ (28)
$\displaystyle=\frac{\bm{U}_{l}^{R_{i}}-\widetilde{\bm{U}}_{l}^{R_{i}}}{\Delta
t_{i}}.$

Substituting equations (26) and (27) into equation (28) yields

$\Delta\bm{F}_{p,l}^{R_{i}}=\sum_{i=1}^{R_{i}}(\Delta\bm{H}_{l}^{i-1}-\Delta\bm{F}_{p,l}^{i-1}-\frac{1}{\rho_{f}}\nabla
P_{l}^{i-1}).$ (29)

Equation (29) includes three terms:
$\Delta\bm{H}_{l}^{i-1}=\widetilde{\bm{H}}_{l}^{i-1}-\bm{H}_{l}^{i-1}$,
$\Delta\bm{F}_{p,l}^{i-1}=\widetilde{\bm{F}}_{p,l}^{i-1}-\bm{F}_{p,l}^{i-1}$,
and $\nabla P_{l}^{i-1}$, representing the flow velocity error, the
hydrodynamic force error, and the pressure error at each fluid-particle
interaction substep, respectively. If
$\Delta\bm{H}_{l}^{i-1}-\Delta\bm{F}_{p,l}^{i-1}-\nabla
P_{l}^{i-1}/\rho_{f}\neq 0$, the error of $\Delta\bm{F}_{p,l}^{R_{i}}$
increases as $R_{i}$ increases, which leads to the error of $e_{n}$ increasing
with $R_{i}$. As shown in figures 5 (a-f) and figures 8 (a-f), $e_{n,MTSA}$
increases with $R_{i}$, which indicates that the MTSA underestimates the
hydrodynamic force as $R_{i}$ increased and
$\Delta\bm{H}_{l}^{i-1}-\Delta\bm{F}_{p,l}^{i-1}-\nabla
P_{l}^{i-1}/\rho_{f}<0$.

Since insufficient fluid-particle interaction substeps can cause an
underestimation of $e_{n}$, and too many substeps can cause an overestimation
of $e_{n}$. Therefore, we may find an optimal value of $R_{i}$ to minimize the
error of $e_{n}$. Figures 5 (g) and (h) show the contours of
$\max\left|e_{n,MTSA}-e_{n,REF}\right|$ for different $St$ and $D_{p}/\Delta
x$ values. It can be found that the maximum error is small for
$R_{i}\in\left\\{4,5,6\right\\}$ with different $St$ and $D_{p}/\Delta x$
values, which indicates that the MTSA has good accuracy with
$R_{i}\in\left\\{4,5,6\right\\}$ in the particle-wall collisions. Based on the
maximum error in the particle-particle collisions shown in figures 8 (g) and
(h), $R_{i}=4$ is determined to be the optimal value. The maximum error for
$R_{i}=4$ is less than 0.031 in the particle-wall collisions, which is small
enough.

After determining that $R_{i}=4$ is optimal, we employed the MTSA with
$R_{f}=1,R_{i}=4,R_{m}=40$ to calculate $e_{n}$ for particle-wall collisions
in the entire range of $St$. As shown in figure 6, the simulation results of
the MTSA agree well with those of case REF and the experimental results over
the entire range of $St$. Therefore, after inserting proper additional fluid-
particle interaction substeps, the prediction results by the MTSA are
significantly improved compared with those of case S2 in particle-wall
collisions. The yellow error bars in figure 6 indicate the prediction
uncertainty of $e_{n}$ with $R_{f}=0.5,1,2$. The short lengths of the error
bars imply that the MTSA has good convergence with these $R_{f}$ values. A
consistent conclusion can be further drawn from the particle trajectories
shown in figure 7. It can be seen that the particle trajectories predicted by
the MTSA show good agreements with those of case REF, and they are
significantly improved compared with those of case S2. The grey area is the
particle trajectory range with $R_{f}=0.5,1,2$. Such small areas demonstrate
that the particle trajectories in the particle-wall collisions will not be
much affected by $\Delta t_{f}$ if it is sufficiently small ($R_{f}\geqslant
0.5$).

Figure 6: Wet coefficient of restitution in normal particle-wall collisions
predicted by the MTSA ($R_{f}=1$, $R_{i}=4$, $R_{m}=40$). The yellow error bar
indicates the prediction uncertainty of $e_{n}$ with $R_{f}=0.5,1,2$.

Figure 7: Particle trajectories with different Stokes numbers: (a) $St=27$ and
(b) $St=152$. The grey area is the particle trajectory range with
$R_{f}=0.5,1,2$. $t_{c}=t\mid_{\delta_{n}=0}$ is the instant of collision, and
$t_{ref}=\sqrt{D_{p}/\left|\bm{g}\right|}$ is the reference time scale.

####  4.2.2 Particle-particle collision

For particle-particle collisions, we adopted the same simulation setup as that
in §4.1. Figure 8 shows the error of the normal coefficient of restitution
$e_{n}$ between the MTSA and case REF. The error curves in each subfigure
exhibit the same trend as those in particle-wall collisions, for the same
reason as that in particle-wall collisions discussed above.

$D_{p}/\Delta x=20$

$D_{p}/\Delta x=30$

Figure 8: Prediction error of the normal coefficient of restitution $e_{n}$
between the MTSA and case REF in particle-particle collisions: (a,b)
$R_{f}=0.5$, (c,d) $R_{f}=1.0$, (e,f) $R_{f}=2.0$, and (g,h) the contours of
the maximum error $\left|e_{n,MTSA}-e_{n,REF}\right|$ with $St=12,34,135,345$.

(a,c,d,g) $D_{p}/\Delta x=20$ and (b,d,f,h) $D_{p}/\Delta x=30$. The dashed
line corresponds to zero error.

Figures 8 (g) and (h) show the contours of
$\max\left|e_{n,MTSA}-e_{n,REF}\right|$ with different $St$ and $D_{p}/\Delta
x$ values. It is seen that the maximum error is small for
$R_{i}\in\left\\{4\right\\}$ with different $St$ and $D_{p}/\Delta x$ values,
which indicates that the MTSA has good accuracy with $R_{i}=4$ in particle-
particle collisions. Based on the combined results of the particle-wall and
particle-particle collisions, $R_{i}=4$ is determined to be the optimal value
and will be used in the following. The absolute maximum error for $R_{i}=4$ is
less than 0.024 in the particle-particle collisions, which is also small
enough.

We further employed the MTSA with $R_{f}=1,R_{i}=4,R_{m}=40$ to calculate
$e_{n}$ for particle-particle collisions over the entire range of $St$. As
shown in figure 9 (a), the simulation results of the MTSA agree well with
those of case REF and the experimental results. Its predictions are much
better than those of case S2. The yellow error bars represent the prediction
uncertainty of $e_{n}$ with $R_{f}=0.5,1,2$. The short lengths of the error
bars demonstrate that the MTSA has good convergence with these $R_{f}$ values.
Figure 9 (b) shows the particle trajectories with $St=34$. The particle
trajectories predicted by the MTSA have excellent agreement with those of case
REF and are much better than those of case S2. The grey area is the particle
trajectory range with $R_{f}=0.5,1,2$. The small grey areas indicate that the
predicted particle trajectories in the particle-particle collisions are not
much affected by $\Delta t_{f}$ if it is sufficiently small ($R_{f}\geqslant
0.5$).

Figure 9: Simulation results of particle-particle collisions by the MTSA: (a)
wet coefficient of restitution in normal particle-particle collisions, where
the yellow error bars indicate the prediction uncertainty of $e_{n}$ with
$R_{f}=0.5,1,2$; (b) the predicted trajectories of the particle contact points
with $St=34$, where the grey area is the particle trajectory range with
$R_{f}=0.5,1,2$. $t_{c}=t\mid_{\delta_{n}=0}$ is the instant of collision, and
$t_{ref}=\sqrt{D_{p}/\left|\bm{g}\right|}$ is the reference time scale.

To summarize, by considering the results of particle-wall and particle-
particle collisions together, $R_{i}=4$ is determined to be optimal for
collisions with different $St$ and $D_{p}/\Delta x$ values. The prediction
results by the MTSA with $R_{i}=4$ have been shown to be in good agreement
with the experimental measurements and the benchmark results (case REF) and
made significant improvements over the results without inserting additional
fluid-particle interaction substeps (case S2).

####  4.2.3 Oblique particle-wall collision

Here we assessed the MTSA ($R_{f}=1$, $R_{i}=4$, $R_{m}=40$) for oblique
particle-wall collisions and compared the results with the experimental data
of Joseph and Hunt (2004). The physical parameters in the simulations are
close to those in the experiment of Yang and Hunt (2006): $D_{p}=2.5$ mm,
$\rho_{f}=998$ kg/m3, $\mu_{f}=1$ cP and $e_{n,d}=0.97$, and other parameters
are listed in table 4. The effective angles of incidence and rebounding are
defined as $\Psi_{in}=u_{in,t}/u_{in,n}$ and $\Psi_{out}=u_{out,t}/u_{in,n}$,
respectively, where $u_{in,n}$ and $u_{in,t}$ are the normal and tangential
components of the particle impact velocity, respectively, and $u_{out,t}$ is
the tangential component of the particle rebounding velocity. The
computational domain and grid resolution are the same as those used for the
previous particle-wall collisions in §4.1. Following the work of Costa et al.
(2015), the particle falling velocity was controlled by an oblique
acceleration with the direction vector
$\bm{e}_{g}=-\sin(\phi_{in})\bm{e}_{y}-\cos(\phi_{in})\bm{e}_{x}$ to yield the
desired incidence angle $\Psi_{in}$. The magnitude of the particle
acceleration was set to $g=98.1$ m/s2.

Table 4: Parameters used in the simulation of oblique particle-wall collisions. $\mu_{c,wet}$ is the wet friction coefficient. Material |  $\rho_{p}$ (kg/m3) |  $e_{n,d}$ |  $e_{t,d}$ |  $\mu_{c}$ |  $\mu_{c,wet}$
---|---|---|---|---|---
Steel |  7800 |  0.97 |  0.34 |  0.11 |  0.02
Glass |  2540 |  0.97 |  0.39 |  0.10 |  0.15

Comparison of the normalized incidence and the rebounding angle of the oblique
collision between the simulation results obtained by the MTSA and the
experimental data for steel and glass particles are presented in figure 10. It
is seen that the simulation results by the MTSA are in good agreement with the
experimental data of Joseph and Hunt (2004) in the entire range of incidence
angles.

Figure 10: Dependence between the effective angles of the incidence and
rebounding: (a) without viscous fluid and (b) with viscous fluid.

###  4.3 Simulation of turbulent flow over an erodible sediment bed

Now we applied the MTSA with $R_{f}=0.5,R_{i}=4,R_{m}=40$ to simulate a
turbulent flow over an erodible sediment bed. The turbulent flow was driven by
a constant pressure gradient that is balanced by the shear stress on the
sediment bed. Periodic boundary conditions were imposed in the streamwise and
spanwise boundaries of the computational domain, no-slip boundary conditions
were applied on the bottom surface of the domain and the particle surfaces,
and free-slip boundary conditions were imposed on the top surface of the
domain. The sediment bed was laid on the smooth bottom wall of the channel
under gravity. The sediment bed consisted of $N_{p}=7200$ particles with four
layers, as shown in Figure 11. In the present study, the sediment bed was
generated by a sedimentation simulation for particles settling under gravity
while turning off the hydrodynamic force, similar to Kidanemariam and Uhlmann
(2014b).

Figure 11: Instantaneous snapshot of the turbulent flow over the erodible
sediment bed. The flow field is colored by the value of the non-dimensional
streamwise velocity $u/u_{\tau}$.

For the parameters of the carrier phase, the size of the computational domain
was $L_{x}\times L_{y}\times L_{z}=(6\times 1\times 3)H$ with a uniform
Cartesian grid of $N_{x}\times N_{y}\times N_{z}=1200\times 200\times 600$
grid numbers, where $H$ is the computational domain height in the $y$
direction. Here, $x$, $y$, and $z$ denote the streamwise, wall-normal and
spanwise directions, respectively. The friction Reynolds number is defined as
$Re_{\tau}=u_{\tau}H_{e}/\nu=94$, where $u_{\tau}=\sqrt{\tau_{b}/\rho_{f}}$ is
the friction velocity, $\tau_{b}$ is the mean shear stress at an effective
sediment-bed height $y_{b}$, and $H_{e}=H-y_{b}$ is the effective boundary-
layer height. The bulk Reynolds number $Re_{b}=U_{b}H_{e}/\nu$ are 569, 582,
and 614, for case MTSA, SCTA1, and SCTA2, respectively, where
$U_{b}=(1/H)\int_{0}^{H}\left\langle u\right\rangle dy$ is the bulk velocity.
The effective bed height $y_{b}$ was determined at the height where the mean
particle porosity $\left\langle\phi_{p}\right\rangle=0.1$ Kidanemariam and
Uhlmann (2014a), which yielded $y_{b}=0.35H$. The Shields number is
${\it\Theta}=u_{\tau}^{2}/((\rho_{p}/\rho_{f}-1)gD_{p})=0.12$, which is above
the critical Shields number ${\it\Theta_{c}}=0.047$ Meyer-Peter and Müller
(1948); Wong and Parker (2006), and the corresponding Galileo number is
$Ga=\sqrt{(\rho_{p}/\rho_{f}-1)gD_{p}^{3}}/\nu=42$. The gravitational
acceleration is $g=9.81$ m/s2.

For the parameters of the dispersed phase, the density ratio between the
particle and the fluid is $\rho_{p}/\rho_{f}=2.65$. The dimensionless particle
diameter is $D_{p}/H=0.1$. Referring to the recent particle-resolved sediment
transport simulations Jain et al. (2021); Kidanemariam et al. (2022), the grid
resolution was chosen as $D_{p}/\Delta x=20$, which is sufficient for sediment
transport. The corresponding Eulerian grid size is $\Delta x^{+}=\Delta
y^{+}=\Delta z^{+}=0.47<1$, which is able to fully resolve the turbulent
channel flow. The superscript $+$ indicates quantities normalized in viscous
units (by $u_{\tau}$ and $\nu_{f}$). The particle diameter in wall units was
$D_{p}^{+}=D_{p}/(\nu_{f}/u_{\tau})=14.5$. The mean volume fraction of the
dispersed phase was 21%.

For the parameters in the collision model, the normal restitution coefficient
is $e_{n,d}=0.97$, the tangential restitution coefficient is $e_{t,d}=0.39$,
and the friction coefficient is $\mu_{c}=0.15$; these were chosen according to
the material properties of sand particles Joseph and Hunt (2004).

The simulations were carried out on the Tianhe-2A supercomputer with 50 nodes
(1200 cores). Each node has two Intel Xeon E5-2692 cores and 64 GB of memory.

####  4.3.1 Effect of particle stiffness

Cases MTSA, SCTA1, and SCTA2 with different particle stiffness listed in table
5 were designed to investigate the effect of particle stiffness on the
turbulence and particle statistics of the sediment transport simulation. For
the MTSA, $\Delta t_{f}=T_{c,min}/R_{f}=7.6\times 10^{-6}$ s was chosen to
ensure that it can resolve mostly all collisions with different impact
velocities, where $T_{c,min}$ is the minimum physical collision time obtained
by equation (17) with the maximum impact velocity $u_{in,max}=8u_{\tau}$ based
on a similar particle-resolved simulation by Vowinckel et al. (2014). The
physical stiffness $E=55$ GPa for the glass particle was applied in the MTSA,
while the particle stiffness was reduced in cases SCTA1 and SCTA2 by
stretching the collision time with the stretching coefficient $N=10$ Kempe and
Fröhlich (2012a). In case SCTA2, the particle stiffness was further reduced by
increasing $T_{c}$ compared with case SCTA1.

Table 5: Parameters used in the simulations of the turbulent flow over the erodible sediment bed. The CFL number was calculated by $CFL=\max_{i=1,3}\left|u_{i}\right|\Delta t/\Delta x$. Case |  $\Delta t_{f}$ (s) |  $T_{c}$ (s) |  $E$ (GPa) | CFL | Notes
---|---|---|---|---|---
MTSA | $7.6\times 10^{-6}$ | Eq.(17) | 55 | 0.113 | Benchmark case
SCTA1 | $7.6\times 10^{-6}$ | $10\Delta t_{f}$ | 0.03 | 0.115 | Same $\Delta t_{f}$ as case MTSA
SCTA2 | $7.6\times 10^{-5}$ | $10\Delta t_{f}$ | $10^{-4}$ | 1.21 | Same $N$ as case SCTA1

In our simulations, once the volume flux of the turbulent flow was
statistically stationary, an averaging procedure described in Appendix C was
applied to the flow and particle variables. In total, 250 snapshots of the
flow field and 2500 snapshots of the particles during a time period of
$15H/u_{\tau}$ were collected for time averaging. In order to distinguish
whether the error came from reducing particle stiffness or statistical
uncertainty, the time period $15H/u_{\tau}$ was equally divided into two
periods of $7.5H/u_{\tau}$. The results averaged over these two time periods
were defined as the lower and upper boundaries (statistical uncertainty range)
of the simulation, as shown by the shaded region in figures 12 to 14. The
lines/symbols in the following figures are the results averaged over a total
time period of $15H/u_{\tau}$. In the following, the turbulence and particle
statistics are compared between cases MTSA, SCTA1 and SCTA2. The wall-normal
coordinate $Y=y-y_{b}$ is adopted and normalized by $H_{e}$.

Figure 12 shows the profiles of the mean particle porosity (planar density)
obtained in each $x-z$ plane. The operator $\left\langle\cdot\right\rangle$
indicates an average over $x-z$ plane (not a slab) and time. With this
averaging operation, the theoretical maximum limit of particle porosity in
dense packing is 0.906 (the largest planar density in the face-centered cubic
structure), which is much larger than 0.74 of the largest particle volume
fraction in three-dimensional space. Following Kidanemariam and Uhlmann
(2014b), the mean particle porosity was calculated by equation (C.4) in
Appendix C. The four peaks under the grey line indicate that the sediment bed
consists of four layers of particles. The peak heights correspond to the mean
center height of each layer of particles Vowinckel et al. (2014). It should be
noted that the particle porosity Kidanemariam and Uhlmann (2014b); Vowinckel
et al. (2014), which is the particle fraction in the $x-z$ plane rather than
the volume fraction Rao and Capecelatro (2019), was employed since the
particles are larger than the grid. If we chose a $y$-slab with a height
$y_{s}$ larger than $D_{p}$, the mean volume fraction $\Phi$ can be calculated
through the integration of the particle porosity as
$\Phi=(1/y_{s})\int_{y_{0}}^{y_{0}+y_{s}}\left\langle\phi_{p}\right\rangle
dy$, where $y_{0}$ is the lower boundary of the $y$-slab. The mean volume
fraction of the entire sediment bed is 0.58 with $y_{0}=-y_{b}/H_{e}$ and
$y_{s}=y_{b}/H_{e}$ in the present simulation, which is lower than the
theoretical maximum limit of dense packing of 0.74 Rao and Capecelatro (2019).
The maximum particle porosity is 0.832, which is also lower than the
theoretical maximum limit of 0.906 for dense packing. As shown in figure 12,
the particle porosity above $Y/H_{e}=0.15$ is approximately 0.01, which
indicates that particles seldom jumped above $Y/H_{e}=0.15$. Therefore,
particle samples above $Y/H_{e}=0.15$ within the averaging time window may not
be sufficiently large to obtain convergent statistics. This problem can also
be seen in figures 9 and 10 of Ji et al. (2014). Therefore, we only focus on
the near-wall region ($0<Y/H_{e}<0.15$) in the following analysis.

Figure 12: Profiles of the mean particle porosity. The grey line indicates the
effective height of the sediment bed. The shaded regions of different colors
indicate the ranges of the statistical errors in different cases.

For comparing the present simulation results with available experimental or
numerical results, we cannot find any particle-resolved simulation studies or
experimental measurements under similar conditions. In many experiments, both
the size of particles and the Reynolds number of turbulent flow are much
larger than the present one. In addition, as the particle volume fraction is
usually very high near a sediment bed, it is very challenging to distinguish
tracers and inertial particles in experiment Brandt and Coletti (2022). As a
result, near-wall turbulence statistics are almost impossible to be extracted.
We cannot find any experimental results of near-sediment-bed turbulence
statistics in the literature also. A measurable quantity of sediment transport
rate can be compared, which is defined as $q=(\pi
D_{p}^{3}/6)\textstyle\left(\sum_{N_{p}}\langle{u}_{p}\rangle\right)/\left(L_{x}L_{y}L_{z}\sqrt{(\rho_{p}/\rho_{f}-1)gD_{p}^{3}}\right)$
in a non-dimensional form, where $\langle{u}_{p}\rangle$ is mean particle
velocity. The non-dimensional sediment transport rates in cases MTSA, SCTA1,
and SCTA2 are $q_{MTSA}=0.0461$, $q_{SCTA1}=0.0568$, and $q_{SCTA2}=0.1396$,
respectively. The well-known empirical formula of Wong and Parker (2006) gives
$q^{*}=4.93\left(\it\Theta-\it\Theta_{c}\right)^{1.6}$, where $\it\Theta_{c}$
is the critical Shields number. Following Wong and Parker (2006), if
${\it\Theta_{c}}=0.047$ is employed, we can get $q^{*}=0.0748$. It seems that
the prediction by SCTA1 is the most accurate, but we argue that the formula of
Wong and Parker (2006) was fitted with data at much larger flow-submerge
($Y/D_{p}\sim 40$) and flow Reynolds numbers. The particle stability increases
with decreasing flow depth under constant shear velocity, which is likely
related to the suppression of the energetic large-scale turbulent motions due
to the limited separation between flow depth and roughness length scales
(smaller flow Reynolds number) Cameron et al. (2020). Therefore, we believe
that the actual sediment transport rate under the present simulation condition
should be much smaller than $q^{*}=0.0748$.

In the near-wall region, the statistical uncertainties of the particle
porosity are negligible compared with the errors caused by reducing particle
stiffness. The present results show that when particle stiffness is reduced,
the mean particle porosity is underestimated, and the underestimation
increases as particle stiffness decreased, as displayed in figure 12. The
relative difference between the results of case SCTA1 and case MTSA and
between case SCTA2 and case MTSA are 24% and 40% at $Y/H_{e}=0.02$,
respectively. Thus, fewer particles are entrained in cases SCTA1 and SCTA2.
This may be caused by excluding the hydrodynamic force during the collisions,
which significantly ignores the effect of turbulence on particle entrainment
with a time scale smaller than $T_{c}$. Vowinckel et al. (2016) found that
96.5% of the entrained particles were triggered by particle collisions.
Therefore, the effects of turbulence may be substantially underestimated if
$T_{c}$ is greatly stretched, and fewer particles would be entrained in the
near-wall region.

Figure 13 (a) shows the streamwise mean flow velocity profile. In the near-
wall region, the statistical uncertainties of the streamwise mean flow
velocity are negligible. The results show that when particle stiffness is
reduced by the SCTA, the streamwise mean flow velocity is overestimated
compared with the MTSA, and the difference increases as particle stiffness
decreases. The relative difference between the results of case SCTA1 and case
MTSA and the results of case SCTA2 and case MTSA are 6% and 35% at
$Y/H_{e}=0.02$, respectively. This is probably because fewer particles are
entrained from the bed in cases SCTA1 and SCTA2, which exert a smaller drag
force on the flow. The profiles of the root-mean-squared (r.m.s) fluctuating
flow velocities in the three directions are shown in figure 13 (b). In the
near-wall region, it can be seen that the wall-normal and spanwise fluctuating
flow velocities are close between cases SCTA1, SCTA2 and MTSA. For the
streamwise fluctuating flow velocity, the statistical uncertainty is smaller
than the uncertainty by the collision model. In fact, the r.m.s streamwise
fluctuating flow velocity is underestimated by the SCTA compared with the
MTSA, and the underestimation increases as particle stiffness decreases. This
may be because particles are more difficult to be entrained in cases SCTA1 and
SCTA2, which imposes fewer disturbances to the fluid in the near-wall region.
Figure 13 (c) shows the profiles of the total shear stress
$\left(-\left\langle u^{\prime}v^{\prime}\right\rangle+\nu_{f}{\rm
d}\langle{u}\rangle/{\rm d}y\right)/u_{\tau}^{2}$, where $u^{\prime}$ and
$v^{\prime}$ are the streamwise and wall-normal flow velocity fluctuations,
respectively. The total shear stress nicely follows a linear profile in the
outer flow indicating that the turbulent flow has been fully developed. The
total shear stress is underestimated by the SCTA compared with the MTSA in the
near-wall region, and the underestimation increases as particle stiffness
decreases.

Figure 13: Turbulence statistics in sediment transport simulations: (a) the
streamwise mean flow velocity and (b) the root-mean-squared fluctuating flow
velocities, and (c) the total shear stress $\left(-\left\langle
u^{\prime}v^{\prime}\right\rangle+\nu_{f}{\rm d}\langle{u}\rangle/{\rm
d}y\right)/u_{\tau}^{2}$. The grey line indicates the effective height of the
sediment bed. The grey dash line indicates $1-Y/H_{e}$. The shaded regions
indicate the ranges of statistical uncertainties.

Figure 14 (a) shows the profiles of the mean particle velocities. It is seen
that the profiles of the wall-normal and spanwise mean particle velocities are
quite close between cases SCTA1, SCTA2 and MTSA. In the near-wall region, the
statistical uncertainties of the streamwise mean particle velocity are much
smaller than those by reducing particle stiffness. This result definitely
shows that when particle stiffness is artificially reduced, the streamwise
mean particle velocity will be overestimated, and the difference increases as
particle stiffness decreases. The relative difference of the streamwise mean
particle velocity between cases SCTA1 and MTSA and between cases SCTA2 and
MTSA can reach up to 87% and 325% at $Y/H_{e}=0.02$, respectively. This may be
because the mean flow velocity of cases SCTA1 and SCTA2 are higher than that
of case MTSA, thus the drag force acting on the entrained particles increases
to accelerate particles to higher speeds in cases SCTA1 and SCTA2. Figure 14
(b) displays the profiles of the mean particle angular velocities. It can be
seen that the profiles of the streamwise and wall-normal mean particle angular
velocities are close between cases SCTA1, SCTA2 and MTSA. In the near-wall
region, the statistical uncertainties of the spanwise mean particle angular
velocity are much smaller than those by reducing particle stiffness. The
results in the near-wall region demonstrate that when particle stiffness is
reduced by the SCTA, the spanwise mean particle angular velocity may be
overestimated, and the difference increases as particle stiffness decreases.
The relative difference of the spanwise mean particle angular velocity between
cases SCTA1 and MTSA and between cases SCTA2 and MTSA can reach up to 62% and
150% at $Y/H_{e}=0.02$, respectively. This may be because the mean flow
velocity gradient is larger in cases SCTA1 and SCTA2.

It should be noted that a version of SCTA similar to Costa et al. (2015) is
adopted here to be compared with the MTSA, which incorporates a two-parameter
lubrication model, a linear spring-dashpot system, and substeps for particle
collisions and motions. It is found that there are evident differences between
the results of the sediment transport case using the SCTA and MTSA for the
statistics especially particle statistics if the particle stiffness is greatly
reduced. The differences increase as particle stiffness decreases (stretched
collision time increases). Using different versions of SCTA with different
model details like the lubrication model, linear vs. non-linear spring-
dashpot, time integration among others may affect the results. However, we
believe that the assumption of excluding the hydrodynamic force during
particle collisions and the much longer collision time in the SCTA are the
major reasons for the differences between the results using the SCTA and MTSA,
as the other modeling and numerical aspects are the same in our code.

Figure 14: Particle statistics in sediment transport simulations: (a) the mean
particle velocity and (b) the mean particle angular velocity. The grey line
indicates the effective height of the sediment bed. The shaded regions
indicate the ranges of the statistical uncertainties.

####  4.3.2 The computational cost of the MTSA

The computational cost of a complete simulation of sediment transport using
the traditional soft-sphere model with a fine time step is very high.
Therefore, we only ran 1000 time steps in case REF to obtain its computation
performance. Here we compare the ratio $(T/\Delta t_{f})/(T_{REF}/\Delta
t_{f,REF})$ between different algorithms, which represents the computational
time ratio for advancing a unit flow time in the simulations, where $T$ is the
elapsed clock time spent over 1000 flow time steps. Furthermore, $\Delta
t_{f}/\Delta t_{f,REF}=16$ (_i.e._ , $R_{f}=0.5$) was used in the following
simulations, and the computational time ratios of $\Delta t_{f}/\Delta
t_{f,REF}=4$ and 8 were also tested. As shown in figure 15, the computational
time ratio of the MTSA is close to the ideal ratio $\Delta t_{f,REF}/\Delta
t_{f}$. The computational efficiency can be increased by an order of magnitude
using the MTSA with $R_{f}=0.5,R_{i}=4,R_{m}=40$ compared with the traditional
soft-sphere model. This is because the pressure Poisson equation is not solved
in the inserted substeps, which significantly reduces the computational cost.
It is also seen that the computational efficiency of the MTSA is not degraded
too much with that of the SCTA.

Figure 15: Variation of the computational time ratio with the time step ratio.
The blue line is the ideal computational time ratio, which equals $\Delta
t_{f,REF}/\Delta t_{f}$.

##  5 Conclusions

In the present study, a multiple-time-step integration algorithm (MTSA) has
been proposed for particle-resolved simulations with physical collision time.
The MTSA adopts three different time steps to resolve fluid flow, fluid-
particle interaction and particle collision. We have successfully demonstrated
the accuracy and convergence of the MTSA in particle-wall and particle-
particle collision problems. It is shown that the MTSA can yield excellent
agreement with experiments and benchmark simulations using the physical
collision time and on a coarse grid ($D_{p}/\Delta x=O(10)$). It also
overcomes the main drawback of a version of SCTA similar to Costa et al.
(2015), in which an artificially stretched collision time is introduced, which
may lead to the collision time being nonphysical and an underestimation of the
hydrodynamic force and particle entrainment in sediment transport.

The impact of the parameters $R_{f}$, $R_{i}$, and $R_{m}$ in the MTSA has
been systematically analyzed and discussed. We found that $R_{i}=4$ and
$R_{m}=40$ are universally suitable for particle-wall and particle-particle
collisions over a wide range of particle Stokes numbers. The absolute error of
$e_{n}$ can be controlled to be within about 0.03 for different $R_{f}$
values. After determining the optimal parameters of the MTSA, we applied the
MTSA to simulate a turbulent flow over an erodible sediment bed. The
statistics of turbulence and particles by the MTSA and a version of SCTA
similar to Costa et al. (2015) were compared to investigate the effect of
particle stiffness. The results demonstrated that if particle stiffness is
artificially reduced by the SCTA, the number of entrained particles over the
erodible sediment bed will be evidently underestimated compared with those by
the MTSA, which will mainly impact the streamwise mean flow velocity,
streamwise fluctuating flow velocity, streamwise mean particle velocity and
spanwise mean particle angular velocity in the near-wall region. The
prediction differences increase as particle stiffness decreases. Furthermore,
the computation efficiencies between the MTSA and the traditional soft-sphere
model were compared. The results indicated that the MTSA can substantially
reduce the computational cost by up to one order of magnitude, and make the
computational cost acceptable for particle-resolved simulations with physical
collision time. As such, the proposed MTSA can be a valuable algorithm for
generating high-fidelity and accurate data efficiently in particle-resolved
simulations.

For future studies, we plan to perform more intensive particle-resolved
simulations of sediment transport, including higher Reynolds numbers and
particle/fluid density ratios. For example, the Reynolds number $Re_{\tau}$ is
an important non-dimensional parameter that may also affect the prediction
difference between the simulations using the SCTA and MTSA, thus simulations
at higher $Re_{\tau}$ should be valuable. Simulating sediment transport using
the adaptive mesh refinement (AMR) technique to reduce computational costs
Zeng et al. (2022a, b) is also within our future plans.

##  6 Acknowledgments

This work was financially supported by grants from the National Natural
Science Foundation of China (Nos. 92052202 and 11972175) and the National
Numerical Windtunnel project. The computation was performed at the Tianhe-2A
supercomputer of the National Supercomputer Center in Guangzhou. The authors
report no conflict of interest.

## Appendix A. Definitions for particle-particle and particle-wall collision

The variables used in particle collision are defined below. Some definitions
depend on whether the interaction is between a particle $p$ and a wall (P-W)
or between a particle $p$ and a particle $q$ (P-P). The variables are defined
as follows:

$\bm{n}$ \- normal unit vector of contact

$\displaystyle\bm{n}=\frac{\bm{x}_{q}-\bm{x}_{p}}{\left|\bm{x}_{q}-\bm{x}_{p}\right|},\quad$
$\displaystyle\text{(P-P)},$ (A.1)
$\displaystyle\bm{n}=\frac{\bm{x}_{w}-\bm{x}_{p}}{\left|\bm{x}_{w}-\bm{x}_{p}\right|},$
$\displaystyle\text{(P-W)},$ (A.2)

$\delta_{n}$ \- distance between two surfaces

$\displaystyle\delta_{n}=\left|\bm{x}_{q}-\bm{x}_{p}\right|-R_{p}-R_{q},\quad$
$\displaystyle\text{(P-P)},$ (A.3)
$\displaystyle\delta_{n}=\left|\bm{x}_{w}-\bm{x}_{p}\right|-R_{p},$
$\displaystyle\text{(P-W)},$ (A.4)

$\bm{u}_{cp}$ \- relative velocity of contact point

$\displaystyle\bm{u}_{cp}=\bm{u}_{p}-\bm{u}_{q}+R_{p}\bm{\omega}_{p}\times\bm{n}+R_{q}\bm{\omega}_{q}\times\bm{n},\quad$
$\displaystyle\text{(P-P)},$ (A.5)
$\displaystyle\bm{u}_{cp}=\bm{u}_{p}+R_{p}\bm{\omega}_{p}\times\bm{n},$
$\displaystyle\text{(P-W)},$ (A.6)

$\bm{u}_{cp,n}$ \- normal component of $\bm{u}_{cp}$

$\displaystyle\bm{u}_{cp,n}=(\bm{u}_{cp}\cdot\bm{n})\bm{n},$ (A.7)

$\bm{u}_{cp,t}$ \- tangential component of $\bm{u}_{cp}$

$\displaystyle\bm{u}_{cp,t}=\bm{u}_{cp}-\bm{u}_{cp,n},$ (A.8)

and $\bm{\delta}_{t}$ is the tangential displacement of the contact point.

The direction of the tangential unit vector changes in every time step.
Therefore, we need to rotate the displacement from the previous time step onto
a plane tangent to $\bm{n}$. Then $\bm{\delta}_{t}$ is calculated to be the
same as Biegert et al. (2017):

$\displaystyle\widetilde{\bm{\delta}_{t}}=\bm{\delta}_{t}^{k-1}-(\bm{\delta}_{t}^{k-1}\cdot\bm{n})\bm{n},$
(A.9)
$\displaystyle\widehat{\bm{\delta}_{t}}=\frac{\left|\bm{\delta}_{t}^{k-1}\right|}{\left|\widetilde{\bm{\delta}_{t}}\right|}\widetilde{\bm{\delta}_{t}},$
(A.10) $\displaystyle\bm{\delta}_{t}^{k}=\widehat{\bm{\delta}_{t}}+\Delta
t\bm{u}_{cp,t}.$ (A.11)

## Appendix B. Discretization schemes

The discretization schemes adopted in the MTSA are as follows: an explicit
second-order Runge–Kutta (RK2) method is used for fluid flow advancement Yang
et al. (2017, 2018); Cui et al. (2018); He et al. (2022). At each substep of
the RK2 method, the fractional-step method of Kim and Moin (1985) is applied
to ensure that the flow field is divergence-free. In order to insert enough
substeps to ensure that the flow can adapt to the change of the particle
motion during collisions and reduce the computational cost, the explicit Euler
method is adopted in the substeps for time advancement. The fluid and particle
are coupled by the multidirect forcing IB method Luo et al. (2007); Breugem
(2012). The discretized equations of fluid flow and particle motion in the
$k$th Runge–Kutta step are given as follows:

$\displaystyle\text{do}\ k=1,2$ $\displaystyle\ \
\hat{\bm{u}}^{k-1,0}={\bm{u}}^{k-1},$ $\displaystyle\ \ {\color[rgb]{1,0,0}do\
i=1,N_{i}}$ $\displaystyle\ \ \ \
\hat{\bm{u}}^{k-1,i}=\hat{\bm{u}}^{k-1,i-1}+{\color[rgb]{1,0,0}\Delta
t_{i}}\left[\alpha_{k}\hat{\bm{H}}^{k-1,i-1}-\beta_{k}\left(\hat{\bm{H}}^{k-2,i-1}-\frac{1}{\rho}\nabla
p^{k-2}\right)\right],$ (B.1) $\displaystyle\ \ \ \
\bm{u}^{*,0}=\hat{\bm{u}}^{k-1,i},$ $\displaystyle\ \ \ \
\bm{F}_{p,l}^{k-1,i,0}=0,$ $\displaystyle\ \ \ \ \text{do}\ s=1,N_{s}$
$\displaystyle\ \ \ \ \ \
\bm{U}_{l}^{*,s-1}=\sum_{ijk}\bm{u}_{ijk}^{*,s-1}\delta_{d}\left(\bm{x}_{ijk}-\bm{X}_{l}^{k-1,i-1,N_{m}}\right)\Delta
x\Delta y\Delta z,$ (B.2) $\displaystyle\ \ \ \ \ \
\bm{F}_{p,l}^{*,s-1}=\frac{\bm{U}_{p}{\left(\bm{X}_{l}^{k-1,i-1,N_{m}}\right)-\bm{U}_{l}^{*,s-1}}}{{\color[rgb]{1,0,0}\Delta
t_{i}}},$ (B.3) $\displaystyle\ \ \ \ \ \
\bm{f}_{ijk}^{*,s-1}=\sum_{l=1}^{N_{l}}\bm{F}_{p,l}^{*,s-1}\delta_{d}\left(\bm{x}_{ijk}-\bm{X}_{l}^{k-1,i-1,N_{m}}\right)\Delta
V_{l},$ (B.4) $\displaystyle\ \ \ \ \ \
\bm{u}^{*,s}=\bm{u}^{*,s-1}+\left(\alpha_{k}-\beta_{k}\right){\color[rgb]{1,0,0}\Delta
t_{i}}\bm{f}^{*,s-1},$ (B.5) $\displaystyle\ \ \ \ \ \
\bm{F}_{p,l}^{k-1,i,s}=\bm{F}_{p,l}^{k-1,i,s-1}+\left(\alpha_{k}-\beta_{k}\right)\bm{F}_{p,l}^{*,s-1},$
(B.6) enddo $\displaystyle\ \ \ \ \hat{\bm{u}}^{k-1,i}=\bm{u}^{*,N_{s}},$
$\displaystyle\ \ \ \ \bm{u}_{p}^{k-1,i,0}=\bm{u}_{p}^{k-1,i-1,N_{m}},$
$\displaystyle\ \ \ \ \bm{x}_{p}^{k-1,i,0}=\bm{x}_{p}^{k-1,i-1,N_{m}},$
$\displaystyle\ \ \ \
\bm{\omega}_{p}^{k-1,i,0}=\bm{\omega}_{p}^{k-1,i-1,N_{m}},$ $\displaystyle\ \
\ \ {\color[rgb]{0,0,1}do\ p=1,N_{m}}$ $\displaystyle\ \ \ \ \ \
\bm{u}_{p}^{k-1,i,p}=\bm{u}_{p}^{k-1,i,p-1}+\left(\alpha_{k}-\beta_{k}\right)\Bigg{\\{}-\frac{{\color[rgb]{0,0,1}\Delta
t_{p}}}{V_{p}}\frac{\rho_{f}}{\rho_{p}}\sum_{l=1}^{N_{l}}\bm{F}_{p,l}^{k-1,i,N_{s}}\Delta
V_{l}+$ $\displaystyle\ \ \ \ \ \
\frac{1}{V_{p}}\frac{\rho_{f}}{\rho_{p}}\left(\left\\{\int_{V_{p}}\bm{u}dV\right\\}^{N}-\left\\{\int_{V_{p}}\bm{u}dV\right\\}^{N-1}\right)+{\color[rgb]{0,0,1}\Delta
t_{p}}\left(1-\frac{\rho_{f}}{\rho_{p}}\right)\bm{g}+$ $\displaystyle\ \ \ \ \
\ {\color[rgb]{0,0,1}\Delta
t_{p}}\frac{\bm{F}_{p,lub}^{k-1,i,p-1}}{\rho_{p}V_{p}}+{\color[rgb]{0,0,1}\Delta
t_{p}}\frac{\bm{F}_{p,col}^{k-1,i,p-1}}{\rho_{p}V_{p}},\Bigg{\\}}$ (B.7)
$\displaystyle\ \ \ \ \ \
\bm{x}_{p}^{k-1,i,p}=\bm{x}_{p}^{k-1,i,p-1}+\left(\alpha_{k}-\beta_{k}\right){\color[rgb]{0,0,1}\Delta
t_{p}}\bm{u}_{p}^{k-1,i,p},$ (B.8) $\displaystyle\ \ \ \ \ \
\bm{\omega}_{p}^{k-1,i,p}=\bm{\omega}_{p}^{k-1,i,p-1}+\left(\alpha_{k}-\beta_{k}\right)\Bigg{\\{}-{\color[rgb]{0,0,1}\Delta
t_{p}}\frac{\rho_{f}}{I_{p}}\sum_{l=1}^{N_{l}}\bm{r}_{l}\times\bm{F}_{p,l}^{k-1,i,N_{s}}\Delta
V_{l}+$ $\displaystyle\ \ \ \ \ \
\frac{\rho_{f}}{I_{p}}\left(\left\\{\int_{V_{p}}\bm{r}\times\bm{u}dV\right\\}^{N}-\left\\{\int_{V_{p}}\bm{r}\times\bm{u}dV\right\\}^{N-1}\right)+{\color[rgb]{0,0,1}\Delta
t_{p}}\frac{\bm{T}_{p,col}^{k-1,i,p-1}}{I_{p}}\Bigg{\\}},$ (B.9)
$\displaystyle\ \ \ \ {\color[rgb]{0,0,1}enddo}$ $\displaystyle\ \ \ \
\bm{U}_{p}\left(\bm{X}_{l}^{k-1,i,N_{m}}\right)=\bm{u}_{p}^{k-1,i,N_{m}}+\bm{\omega}_{p}^{k-1,i,N_{m}}\times(\bm{X}_{l}^{k-1,i,N_{m}}-\bm{x}_{p}^{k-1,i,N_{m}}),$
(B.10) $\displaystyle\ \ {\color[rgb]{1,0,0}enddo}$ $\displaystyle\ \
\nabla^{2}p^{k-1}=\frac{\rho_{f}}{\alpha_{k}\Delta
t_{f}}\nabla\cdot\hat{\bm{u}}^{k-1,N_{i}},$ (B.11) $\displaystyle\ \
\bm{u}^{k}=\hat{\bm{u}}^{k-1,N_{i}}-\frac{\alpha_{k}\Delta
t_{f}}{\rho_{f}}\nabla p^{k-1},$ (B.12) enddo

where the superscripts $N$, $k$, $i$, $p$, and $s$ are the indices of the flow
time step, the Runge–Kutta substep, the fluid–particle interaction substep,
the particle motion substep, and the IB method iteration substep,
respectively. The coefficients in the RK2 scheme are $\alpha_{1}=1$,
$\beta_{1}=0$, and $\alpha_{2}=\beta_{2}=0.5$. $N_{s}$ is the number of
iterations to calculate the IB force in the multidirect forcing immersed
boundary methodLuo et al. (2007); Breugem (2012). The idea of this method is
to iteratively determine the IB force to make the fluid velocity modified by
the IB force better satisfy the no-slip boundary condition on the particle
surface. $N_{s}=3$ is employed following the suggestion of Breugem (2012). The
subscripts $ijk$ refer to the quantities on an Eulerian grid with indices
$(i,j,k)$, and the subscript $l$ refers to the quantity on a Lagrangian point
with index $l$. $\hat{\bm{u}}$ is the fluid velocity in f fluid-particle
interaction substep,
$\hat{\bm{H}}=-\nabla\cdot\left(\bm{\hat{u}\hat{u}}\right)+\nu_{f}\nabla^{2}\bm{\hat{u}}$
is the sum of the convection and viscous terms. $\bm{u}^{*}$, $\bm{f}^{*}$ and
$\bm{x}$ are the intermediate fluid velocity, volume force, and coordinate on
the Eulerian grid, respectively, and $\bm{U}^{*}$, $\bm{F}^{*}$, $\bm{X}$ are
the corresponding quantities on the Lagrangian grid. The quantities between
the Eulerian and Lagrangian grid points are transferred through a regularized
Dirac delta function $\delta_{d}$. Here, we use the three-point regularized
Dirac delta function of Roma et al. (1999). The volume integrals
$\int_{V_{p}}{\bm{u}dV}$ and $\int_{V_{p}}{\bm{r}\times\bm{u}dV}$ are
numerically calculated according to Kempe and Fröhlich (2012b). $\Delta V_{l}$
is the volume of the Lagrangian grid cell, and $\bm{U}_{p}$ is the velocity on
the particle surface.

## Appendix C. Averaging for flow and particle variables

### C.1. Averaging for flow variables

Before averaging the flow variables, an indicator function needs to be defined
as $\phi_{f}(\bm{x},t)$ to distinguish an Eulerian grid point at a position
$\bm{x}$ that is inside or outside of a particle, following Kidanemariam and
Uhlmann (2014b), which is

$\phi_{f}(\bm{x},t)=\begin{cases}1,&\text{ if }\bm{x}\ \text{is outside a
particle,}\\\ 0,&\text{ otherwise.}\end{cases}$ (C.1)

Based on the indicator function $\phi_{f}(\bm{x},t)$, only the flow data
outside of the particle are accounted for averaging as follows:

$n_{f}(y_{j})=\sum\limits_{n=1}^{N_{t}}\sum\limits_{i=1}^{N_{x}}\sum\limits_{k=1}^{N_{z}}\phi_{f}(\bm{x}_{ijk},t^{n}),$
(C.2)

where $n_{f}(y_{j})$ is the total number of grids in the $x-z$ plane over
$N_{t}$ time steps for the flow statistics at a given height $y_{j}$.
Therefore, the ensemble average of the flow variables $\xi_{f}(\bm{x},t)$ can
be defined as

$\left\langle\xi_{f}\right\rangle(y_{j})=\frac{1}{n_{f}(y_{j})}\sum\limits_{n=1}^{N_{t}}\sum\limits_{i=1}^{N_{x}}\sum\limits_{k=1}^{N_{z}}\xi_{f}(\bm{x}_{ijk},t^{n})\phi_{f}(\bm{x}_{ijk},t^{n}),$
(C.3)

where the operator $\left\langle\cdot\right\rangle$ indicates the average over
an $x-z$ plane and time. The particle porosity
$\left\langle\phi_{p}\right\rangle(y_{j})$ can be deduced from
$\phi_{f}(\bm{x},t)$:

$\left\langle\phi_{p}\right\rangle(y_{j})=\frac{1}{N_{t}N_{x}N_{z}}\sum\limits_{n=1}^{N_{t}}\sum\limits_{i=1}^{N_{x}}\sum\limits_{k=1}^{N_{z}}(1-\phi_{f}(\bm{x}_{ijk},t^{n})),$
(C.4)

### C.2. Averaging for particle variables

Particle variables are averaged over all particles within each slab. A slab is
generated by dividing $H$ by a thickness of $\Delta h$. An indicator function
$\varphi_{p}^{j}(y_{p},t)$ is defined to distinguish the center height $y_{p}$
of a particle inside or outside of a slab with index $j$ as follows:

$\varphi_{p}^{j}(y_{p},t)=\begin{cases}1,&\text{ if }(j-1)\Delta h\leqslant
y<j\Delta h,\\\ 0,&\text{ otherwise.}\end{cases}$ (C.5)

Based on the indicator function $\varphi_{p}^{j}(y_{p},t)$, the particle
number in each slab can be calculated as

$n_{p}^{j}=\sum\limits_{n=1}^{N_{t}}\sum\limits_{l=1}^{N_{p}}\varphi_{p}^{j}(y_{p}^{l},t^{n}),$
(C.6)

where $n_{p}^{j}$ is the total particle number in the $j$ slab over $N_{t}$
time steps. Therefore, the averaged particle variable $\xi_{p}$ can be defined
as

$\left\langle\xi_{p}(y^{j})\right\rangle=\frac{1}{n_{p}^{j}}\sum\limits_{n=1}^{N_{t}}\sum\limits_{l=1}^{N_{p}}\varphi_{p}^{j}(y_{p}^{l},t^{n})\xi_{p}^{l}(t^{n}).$
(C.7)

A slab thickness of $\Delta h=D_{p}/4$ was selected in the present study.

## References

  * Akiki et al. (2017) Akiki, G., Jackson, T., Balachandar, S., 2017. Pairwise interaction extended point-particle model for a random array of monodisperse spheres. Journal of Fluid Mechanics 813, 882–928.
  * Ardekani et al. (2016) Ardekani, M. N., Costa, P., Breugem, W. P., Brandt, L., 2016. Numerical study of the sedimentation of spheroidal particles. International Journal of Multiphase Flow 87, 16–34.
  * Bagchi and Balachandar (2003) Bagchi, P., Balachandar, S., 2003. Effect of turbulence on the drag and lift of a particle. Physics of fluids 15 (11), 3496–3513.
  * Balachandar and Eaton (2010) Balachandar, S., Eaton, J. K., 2010. Turbulent dispersed multiphase flow. Annual Review of Fluid Mechanics 42, 111–133.
  * Biegert et al. (2017) Biegert, E., Vowinckel, B., Meiburg, E., 2017. A collision model for grain-resolving simulations of flows over dense, mobile, polydisperse granular sediment beds. Journal of Computational Physics 340, 105–127.
  * Brandt and Coletti (2022) Brandt, L., Coletti, F., 2022. Particle-laden turbulence: progress and perspectives. Annual Review of Fluid Mechanics 54, 159–189.
  * Brenner (1961) Brenner, H., 1961. The slow motion of a sphere through a viscous fluid towards a plane surface. Chemical Engineering Science 16 (3-4), 242–251.
  * Breugem (2010) Breugem, W.-P., 2010. A combined soft-sphere collision/immersed boundary method for resolved simulations of particulate flows. In: Fluids Engineering Division Summer Meeting. No. FEDSM-ICNMM2010-30634. pp. 2381–2392.
  * Breugem (2012) Breugem, W.-P., 2012. A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. Journal of Computational Physics 231 (13), 4469–4498.
  * Cameron et al. (2020) Cameron, S., Nikora, V., Witz, M., 2020. Entrainment of sediment particles by very large-scale motions. Journal of Fluid Mechanics 888, A7.
  * Capecelatro and Desjardins (2013) Capecelatro, J., Desjardins, O., 2013. An Euler–Lagrange strategy for simulating particle-laden flows. Journal of Computational Physics 238, 1–31.
  * Cooley and O’Neill (1969) Cooley, M., O’Neill, M., 1969. On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere. Mathematika 16 (1), 37–49.
  * Costa et al. (2015) Costa, P., Boersma, B. J., Westerweel, J., Breugem, W.-P., 2015. Collision model for fully resolved simulations of flows laden with finite-size particles. Physical Review E 92 (5), 053012.
  * Costa et al. (2020) Costa, P., Brandt, L., Picano, F., 2020. Interface-resolved simulations of small inertial particles in turbulent channel flow. Journal of Fluid Mechanics 883, A54.
  * Costa et al. (2021) Costa, P., Brandt, L., Picano, F., 2021. Near-wall turbulence modulation by small inertial particles. Journal of Fluid Mechanics 922, A9.
  * Cox and Brenner (1967) Cox, R. G., Brenner, H., 1967. The slow motion of a sphere through a viscous fluid towards a plane surface—ii small gap widths, including inertial effects. Chemical Engineering Science 22 (12), 1753–1777.
  * Cui et al. (2018) Cui, Z., Yang, Z., Jiang, H.-Z., Huang, W.-X., Shen, L., 2018. A sharp-interface immersed boundary method for simulating incompressible flows with arbitrarily deforming smooth boundaries. International Journal of Computational Methods 15 (01), 1750080.
  * Darmana et al. (2006) Darmana, D., Deen, N. G., Kuipers, J., 2006. Parallelization of an Euler–Lagrange model using mixed domain decomposition and a mirror domain technique: Application to dispersed gas–liquid two-phase flow. Journal of Computational Physics 220 (1), 216–248.
  * Deen et al. (2007) Deen, N., Annaland, M. V. S., Van der Hoef, M. A., Kuipers, J., 2007. Review of discrete particle modeling of fluidized beds. Chemical Engineering Science 62 (1-2), 28–44.
  * Eaton (2009) Eaton, J. K., 2009. Two-way coupled turbulence simulations of gas-particle flows using point-particle tracking. International Journal of Multiphase Flow 35 (9), 792–800.
  * Feng et al. (2010) Feng, Z.-G., Michaelides, E. E., Mao, S., 2010. A three-dimensional resolved discrete particle method for studying particle-wall collision in a viscous fluid. Journal of Fluids Engineering 132 (9), 091302.
  * Finn et al. (2016) Finn, J. R., Li, M., Apte, S. V., 2016. Particle based modelling and simulation of natural sand dynamics in the wave bottom boundary layer. Journal of Fluid Mechanics 796, 340–385.
  * Glowinski et al. (1999) Glowinski, R., Pan, T.-W., Hesla, T. I., Joseph, D. D., 1999. A distributed Lagrange multiplier/fictitious domain method for particulate flows. International Journal of Multiphase Flow 25 (5), 755–794.
  * Glowinski et al. (2001) Glowinski, R., Pan, T.-W., Hesla, T. I., Joseph, D. D., Periaux, J., 2001. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow. Journal of Computational Physics 169 (2), 363–426.
  * Gondret et al. (2002) Gondret, P., Lance, M., Petit, L., 2002. Bouncing motion of spherical particles in fluids. Physics of Fluids 14 (2), 643–652.
  * He et al. (2022) He, S., Yang, Z., Sotiropoulos, F., Shen, L., 2022. Numerical simulation of interaction between multiphase flows and thin flexible structures. Journal of Computational Physics 448, 110691.
  * Hertz (1882) Hertz, H., 1882. Über die berührung fester elastischer körper. Journal für die reine und angewandte Mathematik 92 (22), 156–171.
  * Homann et al. (2013) Homann, H., Bec, J., Grauer, R., 2013. Effect of turbulent fluctuations on the drag and lift forces on a towed sphere and its boundary layer. Journal of Fluid Mechanics 721, 155–179.
  * Huang and Tian (2019) Huang, W.-X., Tian, F.-B., 2019. Recent trends and progress in the immersed boundary method. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 233 (23-24), 7617–7636.
  * Izard et al. (2014) Izard, E., Bonometti, T., Lacaze, L., 2014. Simulation of an avalanche in a fluid with a soft-sphere/immersed boundary method including a lubrication force. Journal of Computational Multiphase Flows 6 (4), 391–406.
  * Jain et al. (2019) Jain, R., Tschisgale, S., Fröhlich, J., 2019. A collision model for dns with ellipsoidal particles in viscous fluid. International Journal of Multiphase Flow 120, 103087.
  * Jain et al. (2021) Jain, R., Tschisgale, S., Fröhlich, J., 2021. Impact of shape: DNS of sediment transport with non-spherical particles. Journal of Fluid Mechanics 916, A38.
  * Ji et al. (2014) Ji, C., Munjiza, A., Avital, E., Xu, D., Williams, J., 2014. Saltation of particles in turbulent channel flow. Physical Review E 89 (5), 052202.
  * Joseph and Hunt (2004) Joseph, G., Hunt, M., 2004. Oblique particle-wall collisions in a liquid. Journal of Fluid Mechanics 510, 71–93.
  * Joseph et al. (2001) Joseph, G., Zenit, R., Hunt, M., Rosenwinkel, A., 2001. Particle-wall collisions in a viscous fluid. Journal of Fluid Mechanics 433, 329–346.
  * Kempe and Fröhlich (2012a) Kempe, T., Fröhlich, J., 2012a. Collision modelling for the interface-resolved simulation of spherical particles in viscous fluids. Journal of Fluid Mechanics 709, 445–489.
  * Kempe and Fröhlich (2012b) Kempe, T., Fröhlich, J., 2012b. An improved immersed boundary method with direct forcing for the simulation of particle laden flows. Journal of Computational Physics 231 (9), 3663–3684.
  * Kidanemariam et al. (2022) Kidanemariam, A. G., Scherer, M., Uhlmann, M., 2022. Open-channel flow over evolving subaqueous ripples. Journal of Fluid Mechanics 937, A26.
  * Kidanemariam and Uhlmann (2014a) Kidanemariam, A. G., Uhlmann, M., 2014a. Direct numerical simulation of pattern formation in subaqueous sediment. Journal of Fluid Mechanics 750.
  * Kidanemariam and Uhlmann (2014b) Kidanemariam, A. G., Uhlmann, M., 2014b. Interface-resolved direct numerical simulation of the erosion of a sediment bed sheared by laminar channel flow. International Journal of Multiphase Flow 67, 174–188.
  * Kim and Moin (1985) Kim, J., Moin, P., 1985. Application of a fractional-step method to incompressible navier-stokes equations. Journal of Computational Physics 59 (2), 308–323.
  * Li et al. (2020) Li, Q., Abbas, M., Morris, J. F., 2020. Particle approach to a stagnation point at a wall: Viscous damping and collision dynamics. Physical Review Fluids 5 (10), 104301.
  * Luo et al. (2017) Luo, K., Hu, C., Wu, F., Fan, J., 2017. Direct numerical simulation of turbulent boundary layer with fully resolved particles at low volume fraction. Physics of Fluids 29 (5), 053301.
  * Luo et al. (2007) Luo, K., Wang, Z., Fan, J., Cen, K., 2007. Full-scale solutions to particle-laden flows: Multidirect forcing and immersed boundary method. Physical Review E 76 (6), 066709.
  * Marshall (2009) Marshall, J. S., 2009. Discrete-element modeling of particulate aerosol flows. Journal of Computational Physics 228 (5), 1541–1561.
  * Meyer-Peter and Müller (1948) Meyer-Peter, E., Müller, R., 1948. Formulas for bed-load transport. In: IAHSR 2nd Meeting, Stockholm, appendix 2. IAHR.
  * Mindlin and Deresiewicz (1953) Mindlin, R. D., Deresiewicz, H., 1953. Elastic spheres in contact under varying oblique forces. Journal of Applied Mechanics 20 (3), 327–344.
  * Mittal and Iaccarino (2005) Mittal, R., Iaccarino, G., 2005. Immersed boundary methods. Annual Review of Fluid Mechanics 37, 239–261.
  * Papista et al. (2011) Papista, E., Dimitrakis, D., Yiantsios, S., 2011. Direct numerical simulation of incipient sediment motion and hydraulic conveying. Industrial & Engineering Chemistry Research 50 (2), 630–638.
  * Peng et al. (2020) Peng, C., Ayala, O. M., Wang, L.-P., 2020. Flow modulation by a few fixed spherical particles in a turbulent channel flow. Journal of Fluid Mechanics 884, A15.
  * Picano et al. (2015) Picano, F., Breugem, W.-P., Brandt, L., 2015. Turbulent channel flow of dense suspensions of neutrally buoyant spheres. Journal of Fluid Mechanics 764, 463–487.
  * Rao and Capecelatro (2019) Rao, A. A., Capecelatro, J., 2019. Coarse-grained modeling of sheared granular beds. International Journal of Multiphase Flow 114, 258–267.
  * Ray et al. (2015) Ray, S., Kempe, T., Fröhlich, J., 2015. Efficient modelling of particle collisions using a non-linear viscoelastic contact force. International Journal of Multiphase Flow 76, 101–110.
  * Rettinger and Rüde (2022) Rettinger, C., Rüde, U., 2022. An efficient four-way coupled lattice Boltzmann–discrete element method for fully resolved simulations of particle-laden flows. Journal of Computational Physics 453, 110942.
  * Roma et al. (1999) Roma, A. M., Peskin, C. S., Berger, M. J., 1999. An adaptive version of the immersed boundary method. Journal of Computational Physics 153 (2), 509–534.
  * Seyed-Ahmadi and Wachs (2020) Seyed-Ahmadi, A., Wachs, A., 2020. Microstructure-informed probability-driven point-particle model for hydrodynamic forces and torques in particle-laden flows. Journal of Fluid Mechanics 900, A21.
  * Shao et al. (2012) Shao, X., Wu, T., Yu, Z., 2012. Fully resolved numerical simulation of particle-laden turbulent flow in a horizontal channel at a low Reynolds number. Journal of Fluid Mechanics 693, 319–344.
  * Tao et al. (2018) Tao, S., Zhang, H., Guo, Z., Wang, L.-P., 2018. A combined immersed boundary and discrete unified gas kinetic scheme for particle–fluid flows. Journal of Computational Physics 375, 498–518.
  * Tenneti and Subramaniam (2014) Tenneti, S., Subramaniam, S., 2014. Particle-resolved direct numerical simulation for gas-solid flow model development. Annual Review of Fluid Mechanics 46, 199–230.
  * Tschisgale et al. (2018) Tschisgale, S., Kempe, T., Fröhlich, J., 2018. A general implicit direct forcing immersed boundary method for rigid particles. Computers & Fluids 170, 285–298.
  * Uhlmann (2005) Uhlmann, M., 2005. An immersed boundary method with direct forcing for the simulation of particulate flows. Journal of Computational Physics 209 (2), 448–476.
  * Uhlmann (2008) Uhlmann, M., 2008. Interface-resolved direct numerical simulation of vertical particulate channel flow in the turbulent regime. Physics of Fluids 20 (5), 053305\.
  * Verzicco (2023) Verzicco, R., 2023. Immersed boundary methods: Historical perspective and future outlook. Annual Review of Fluid Mechanics 55, 129–155.
  * Vowinckel et al. (2016) Vowinckel, B., Jain, R., Kempe, T., Fröhlich, J., 2016. Entrainment of single particles in a turbulent open-channel flow: A numerical study. Journal of Hydraulic Research 54 (2), 158–171.
  * Vowinckel et al. (2014) Vowinckel, B., Kempe, T., Fröhlich, J., 2014. Fluid–particle interaction in turbulent open channel flow with fully-resolved mobile beds. Advances in Water Resources 72, 32–44.
  * Vreman (2016) Vreman, A., 2016. Particle-resolved direct numerical simulation of homogeneous isotropic turbulence modified by small fixed spheres. Journal of Fluid Mechanics 796, 40–85.
  * Wang et al. (2019) Wang, S., Vanella, M., Balaras, E., 2019. A hydrodynamic stress model for simulating turbulence/particle interactions with immersed boundary methods. Journal of Computational Physics 382, 240–263.
  * Wang et al. (2022) Wang, Y., Zhu, Z., Hu, R., Shen, L., 2022. Direct numerical simulation of a stationary spherical particle in fluctuating inflows. AIP Advances 12 (2), 025019\.
  * Wong and Parker (2006) Wong, M., Parker, G., 2006. Reanalysis and correction of bed-load relation of Meyer-Peter and Müller using their own database. Journal of Hydraulic Engineering 132 (11), 1159–1168.
  * Xia et al. (2020) Xia, Y., Xiong, H., Yu, Z., Zhu, C., 2020. Effects of the collision model in interface-resolved simulations of particle-laden turbulent channel flows. Physics of Fluids 32 (10), 103303.
  * Xia et al. (2022a) Xia, Y., Yu, Z., Lin, Z., Guo, Y., 2022a. Improved modelling of interfacial terms in the second-moment closure for particle-laden flows based on interface-resolved simulation data. Journal of Fluid Mechanics 952, A25.
  * Xia et al. (2022b) Xia, Y., Yu, Z., Pan, D., Lin, Z., Guo, Y., 2022b. Drag model from interface-resolved simulations of particle sedimentation in a periodic domain and vertical turbulent channel flows. Journal of Fluid Mechanics 944, A25.
  * Yang and Hunt (2006) Yang, F.-L., Hunt, M., 2006. Dynamics of particle-particle collisions in a viscous liquid. Physics of Fluids 18 (12), 121506.
  * Yang et al. (2018) Yang, Z., Deng, B.-Q., Shen, L., 2018. Direct numerical simulation of wind turbulence over breaking waves. Journal of Fluid Mechanics 850, 120–155.
  * Yang et al. (2017) Yang, Z., Lu, X.-H., Guo, X., Liu, Y., Shen, L., 2017. Numerical simulation of sediment suspension and transport under plunging breaking waves. Computers & Fluids 158, 57–71.
  * Yu and Shao (2007) Yu, Z., Shao, X., 2007. A direct-forcing fictitious domain method for particulate flows. Journal of Computational Physics 227 (1), 292–314.
  * Yu and Shao (2010) Yu, Z., Shao, X., 2010. Direct numerical simulation of particulate flows with a fictitious domain method. International Journal of Multiphase Flow 36 (2), 127–134.
  * Yu et al. (2021) Yu, Z., Xia, Y., Guo, Y., Lin, J., 2021. Modulation of turbulence intensity by heavy finite-size particles in upward channel flow. Journal of Fluid Mechanics 913, A3.
  * Zeng et al. (2022a) Zeng, Y., Bhalla, A. P. S., Shen, L., 2022a. A subcycling/non-subcycling time advancement scheme-based dlm immersed boundary method framework for solving single and multiphase fluid–structure interaction problems on dynamically adaptive grids. Computers & Fluids 238, 105358\.
  * Zeng et al. (2022b) Zeng, Y., Xuan, A., Blaschke, J., Shen, L., 2022b. A parallel cell-centered adaptive level set framework for efficient simulation of two-phase flows with subcycling and non-subcycling. Journal of Computational Physics 448, 110740.
  * Zenit et al. (1997) Zenit, R., Hunt, M. L., Brennen, C. E., 1997. Collisional particle pressure measurements in solid-liquid flows. Journal of Fluid Mechanics 353, 261–283.
  * Zhang et al. (2020) Zhang, Y., Jiang, M., Chen, X., Yu, Y., Zhou, Q., 2020. Modeling of the filtered drag force in gas–solid flows via a deep learning approach. Chemical Engineering Science 225, 115835.
  * Zhou and Fan (2014) Zhou, Q., Fan, L.-S., 2014. A second-order accurate immersed boundary-lattice boltzmann method for particle-laden flows. Journal of Computational Physics 268, 269–301.
  * Zhou and Fan (2015) Zhou, Q., Fan, L.-S., 2015. Direct numerical simulation of low-reynolds-number flow past arrays of rotating spheres. Journal of Fluid Mechanics 765, 396–423.
  * Zhu et al. (2022) Zhu, Z., Hu, R., Lei, Y., Shen, L., Zheng, X., 2022. Particle resolved simulation of sediment transport by a hybrid parallel approach. International Journal of Multiphase Flow 152, 104072.