A collision model for grain-resolving simulations of flows over dense, mobile, polydisperse granular sediment beds

1. Introduction

• The critical Shields parameter poorly predicts particle-erosion onset, and experimental difficulty measuring dense particle-laden flows has slowed progress.Numerical simulations, particularly DNS with the Immersed Boundary Method, offer an alternative for studying fluid-particle interactions in these conditions.Source: §1 · p. 2
• Particle-particle contact is ubiquitous in dense sediment beds, yet hard-sphere models cannot represent simultaneous collisions or enduring contacts and instead leave particles in constant minute motion.These limitations motivate collision models able to represent enduring contact among multiple interacting particles.Source: §1 · p. 3
• The repulsive-potential model creates an artificial gap between colliding particles that protrudes them into the flow, a critical issue because particle protrusion affects mobilization.In IBM, the gap is introduced so fluid between particle surfaces remains resolvable.Source: §1 · pp. 3–4
• Choosing repulsive stiffness is problematic: values high enough for impacts can cause excessive rebound, while stiffness too low lengthens collisions and can dampen bed motion unphysically.The introduction notes sediment transport spans enduring contacts as well as high-impact collisions, making a single stable stiffness difficult to select.Source: §1 · p. 4
• Prior sediment-flow simulations using Glowinski et al. (1999) were limited to small volume fractions, while Kidanemariam and Uhlmann (2014) avoided stiffness calibration but retained an artificial gap.The former did not permit conclusions about particle-particle interaction.Source: §1 · p. 4
• A more consistent approach combines lubrication forces for near-contact fluid effects with normal and tangential contact forces for material deformation and friction (Simeonov and Calantoni, 2012; Kempe and Fröhlich, 2012a; Izard et al., 2014; Costa et al., 2015; Sierakowski and Prosperetti, 2016).Lubrication accounts for fluid forces that the computational mesh cannot resolve.Source: §1 · p. 5
• The model extends Kempe and Fröhlich (2012a) to dense packing, addressing neglected hydrodynamic forces during contact and a tangential force model that does not converge for enduring contact.The latter model was designed to enforce zero slip between particles.Source: §1 · pp. 5–6
• Prior work had not addressed enduring contact or large-scale collective validation, while polydisperse sediment studies lacked particle-flow feedback and benchmark validation for settling or wall collisions.The cited Fukuoka et al. (2014) study used front tracking for polydisperse sediment but did not account for particle feedback on the flow.Source: §1 · p. 6
• The work targets collision models for polydisperse sediment, avoiding artificial gaps and adaptively calibrating stiffness across Stokes numbers.These are among the paper’s stated key challenges.Source: §1 · p. 7
• It also develops criteria for enduring normal and oblique contact while minimizing tunable parameters.For enduring contact, the authors retain hydrodynamic forces, buoyant weight, and collision forces in the particle momentum balance.Source: §1 · p. 7
• The implementation uses Kempe and Fröhlich (2012a)’s adaptive procedure for normal forces and Thornton et al. (2013)’s DEM-derived tangential model.Both approaches are extended to enduring contact.Source: §1 · p. 7
• The enhancements reproduce laboratory benchmark cases for binary collisions, and the study validates fluid and particle velocity profiles and bulk flow quantities against Aussillous et al.’s data.The supplied passage describes the validation scope but does not state its quantitative outcome.Source: §1 · p. 7

2.1. Fluid solver

• The incompressible Newtonian fluid is solved on a uniform grid, with IBM forcing enforcing no-slip at particle surfaces through interpolation and spreading between Eulerian and Lagrangian points.Source: §2.1 · p. 8

2.2. Computation of particle motion

• The IBM fluid-particle coupling is checked against sphere-settling experiments; the unbounded case reaches terminal velocity, while the wall-bounded case decelerates as it approaches the wall.The unbounded and wall-bounded cases use Rep = 41 and Rep = 12, respectively.Source: §2.2 · pp. 9–10, Figure 1 · p. 10

2.3. Structure of the collision model

• Near a wall, overlapping IBM marker supports are disabled to avoid out-of-domain information and competition in enforcing no-slip.Source: §2.3 · pp. 10–11
• Lubrication models subgrid and omitted-marker fluid forces within 2h; contact begins at ζn ≤0 to limit overlap and represent momentum transfer and energy loss.Source: §2.3 · p. 11
• For polydisperse beds, each particle’s collision force sums normal and tangential interactions with other particles and the wall, allowing arbitrary radii for colliding particles.Source: §2.3 · pp. 12–13
• Appendix A, equation (A.6), defines a quantity that accounts for surface overlap.Source: §2.3 · p. 13
• The text says mathematical descriptions of the models used will follow in the next section.Source: §2.3 · p. 13

3.1. Lubrication model

• Lubrication forces oppose relative motion during approach and rebound, and remain constant below ζn,min to provide continuous forcing near sustained contact.The gap process is unresolved below ζn < 2h; ζn,min can be interpreted as particle-surface micro-texture.Source: §3.1 · pp. 13–14

3.2. Normal contact model

• The explicit ACTM formulation's uout error stays below 1.3% for edry > 0.7 and 3% for edry > 0.4.ACTM adaptively sets stiffness and damping from impact velocity, desired restitution, and collision time; stretching collisions accommodates the fluid solver's longer timestep.Source: §3.2 · pp. 14–15

3.3. Tangential collision model

• The tangential spring-dashpot model uses relative tangential velocity and accumulated spring displacement to represent frictional contact.The accumulated displacement provides frictional memory, permitting a steady-state frictional bed configuration.Source: §3.3 · p. 16

• The tangential force is limited according to Coulomb's friction criterion.Source: Appendix C · p. 16
• Stored tangential displacement provides frictional memory that permits steady-state bed configurations, unlike slip-only friction response.The contrast is with models using only gt,cp, including the model proposed by Kempe and Fröhlich (2012a).Source: Appendix C · p. 16
• The selected kt and dt values yield correct rebound characteristics for oblique impacts after creating a uniform collision time Tc.The model's parameters can be set from experimentally obtained material properties rather than calibration.Source: Appendix C · p. 17
• The model requires no calibration, using experimentally obtained material properties to parameterize its frictional contact.Source: Appendix C · p. 17

4.1. Motivation

• Higher-accuracy time integration and more particle-motion timesteps, without changing the fluid timestep, enabled robust reproduction of Gondret et al. (2002) collision trajectories.Source: §4.1 · p. 17

4.2. Improved time integration

• The collision update uses a predictor-corrector with third-order velocity prediction and second-order Crank-Nicolson remaining steps, improving accuracy at minimal computational cost.It reuses the three-step Runge-Kutta scheme used for the Navier-Stokes equations.Source: §4.2 · pp. 18–19

4.3. Temporal substepping

• For St = 27, changing only the initial particle position within one grid cell produced rebound-height variation up to 83%, showing sensitivity to initial conditions.Five simulations varied the initial position y0 within an interval of one grid cell h.Source: §4.3 · p. 20, Figure 5 · p. 22
• Lubrication remains under-resolved near ζn = 0 before and after contact, so timestep placement changes total impulse and simulation variability, with the effect depending strongly on St.The normal contact force evolves smoothly at the fluid solver timestep, but lubrication does not.Source: §4.3 · p. 22
• The method uses 15 equal substeps per fluid timestep (∆tsub = ∆t/15), solving particle motion with three-step RK while re-evaluating lubrication and contact forces and holding hydrodynamic forces constant.Source: §4.3 · pp. 22–23
• Using a collision timestep 15 times smaller than the fluid timestep, the method almost eliminates variability in particle-wall rebound trajectories; the trajectory resolved with fifteen times more timesteps better approximates the implemented model’s exact solution.Source: §4.3 · pp. 23–24

4.5. Particle momentum balance for high Stokes number collisions

• Including gravitational force during contact negligibly changes the desired Tc and edry, whereas neglecting fluid forces can be unphysical for enduring contact, defined by contact timescales matching or exceeding the flow timescale.The authors include buoyant weight during contact; enduring contact requires a different treatment from brief high-Stokes-number collisions.Source: §4.5 · pp. 26–27
• Excluding fluid forces during high-Stokes-number contact agrees with experiments, while including them causes excessive damping because the surrounding fluid cannot adapt to the particle’s changed motion.Neglecting fluid forces has minimal effect on the general flow when contact is much shorter than the flow timescale.Source: §4.5 · pp. 26–27
• Under Eq. (21), a particle in contact with the wall does not experience hydrodynamic forces.Source: §4.5 · p. 27
• It can therefore remain on the wall until another particle collides with it, rather than being swept into the flow or carried downstream.Source: §4.5 · p. 27
• The authors address this problem in Section 5 by introducing a threshold for including hydrodynamic forces in Eqs. (3) and (4).Source: §4.5 · p. 27

5.1. Accounting for fluid forces

• The model includes fluid forces below Stcrit and retains hydrodynamic stresses and buoyant weight in enduring-contact momentum balance, making Shields-parameter considerations applicable.Source: §5.1 · pp. 27, 29, 30
• At higher Stokes numbers, including hydrodynamic forces during contact significantly undershoots the rebound trajectory.The undershooting decreases as the Stokes number decreases.Source: §5.1 · p. 28
• For St < 5, collisions show no appreciable rebound, so the model sets Stcrit = 5 for enduring contact.The comparison uses particle velocity at contact.Source: §5.1 · p. 28, Figure 9 · p. 29

5.2. Optimizing enduring particle overlap

• To limit overlap as low impact velocities drive computed stiffness upward, stiffness uses critical-velocity and gravity-overlap constraints with ϵ = 10^-3.The prescription aims for a uniform set of collision stiffnesses that minimizes particle overlap.Source: §5.2 · pp. 30–31

5.3. Rolling and sliding motion

• The friction scheme tests sticking with static friction and switches to kinetic friction once slipping occurs, limiting competing no-slip contacts that can destabilize calculations.Source: §5.3 · pp. 31–32
• The simulations use kinetic friction coefficient µk = 0.15 and static coefficient µs = 0.8, based on prior measurements for glass spheres and rock materials.The cited sources are Joseph and Hunt (2004) for µk and Dieterich (1972) for µs.Source: §5.3 · p. 32
• The model reproduces Foerster et al. (1994)'s experimental rebound-angle dependence on impact angle in dry oblique particle-wall collisions.Figure 10 shows rebound versus impact angle.Source: §5.3 · pp. 32–33, Figure 10 · p. 33
• At low impact angles, the linear-spring tangential model captures negative rebound angles, which occur when release leaves |up| < |ωp,z|.Perfect sticking at release instead gives a zero rebound angle.Source: §5.3 · p. 33
• In Couette flow, ReH=50 produces perfect rolling without slip, whereas ReH=10 produces slipping.Figure 11 shows the translational and rotational velocities for the two cases.Source: Figure 11 · p. 34, §5.3 · p. 34
• Fluid forces let a wall-contacting particle reach steady speed; at lower Reynolds number, increased viscous drag overcomes friction from its weight, causing slip.The Couette-flow simulation accounts for fluid forces while the particle is in contact with the wall.Source: §5.3 · p. 34

6.1. Physical setup

• With periodic boundary conditions, the simulations aimed only to reproduce long-term steady-state flow conditions; obtained fluid heights hf were larger than experimental values in every reported case except A10M.Source: §6.1 · pp. 36, 39
• A low CFL of 0.1 was required for coupling stability amid many simultaneous particle interactions; coarser grid resolution showed no appreciable bulk-flow difference.The simulations used Dp/h = 22.7 to resolve interstitial flow.Source: §6.1 · p. 38

6.2. Comparison of wall-normal profiles

• Particle-velocity profiles compare very well with experiments for cases A9 and A10, and reasonably well for cases A1 and A2.Figure 13 shows the wall-normal profile comparisons.Source: Figure 13 · p. 39, §6.2 · p. 40
• Fluid velocity is parabolic above the bed; at the interface it slips relative to particles, while the two profiles are similar within the bed.The slight in-bed difference is attributed to flow between particles.Source: §6.2 · p. 40
• Profile discrepancies partly reflect bed-height and Shields-number differences, which are largest for A1 and A2 and may explain their larger velocity-profile deviations.The simulations’ fluid heights exceed the experimental values except in case A10M.Source: §6.2 · p. 40
• The wall-normal-profile comparison notes that results can be sensitive to fluid height hf.Source: §6.2 · p. 40

6.3. Comparison of bulk quantities

• Simulated particle velocity fluxes agree well with experiments but generally underestimate mean values; cases A9 and A10 match better, and results remain within experimental error.Table 5 compares simulated and experimental velocity fluxes.Source: §6.3 · p. 41, Table 5 · p. 41

6.4. Polydisperse flow

• A10M shows no appreciable change in bed volume fractions versus A10, and its average bed volume fraction is φ = 0.609, consistent with random sphere packing.Figure 14 presents wall-normal profiles of average particle volume fractions.Source: §6.4 · p. 42, Figure 14 · p. 42
• A10M has a slightly higher velocity profile than A10, likely because its lower fluid height hf produces a higher Shields number and an overpredicted mean experimental flux.The flux comparison is reported in Table 5.Source: Table 5 · p. 41, §6.4 · p. 42
• Despite the flux overprediction, the polydisperse particle velocity profile and flux agree very well with experiments, supporting monodisperse spheres as a valid approximation for this setup.The authors attribute the slightly increased velocity profile to the lower fluid height in A10M.Source: §6.4 · p. 42

7. Conclusions

• The model is intended for phase-resolved simulations in which the dispersed phase is represented by the Immersed Boundary Method.The authors state that the achieved accuracy will enable detailed analysis of phase-resolved simulation data.Source: §7 · pp. 42, 44
• The model accounts for lubrication, normal repulsion, and tangential friction without parameters requiring arbitrary calibration.It also supports rolling and sliding contact, distinguishing sticking from sliding.Source: §7 · pp. 42–43
• An improved time-integration scheme is necessary for consistent particle-wall collision results.The model also extends to flows over dense granular sediments.Source: §7 · p. 43
• The enhancements are crucial for thick sediment packings, which can reach several particle diameters and remain numerically stable at steady state.The reported simulations retain the full momentum balance of enduringly contacting particles, including hydrodynamic forces and buoyant weight.Source: §7 · p. 43
• The model achieves satisfactory agreement in binary-collision benchmarks and collective motion over a dense granular packing, and presents a first polydisperse-sediment test.The collective case is a horizontal flow over a dense packing.Source: §7 · p. 43
• The authors suggest, without explicitly establishing, possible applicability to vertical channel flows and neutrally buoyant particles in laminar and turbulent conditions.They describe the approach as a tool for high-fidelity, grain-scale data in multiphase flows involving rigid spheres.Source: §7 · p. 44

Appendix A. Definitions for particle-particle and particle-wall collisions

• P-P denotes particle-particle interaction, P-W particle-wall interaction, and P-F particle-fixed interaction; fixed-particle collisions are generally handled like P-P collisions unless stated otherwise.The definitions distinguish wall and particle interactions where needed.Source: Appendix A · p. 51
• The contact notation includes the closest wall point, contact normal, surface gap, contact location, and particle radius relative to the contact point.Effective radius and effective mass are also defined.Source: Appendix A · pp. 51–52
• Relative center-of-mass velocity and contact-point velocity are distinguished, with the latter resolved into normal and tangential components.The appendix gives a relation for contact-point relative velocity in particle-particle interactions.Source: Appendix A · p. 52

Appendix B. Calculating the normal contact model coefficients

• Ray et al. (2015) derive algebraic expressions for the normal-contact stiffness and damping coefficients using nonlinear transformations and a series expansion of equation (11).Source: Appendix B · p. 52
• The expressions use A = 0.716, B = 0.830, C = 0.744, α = 1.111, and τc,0 = 3.218.Source: Appendix B · p. 53
• The impact velocity is measured as uin = gn,cp · n at the first occurrence.Source: Appendix B · p. 53

Appendix C. The tangential displacement vector

• The tangential spring model uses displacement representing accumulated relative motion between two surfaces.Source: Appendix C · p. 53
• The discrete update rotates the previous timestep’s displacement onto the plane tangent to the surfaces, accounting for a changing reference frame.Without this rotation, the linear spring could contribute to the normal force between particles.Source: Appendix C · p. 53
• During Coulomb slip, the displacement is reset to the value that achieves the Coulomb friction force, preventing unbounded growth as surfaces continue sliding.Source: Appendix C · p. 53