A collision model for grain-resolving simulations of flows over dense, mobile, polydisperse granular sediment beds
AA collision model for grain-resolving simulations of flowsover dense, mobile, polydisperse granular sediment beds
Edward Biegert , Bernhard Vowinckel , Eckart Meiburg a Department of Mechanical Engineering, University of California, Santa Barbara, CA, USA
Abstract
We present a collision model for phase-resolved Direct Numerical Simulations ofsediment transport that couple the fluid and particles by the Immersed Bound-ary Method. Typically, a contact model for these types of simulations comprisesa lubrication force for particles in close proximity to another solid object, a nor-mal contact force to prevent particles from overlapping, and a tangential contactforce to account for friction. Our model extends the work of previous authorsto improve upon the time integration scheme to obtain consistent results forparticle-wall collisions. Furthermore, we account for polydisperse spherical par-ticles and introduce new criteria to account for enduring contact, which occursin many sediment transport situations. This is done without using arbitrary val-ues for physically-defined parameters and by maintaining the full momentumbalance of a particle in enduring contact. We validate our model against severaltest cases for binary particle-wall collisions as well as the collective motion ofa sediment bed sheared by a viscous flow, yielding satisfactory agreement withexperimental data by various authors.
Keywords:
Direct Numerical Simulations, Immersed Boundary Method,contact modeling, particle-laden flow, sediment transport [email protected] [email protected] [email protected] Preprint submitted to Journal of Computational Physics September 30, 2018 c (cid:13) http: // creativecommons. org/ licenses/ by-nc-nd/ 4. 0/ a r X i v : . [ phy s i c s . f l u - dyn ] A p r . Introduction The flow over dense, mobile granular beds plays a central role in multiple ap-plications in environmental, mechanical, and process engineering. Prime exam-ples of this type of problem are turbidity currents and powder snow avalanches(Meiburg and Kneller, 2010), for which resuspension of particles essentially de-termines the dynamics of the flow. The resuspension threshold is quantified bythe ratio of hydrodynamic drag and lift forces to the weight of the particles,known as the Shields number (Shields, 1936). The critical Shields parameter,however, has proven to be a poor predictor for the onset of particle erosion(Garcia, 2008), and substantial efforts have been made in hydraulic engineeringto overcome this difficulty (Seminara, 2010). To date, progress has been slowdue to the experimental difficulty of measuring dense particle-laden flows in alaboratory. In recent years, numerical simulations have provided an alterna-tive way to study fluid-particle interactions under these complex conditions. Inparticular, the Direct Numerical Simulation (DNS) of particle-laden flows usingthe Immersed Boundary Method (IBM) has become a very attractive option.A main advantage of this approach is that it allows for a fully-coupled systemby accounting for particle-fluid, fluid-particle, and particle-particle interactions.While coupling between the particle and the fluid can be realized by a straight-forward implementation of an IBM, particle-particle coupling has to be modeledby suitable expressions for the collision processes involved.A number of different phase-resolving methods to simulate mobile particlesin a viscous flow have been developed in the past two decades. Glowinski et al.(1999) and Patankar et al. (2000) developed a Distributed Lagrange Multiplier(DLM)/fictitious domain approach that forces the fluid throughout the volumeof the particle to move in a rigid body motion with the particle. While thismethod was designed for a finite element framework, Kajishima et al. (2001)and Apte et al. (2009) later developed different versions to be used in a finitedifference framework. More recently, other methods that enforce the no-slipcondition on the particle surface have been developed. Zhang and Prosperetti22005) developed PHYSALIS, which uses the analytical solution for Stokes flowaround a sphere. Uhlmann (2005) developed an IBM that enforces the no-slipcondition using interpolation and spreading operations via Dirac delta functions.Kempe and Fr¨ohlich (2012b) extended this method to make it stable for a largerrange of particle/fluid density ratios. Although there exists a variety of phase-resolving methods, the choice of collision models, on the other hand, has notbeen as diverse. Derksen (2011) and Derksen (2015) used a hard-sphere model,which resolves collisions instantaneously. Glowinski et al. (1999) have developeda repulsive potential (RP) model that prevents particles from overlapping byapplying a repulsive force at some small distance before the particles come incontact. Many other authors have adopted this model for simulations involvingdilute suspensions of particles. For example, Uhlmann (2008) and Santarelliand Fr¨ohlich (2015) investigated particles in a vertical turbulent channel flow,Lucci et al. (2010) studied the impact of finite size particles on isotropic turbu-lence, and Breugem (2012) and Picano et al. (2015) have presented results fora horizontal flow laden with neutrally-buoyant particles. In these simulations,particles rarely came in contact, and thus were successfully governed mostly bythe IBM, using the repulsive potential only to prevent overlap.For the situation involving shear flow over a densely-packed sediment bed,however, particle-particle contact becomes ubiquitous (Balachandar and Eaton,2010). Hard sphere models cannot maintain simultaneous collisions or enduringcontacts between multiple interacting particles, but instead represent sedimentbeds as particles in constant, minute motion (Derksen, 2015). However, theyhave been used to reproduce critical erosion conditions for a laminar shear flow(Derksen, 2011). The drawbacks of the RP model for the situation of sedimenttransport have been clearly-elaborated by Kempe et al. (2014). Using the RPin the framework of the IBM introduces an artificial gap of two times the gridcell size between colliding particles so that the fluid in the gap between theparticle surfaces can still be resolved (Glowinski et al., 2001). For the situationof sediment transport, however, the artificial gap also introduces an unphysi-cal protrusion of the particles into the horizontal flow, which is critical as the3rotrusion has been acknowledged to be a very sensitive parameter for particlemobilization (Fenton and Abbott, 1977). In addition, the RP model introducesa material stiffness k n which has to be calculated a priori to design a collisionmodel that is numerically stable. If the value of k n is chosen too high, the repul-sive force is overestimated and the particle would experience an unphysical highrebound velocity. On the other hand, if k n is too low, the duration of the particlecollision would be too large. The resuspension mechanisms, however, generatea high variety of particle impact velocities u in for the collisions, typically char-acterized by the nondimensional Stokes number St = ρ p u in D p / (9 ρ f ν f ), where ρ p and ρ f are the particle and fluid density, respectively, D p is the particlediameter, and ν f is the kinematic viscosity of the fluid. In fact, for bed-loadtransport in water the variety can span from St (cid:29)
10 (saltating particles, Ni˜noand Garc´ıa, 1998) to St = O (10) (rolling particles, Lajeunesse et al., 2010) toSt (cid:28) k n , but still relies on anartificial gap distance. Thus, in order to obtain appropriate bulk sediment trans-port quantities, they calibrated the dry restitution coefficient e dry = − u out /u in ,where u out is the rebound velocity as soon as the collision process is finished,although this parameter can be set exactly as a material property. It describesthe dissipation of kinetic energy due to the inelastic mechanics of the dry con-tact and is typically in the range of 0 . ≤ e dry < e dry = 0 . < ζ n ≤ h , where ζ n is thedistance between the two surfaces), and a contact force when the surfaces comein contact and slightly overlap ( ζ n ≤ ψ in . In this reference, however, they neitherconsidered the situation of enduring contact nor conducted a validation on alarger scale addressing the collective effects of particle motion.Finally, another aspect that has received far less attention so far is the sit-uation of a sediment consisting of particles with varying particle diameters.Interestingly, all of the references cited so far deal with spherical monodisperseparticles. To the knowledge of the authors the only study considering horizontalchannel flows laden with polydisperse sediment has been performed by Fukuokaet al. (2014) using a front-tracking technique, but neither did the authors ac-count for the feedback of the particles on the flow nor have they provided avalidation for the experimental standard benchmark test cases such as particlessettling in an ambient fluid or colliding with a wall. The absence of studies ad-dressing polydisperse sediment with fully coupled IBM simulations is ever themore surprising, since its impact has been acknowledged as a key issue in thedevelopment and evolution of bedforms by segregation effects as reviewed byCharru et al. (2013). 6s a consequence, the present work aims to resolve the problems mentionedabove. Among the key challenges identified are i) deriving collision models forpolydisperse sediment, ii) avoiding the introduction of an artificial gap betweencolliding particles, iii) adaptively-calibrating the particle stiffness to simulate awide range of Stokes numbers in a consistent manner, iv) introducing suitablecriteria to extend existing models towards the numerically-challenging situationof enduring contact for both normal and oblique collisions, and v) minimizingthe number of tunable parameters within the model framework. We achieveour goals by presenting an implementation of collision models for polydispersesediment. We use the adaptive procedure proposed by Kempe and Fr¨ohlich(2012a) for normal forces and the tangential model of Thornton et al. (2013),which stems from DEM. Furthermore, we extend both of these approaches forthe situation of enduring contact. In particular, for enduring contact, we tookcare to retain all the governing terms of the momentum balance of a particle,i.e. hydrodynamic forces, buoyant weight, and collision forces. This measureturns out to be crucial when simulating flows over sediment beds, as the Shieldsparameter is based on the ratio of hydrodynamic to buoyant forces. The pro-posed enhancements allow us to reproduce several laboratory benchmark testcases for binary collisions. In addition, we present a detailed validation of oursimulation results with wall-normal profiles of the fluid and particle velocitiesas well as bulk flow quantities using the experimental data of Aussillous et al.(2013).The paper is structured as follows. We briefly recall the numerical method,including the fluid solver, IBM, and the structure of the collision model, inSection 2, followed by the mathematical description of the collision model em-ployed in Section 3. We then present necessary enhancements to the collisionmodel to deal with small Stokes numbers (Section 4) and to simulate densegranular packings with the gross of the particles in enduring contact (Section5). Subsequently, the enhanced model is validated for the collective motion ofpolydisperse sediment sheared by a laminar flow in Section 6.7 . Particle motion and four-way coupling in the framework of theImmersed Boundary Method For the present simulations, we solve the unsteady Navier-Stokes equationsfor an incompressible Newtonian fluid, given by ∂ u ∂t + ∇ · ( uu ) = − ρ f ∇ p + ν f ∇ u + f IBM , (1)and the continuity equation, given by ∇ · u = 0 , (2)on a uniform rectangular grid with grid cell size ∆ x = ∆ y = ∆ z = h . Here, u = ( u, v, w ) T designates the fluid velocity vector in Cartesian components, p the pressure, ν f the kinematic viscosity, t the time, and f IBM an artificialvolume force introduced by the IBM (Mittal and Iaccarino, 2005). This volumeforce, which acts on the right-hand side of (1) in the vicinity of the inter-phase boundaries, connects the motion of the particles to the fluid phase. Thetransfer of quantities, such as force and velocity, between Eulerian points andLagrangian points, i.e. between fluid points of the regular background grid andpoints on the surface of the particle, is performed by interpolation and spreadingoperations via a weighted sum of regularized Dirac delta functions (Uhlmann,2005), of which we use the 3-point stencil function of Roma et al. (1999). Thesource term f IBM is computed in such a way that the no-slip condition at theparticle surface is satisfied. Time advancement is achieved by a fractional stepmethod, a third-order explicit low-storage three-step Runge-Kutta (RK) schemeis employed for the convective terms, and a second-order semi-implicit Crank-Nicolson scheme is used for the viscous terms (Fadlun et al., 2000), which aresolved with the conjugate-gradient method. Spatial derivatives are evaluatedusing second-order central-differencing. The pressure is treated with a directsolver based on Fast Fourier Transformations (FFT). Our code can handle a8ariety of boundary conditions at the different walls, including no-slip, slip,periodic, and inflow/outflow.
Within the framework of the IBM, we calculate the motion of each individualspherical particle by solving an ordinary differential equation for its translationalvelocity u p = ( u p , v p , w p ) T m p d u p d t = (cid:73) Γ p τ · n d A (cid:124) (cid:123)(cid:122) (cid:125) = F h,p + V p ( ρ p − ρ f ) g (cid:124) (cid:123)(cid:122) (cid:125) = F g,p + F c,p , (3)and its angular velocity ω p = ( ω p,x , ω p,y , ω p,z ) T I p d ω p d t = (cid:73) Γ p r × ( τ · n ) d A (cid:124) (cid:123)(cid:122) (cid:125) = T h,p + T c,p . (4)Here, m p is the particle mass, Γ p the fluid-particle interface, τ the hydrodynamicstress tensor, ρ p the particle density, V p the particle volume, g the gravitationalacceleration, I p = 8 πρ p R p /
15 the moment of inertia, and R p the particle radius.Furthermore, the vector n is the outward-pointing normal on the interface Γ p , r = x − x p is the position vector of the surface point with respect to the center ofmass x p of a particle, and F c,p and T c,p are the force and torque due to particlecollisions, respectively. Furthermore, note the designation of the hydrodynamicforce and torque as F h,p and T h,p , respectively, as well as F g,p the force due togravity, which will be used in the following for brevity.We employ the approach of Kempe and Fr¨ohlich (2012b) for evaluating theIBM forces and solving (3) and (4). We validated the fluid-particle couplingof the method against experimental data of a sphere settling in an unboundedquiescent fluid (Mordant and Pinton, 2000) as well as towards a wall (Ten Cateet al., 2002). The particle was initially at rest for both setups and accelerateddownwards due to gravity. For the unbounded case, the particle reached aconstant terminal velocity u ∞ . The Reynolds number Re p = u ∞ D p /ν f basedon this settling velocity is Re p = 41. For the wall-bounded case, the particle9 t/t ref -1.2-1.0-0.8-0.6-0.4-0.20.0 v / v r e f Mordant & PintonPresent t/t ref -0.25-0.20-0.15-0.10-0.050.00 v / v r e f ten Cate et al.Present a) b) Figure 1:
Comparison of the present numerical method against experimental data with v ref = (cid:112) gD p and t ref = (cid:112) D p /g . a) Settling sphere in an infinite medium at Re p =41 (Mordant and Pinton, 2000) and b) settling sphere in a wall-bounded medium at Re p = 12 (Ten Cate et al., 2002). reached a terminal velocity but then decelerated as it approached the wall. Fornow, we only consider the particle motion before impact, which is governedby the IBM. The collision model governing the impact will be described andvalidated further below in the text. The Reynolds number based on the velocitybefore deceleration is Re p = 12. In both cases, the particle was discretized with20 grid cells per diameter. The respective data are plotted in Figure 1, showingexcellent agreement. As mentioned in the introduction, one of the major advantages of the IBMis the direct computation of long-range interactions between the particles. Onlyshort-range interactions and collisions need to be modeled. For example, con-sider a particle approaching and colliding with a wall, as shown in Figure 2. Asthe particle comes close to the wall, two problems need to be dealt with: first,the smoothed Dirac delta functions used for the IBM overlap with the wall, andsecond, the discrete mesh can no longer resolve the fluid being squeezed outfrom between the two surfaces.We can solve the first problem by disabling Lagrangian marker points whosesupports overlap with the wall (red dashed circle in Figure 2a), as was done10) b)
Figure 2:
Regions where short-range interactions occur. Points on particle surfacerepresent Lagrangian markers. Red markers have been turned off. a) Lubricationregion where red dashed circle illustrates support of Dirac delta function and b) Contactregion. by Kempe and Fr¨ohlich (2012b). This means that the forcing by these selectmarkers on both the fluid and the particles is ignored, preventing the particlefrom using undefined information from outside of the domain and from com-peting with the wall for enforcing the no-slip condition. Figure 2 illustrates thered markers that have been disabled. Similarly, overlapping markers betweentwo particles are disabled. We solve the second problem by adding a lubrica-tion force, which models the subgrid forces on the particle due to the narrowgap and also accounts for some of the fluid forces from the disabled Lagrangianmarkers. We apply this force when the particle-wall distance is less than twogrid cells (0 < ζ n ≤ h ), illustrated by the blue region in Figure 2a. Once theparticle comes into contact with the wall ( ζ n ≤ igure 3: Sketch of polydisperse particles in a mobile granular bed and the resultinginteractions due to collisions. forces F n = ζ n > h lubrication model (9) 0 < ζ n ≤ h normal contact model (10) ζ n ≤ F t = ζ n > ζ n ≤ . (6)The interactions of a single particle with its environment, however, are ingeneral more diverse. Let us consider a particle p embedded in a mobile granularbed of polydisperse, spherical particles (Figure 3). The dynamics of this particleare mainly determined by all the collision forces exerted upon it by particles q, q (cid:54) = p as well as the wall. The total force F c,p acting on a particle p duringthe collision process may be decomposed as F c,p = N p (cid:88) q, q (cid:54) = p ( F n,q + F t,q ) + F n,w + F t,w , (7)where N p is the number of particles simulated, F n,q and F n,w are the normalcollision forces described by (5) with particle q and the wall, respectively, and12 t,q and F t,w are the tangential collision forces described by (6) with particle q and the wall, respectively. In what follows, the mathematical expressions areformulated for both particle-wall and particle-particle collisions, where the radiiof the two colliding particles can be arbitrary. Whenever a distinction betweenparticle-wall and particle-particle has to be made, a reference to Appendix Ais given for brevity, providing all definitions and nomenclature needed to distin-guish between the two different cases.The torque T c,p on a spherical particle p generated by the tangential contactforces is T c,p = N p (cid:88) q, q (cid:54) = p R p,cp n p,q × F t,q + R p,cp n p,w × F t,w (8)where n p,q and n p,w are the unit vectors pointing to the collision partner q orthe wall, respectively, and R p,cp is the particle radius at the contact point asdefined per Appendix A in (A.6), which accounts for surface overlap. In thenext section, we will provide the mathematical description of the models usedin the present study.
3. Collision modeling
When the distance between the surfaces of two approaching particles be-comes small, the fluid is squeezed out of the gap. The fluid grid cannot resolvethis process as soon as ζ n < h , where h is the grid cell size. Hence, we employa lubrication model, which also acts on particles rebounding after the collision,when fluid is drawn into the gap. The lubrication force is dissipative, since itis always directed opposite to the relative velocity. The model is based on theanalytical derivation of Cox and Brenner (1967), who solved for the force underStokes flow conditions F n = − πρ f ν f R eff max ( ζ n , ζ n,min ) g n,cp , (9)where R eff is the effective radius accounting for polydisperse sediment and g n,cp the normal component of the relative particle velocity as defined per Appendix13 in (A.1) and (A.9) respectively. The original model scales as 1 /ζ n whichintroduces a singularity as ζ n →
0. This has been addressed by Simeonovand Calantoni (2012) and Kempe and Fr¨ohlich (2012a), who set F n = 0 for0 < ζ n < ζ n,min , and Izard et al. (2014), who have shifted the denominator of(9) to ζ n + ζ n,min . In the present approach, as in the approach of Costa et al.(2015), lubrication forces are held constant as soon as the gap size becomessmaller than the critical value ζ n,min , which provides a continuous forcing onthe particles when they are close to sustained contact. This parameter can beinterpreted as the micro-texture of the particle surface, which acts as a surfaceroughness, as will be discussed further in Section 4.4 below. To account for normal contact forces, we implemented the Adaptive CollisionTime Model (ACTM) proposed by Kempe and Fr¨ohlich (2012a). The main ideaof the ACTM is to use an adaptive procedure to obtain the desired restitutioncoefficient e dry and to resolve the collision on the timescale of the fluid solver.The ACTM is based on a nonlinear spring-dashpot system F n = − k n | ζ n | / n − d n g n,cp , (10)which involves empirical parameters for the coefficients of stiffness k n and damp-ing d n . Here, n is the normal vector pointing either towards the collision partneror towards the wall as defined per Appendix A in (A.3). The nonlinear term | ζ n | / arises from Hertzian contact theory (Hertz, 1882).Since the timescale of a collision according to Hertzian contact theory isseveral orders of magnitude smaller than the typical temporal discretizationof the fluid solver, the collision event has to be stretched in time to maintainthe efficiency of the numerical procedure. This measure is also needed for thefluid to adapt to the sudden change in the particle trajectory (Kempe andFr¨ohlich, 2012a; Costa et al., 2015). However, the duration of contact T c ismainly determined by the stiffness parameter k n . This becomes of particularimportance for the complex situation of bed-load transport, where a broad range14f impact velocities is encountered, ranging from high-impact collisions at thetop of the bed to enduring contact within the bed.The ACTM fixes this problem by adaptively calibrating the parameters k n and d n depending on the impact velocity u in , the desired restitution coeffi-cient e dry , and the desired collision time T c . The latter is a parameter of themodel and should be minimized to avoid excessive particle overlaps and tempo-ral stretching. Kempe and Fr¨ohlich (2012a) demonstrated that T c = 10∆ t is asuitable choice for the collision time given that all timescales related to fluid andparticle motion are significantly larger than the timescale of particle contact.For glass and hard metals, e dry = 0 .
97 is a typical value (e.g Foerster et al.,1994; Joseph et al., 2001; Gondret et al., 2002). For immersed collisions, therestitution coefficient e wet , measured some small distance away from the wall,becomes a function of the Stokes number (Joseph et al., 2001). The ACTM,however, uses the IBM and lubrication model to account for e wet through ad-ditional dissipative fluid effects.In order to find values for k n and d n , we first neglect all non-contact forcesacting on the particle so that (3) and (10) together give the following nonlinearordinary differential equation: m eff d ζ n d t + d n d ζ n d t + k n ζ / n = 0 , (11)where m eff is the effective mass accounting for polydisperse sediment as definedper Appendix A in (A.2). Note that d ζ n / d t = − g n,cp · n . Together with (11) andinitial and final conditions, the constraints e dry and T c allow for determinationof k n and d n using either an iterative procedure, as was done by Kempe andFr¨ohlich (2012a), or an explicit formulation, as was proposed by Ray et al.(2015). In the present study, we implemented the explicit formulation, whichis provided in Appendix B. According to Ray et al. (2015), the error in u out increases with decreasing e dry , but does not exceed 1.3% for e dry > . e dry > .
4, making this method useful for most sediment materials such assilicate, glass, or even metal. 15 .3. Tangential collision model
To account for frictional contact between the particles, we implemented atangential contact model based on the linear spring-dashpot model described inthe review paper of Thornton et al. (2013): F t,LS = − k t ζ t − d t g t,cp , (12)which has stiffness and damping coefficients k t and d t . This model uses g t,cp , thetangential component of the relative surface velocities as described in (A.10) ofAppendix A, as well as ζ t , the tangential spring displacement, which representsthe accumulated relative tangential motion between the two surfaces: ζ t = (cid:90) tt i g t,cp ( t (cid:48) ) dt (cid:48) , (13)where t i is the time of impact. The discretized form of (13) is described inAppendix C.This model limits the maximum force based on Coulomb’s friction criterion: F t = min ( || F t,LS || , || µ F n || ) t , (14)where µ represents the coefficient of friction between the two surfaces (describedfurther in Section 5.3) and t = F t,LS / || F t,LS || points in the direction of thetangential force.This model has two important features for simulating densely-packed beds.First, the spring allows many particles to interact in a smooth, stable manner,provided the stiffness is chosen properly. Second, the model has a memory of thefriction force via the tangential displacement ζ t , which permits a steady-statefrictional bed configuration. In contrast, a model that only uses g t,cp , such asthe one proposed by Kempe and Fr¨ohlich (2012a), can only react to slip, notpredict it.Similarly to the ACTM, we can adaptively compute k t and d t for each col-lision. According to Thornton et al. (2011), the stiffness can be set to k t = κ m eff π T c . (15)16ere, κ is based on Poisson’s ratio ν : κ = 2(1 − ν )2 − ν , (16)which is a well-studied material property typically ranging between 0 . < ν < .
30 (e.g Foerster et al., 1994; Gondret et al., 2002; Joseph and Hunt, 2004).Hence, a value of ν = 0 .
22 was used in the present study.In addition, the damping is computed according to Thornton et al. (2013)to account for the inelasticity of the collisions d t = 2 (cid:112) m eff k t − ln e dry (cid:113) π + ln e dry . (17)Having created a uniform collision time T c with the normal contact model,we obtain the correct rebound characteristics for oblique impacts using thesevalues for k t and d t , as shown in Section 5.3. Consistently, the model doesnot require any calibration but instead can be parameterized using materialproperties obtained from experiments.
4. Enhancements to the normal contact model
In order to obtain a good agreement with immersed collision experiments, wehad to implement a few enhancements to the normal contact model describedin Section 3.2. Both changing the time integration to a scheme of higher accu-racy and adding more timesteps to the integration of particle motion withoutchanging the fluid timestep allowed us to reproduce the collision trajectories ofGondret et al. (2002) in a robust manner.
The ACTM normal contact force F c,p is a function of the surface distance ζ n and the relative velocity g n,cp , which in turn depend on the particle po-sition x k − p and velocity u k − p at the previous substep k −
1. We can write17his functional dependence as F c,p ( x k − p , u k − p ). Integrating the particle equa-tion of motion with a Forward Euler/Crank Nicholson scheme for the particle’svelocity/position, we obtain: u kp = u k − p + 2∆ tα k m p F c,p (cid:0) x k − p , u k − p (cid:1) (18a) x kp = x k − p + ∆ tα k (cid:0) u kp + u k − p (cid:1) , (18b)where k is the number of the RK-substep and α k is the RK-coefficient (Rai andMoin, 1991). For now, we ignore the hydrodynamic, gravitational, and lubrica-tion forces in order to focus on the contact forces alone. We conducted a simpletest to analyze the accuracy of this scheme. A particle of density ρ p /ρ f = 7 . R p = 10 u ∞ ∆ t was initialized with a velocity of u p = (0 , u ∞ , T at a position y p > R p above a wall at y = 0. Subsequently, the particle wasreleased and eventually collided with the wall. Neglecting hydrodynamic ef-fects as well as gravity yields an impact velocity of u in /u ∞ = 1. Choosing thecollision time to be T c = 10∆ t as suggested by Kempe and Fr¨ohlich (2012a),gave good results for the duration of the desired contact phase T c , but ratherlarge errors of the rebound velocity were observed compared to the prescribed e dry = 1. The value of u out = − e dry u in was overestimated by more than 12%.Turning our attention to Figure 4, we can see that the discretization of(18) leads to a poor estimation of the collision force || F c,p || when compared tothe simulation in which 10 timesteps were used to resolve the collision, whichcan be taken as the exact solution. This inaccuracy in the collision force wasobserved for a variety of simulations using different R p , ρ p , u in , and e dry . Inorder to reduce the error to 0 . T c = 1000∆ t would be required, which is not feasible for simulations of sediment transport.Hence, we implemented a temporal discretization scheme with a higher order ofaccuracy. Utilizing the same three-step RK scheme that integrates the Navier-Stokes equations, we reformulated the collision integration with a predictor-18orrector scheme: (cid:101) u p = u k − p + ∆ tm p (cid:2) γ k F c,p (cid:0) x k − p , u k − p (cid:1) + ζ k F c,p (cid:0) x k − p , u k − p (cid:1)(cid:3) (19a) (cid:101) x p = x k − p + ∆ t α k (cid:0)(cid:101) u p + u k − p (cid:1) (19b) u kp = u k − p + ∆ t α k m p (cid:2) F c,p ( (cid:101) x p , (cid:101) u p ) + F c,p (cid:0) x k − p , u k − p (cid:1)(cid:3) (19c) x kp = x k − p + ∆ t α k (cid:0) u kp + u k − p (cid:1) . (19d)Here, tilde indicates predicted values, and γ k and ζ k are the RK coefficientsfor the explicit third-order scheme according to Rai and Moin (1991). Hence,the velocity predictor step (19a) is third-order accurate while the other stepsuse second-order Crank-Nicholson schemes. A similar approach was taken byCosta et al. (2015), but in this reference, the predicted value is determined by aniterative scheme, which is computationally more costly than the present scheme.In Figure 4 we can see that this predictor-corrector scheme yields a much betterapproximation of || F c,p || compared to the Forward Euler method, reducing theerror of u out by almost two orders of magnitude to 0 .
17% for T c / ∆ t = 10. Forcompleteness, we have also included a solution that uses the Backward Eulermethod, which underestimates the rebound velocity by 11%. This improvementhas been achieved by a minimal increase of the computational costs, as the mostexpensive part of (3) is the computation of the hydrodynamic forces F h,p . Having improved the accuracy of the contact model, we carried out simu-lations of particle-wall collisions in a fluid to compare to the experiments ofGondret et al. (2002). The details of the simulations, including the materialproperties as well as the physical and numerical parameters, are summarized inTable 1. Gondret et al. (2002) released particles from heights large enough toaccelerate to their terminal velocities before colliding with the wall. For thesesimulations, the horizontal wall and vertical particle trajectories allow us toonly consider normal collision forces. To control the impact velocity u in , weaccelerated the particle in the numerical simulations according to the relation u ( t ) = u in (cid:0) e − t − (cid:1) , ζ n > R p . (20)19 .
000 0 .
002 0 .
004 0 .
006 0 .
008 0 .
010 0 .
012 0 . t || F c || Exact solutionForward EulerBackward EulerPredictor-corrector
Figure 4:
Collision forces vs. time for Forward Euler, Backward Euler, and predictor-corrector schemes.
In other words, we prescribed the falling velocity of the particle so that it accel-erated in a smooth manner so that u in matched the Stokes number reported inGondret et al. (2002) as shown in Table 1. Two scenarios were considered: onewith a rather high Stokes number St = 152 and one with a lower Stokes numberof St = 27, the latter of which is within the range of Stokes numbers that havebeen reported for the numerical simulations of Kempe et al. (2014). Once theparticle reached a distance of ζ n = R p , we turned off the prescribed velocity, al-lowing the particle to move on its own volition according to the hydrodynamic,buoyant, and collision forces acting on it.While attempting to reproduce the experimental trajectories, the simula-tions produced large variations in the results from small changes to the initialconditions. To show this, we executed five simulations for St = 27, varyingonly the initial position of the particle y from the value recorded in Table 1within the interval of one grid cell h . Figure 5a shows the range of trajectoriesencountered. For St = 27, a substantial variation in the rebound height of upto 83% can be observed. We would thus expect the collision model to producea variety of incorrect trajectories for the simulation of sediment transport in ahorizontal channel flow. Even the mean of the variety of trajectories is not able20t 27 152Re p
30 164 R p (m) 0.003 0.0015 u in (m/s) 0.518 0.585 ρ p /ρ f ν f (m /s) 1 . × − . × − e dry g (m/s ) 9.81 9.81Domain size (m) ( L x × L y × L z ) 0 . × . × .
08 0 . × . × . × ns × p p × ns × pInitial position of sphere center (m) 0.075 0.197Grid cells in x -direction 256 128Grid cells per diameter 19 19Timestep ∆ t = 2 . t = 8 . Table 1:
Simulation parameters to match the experiments of Gondret et al. (2002).Boundary conditions can be periodic (p), slip (s), or no-slip (ns). to fully reproduce the experimental trajectory.To better understand the observed variability, we plot the time evolution ofthe collision forces, i.e. lubrication and contact forces, for the low Stokes numbercase St = 27 in Figure 6. In this plot, we can see the particle approaching thewall with the lubrication force growing as 1 /ζ n (phase I). Subsequently, thelubrication forces become zero during the contact phase starting at t/t ref = − .
1. During this phase (phase II), the contact force grows and then decays withthe particle-wall overlap as the particle changes direction to rebound. Finally,the particle experiences the lubrication force again during the rebound phasestarting at t/t ref = 0 (phase III). At this time, lubrication is acting in theopposite direction because lubrication is dissipative. The dotted line in Figure 6shows the forces acting on the particle for a time discretization based on CFL =0 . . . . . . . . t/t ref − . . . . . . . . . ζ n / D p Gondret et al.MeanBounds − . . . . . . . t/t ref − . . . . . . . . . ζ n / D p Gondret et al.MeanBounds a) b)
Figure 5:
Sensitivity of rebound trajectories to initial position y for St = 27 . a) Tra-jectories computed without particle substeps and b) trajectories computed with particlesubsteps. modifications described in Section 4.2 above is able to give a smooth evolutionof contact forces with the time step size of the fluid solver. However, it turns outthat the lubrication forces remain under-resolved during approach and rebound,especially as ζ n approaches zero directly before and after the contact phase.This leads to either more or less total impulse acting on the particle, dependingon where the timestep happens to land, which in turn results in variabilitybetween simulations. This effect strongly depends on the Stokes number, sincethe lubrication force decreases with increasing St. Hence, the ratio of the normalcontact force to the lubrication force increases when approaching dry contactconditions.Since the lubrication model used is an algebraic relation that does not dependon the surrounding hydrodynamics, we have implemented a substepping methodthat integrates the particle motion with smaller timesteps than the fluid motion.This method works as follows:1. We solve the fluid equations of motion, IBM, and hydrodynamic forcesacting on the particle as normal.2. We divide the fluid RK substep k into a number of substeps N sub,k = { , , } . This choice results in a total of 15 substeps of constant sizeper fluid timestep (∆ t sub = ∆ t/ α k =22 . − . − . − . − . − .
02 0 .
00 0 . t/t ref − − F / W b u o y I II III
With substepsWithout substeps
Figure 6:
Collision forces acting on the particle at
St = 27 including: the lubricationforce during approach (phase I), the normal force during contact (phase II), and thelubrication force during rebound (phase III). Vertical dashed lines indicate a change inphase. W buoy = (1 − ρ f /ρ p ) m p g is the buoyant weight of the particle. { / , / , / } as used in (19).3. For each of the substeps, we solve the particle equations of motion withthe three-step RK method. As we update the particle velocities and po-sitions, we re-evaluate the collision (lubrication and contact) forces, butthe hydrodynamic forces remain constant. This compromise makes thepresent approach very efficient.4. At the end of the 8, 2, or 5 substeps, we use the final particle position andvelocity for the next fluid RK substep.This measure effectively increases the resolution of a collision to a timestep 15times smaller than the fluid timestep to integrate particle motion, allowing us tocompute the lubrication forces with higher accuracy. Since the contact durationof T c = 10∆ t is maintained, the contact phase is now resolved with a total of 150timesteps with only a marginal increase to the computational cost. Substeppinghas also been used by Kidanemariam and Uhlmann (2014) and Costa et al.(2015) but the authors did not illustrate the variability we have observed for thetrajectories of particle-wall collisions. Meanwhile, Kidanemariam and Uhlmann232014) do not provide a comparison with the data of Gondret et al. (2002) at all.The scheme presented by Costa et al. (2015) still relies on an iterative proceduresubdividing every fluid timestep into 50 substeps, which is less efficient thanthe scheme presented here. The results of our approach can be appreciated inFigure 6. The solid line, which was resolved with fifteen times more timesteps,can be viewed as a better approximation of the exact solution to the modelwe have implemented. Figure 5b shows how this method almost eliminates thevariability in the rebound trajectories of the particle-wall collisions discussedabove. As shown in Figure 5, the improved integration scheme described in Sections4.2 and 4.3 yields excellent results in reproducing the rebound trajectory ofthe St = 27 experiment of Gondret et al. (2002). Having obtained consistentresults that are insensitive to the initial condition, we can use the same setupof particle-wall collisions to select the most suitable surface roughness ζ n,min for the lubrication model (9). This is the only parameter involved that requirescalibration as an inverse problem. However, the range of values that can beassigned to ζ n,min should neither fall below the surface roughness of the actualsimulated particle nor exceed the length of a grid cell in order for the lubricationmodel to make physical sense.The impact of ζ n,min on particle rebound trajectories for St = 27 is illus-trated in Figure 7. A clear trend can be identified: decreasing the value of ζ n,min also decreases the rebound height due to more damping within the lu-brication layer. However, the results are moderately sensitive to the roughnessvalue. For instance, note that changing the roughness by an order of magnitudehas a similar effect to excluding substeps (as shown in Figure 5). Based on thepresent results, we selected ζ n,min = 3e-3 R p to optimize agreement with theexperimental data. We have used this value for all simulations in the presentwork. Note that a surface roughness of 1e-4 R p has been reported by Gondretet al. (2002), and other authors have used the physical particle roughness length24 . . . . . . . t/t ref − . . . . . . . . . ζ n / D p Gondret et al. ζ n,min = 2e-3 R p ζ n,min = 3e-3 R p ζ n,min = 4e-3 R p ζ n,min = 6e-3 R p ζ n,min = 1e-2 R p Figure 7:
Effect of changing ζ n,min on rebound trajectories for St = 27 . for this parameter (Kempe and Fr¨ohlich, 2012a; Costa et al., 2015) to avoid thesingularity in the lubrication force. Thus, we do not consider this parameter tobe an exact physical representation of the actual surface roughness, but ratheras a parameter to be calibrated within a reasonable range (small enough to bemeaningful relative to the particle size and large enough to be resolved by thesubsteps). Finally, we present a clarification to the ACTM as written by Kempe andFr¨ohlich (2012a). As already mentioned in Section 3.2, the ACTM assumesthat (11) represents the equation of motion for the particle in determining thecoefficients k n and d n . In other words, no fluid or gravitational forces act on theparticle during the contact phase. Though not stated in their paper, Kempe andFr¨ohlich (2012a) excluded hydrodynamic and buoyant weight forces in order toreproduce the trajectories of Gondret et al. (2002) (Kempe & Fr¨ohlich, 2016,private communication). Thus, during contact the non-disabled Lagrangianmarkers still affect the fluid, but not the particle momentum. This procedure issomewhat delicate for the situation of sediment transport in a horizontal channelconsidering the fact that the governing nondimensional number is the ratio of the25 t/t ref . . . . . . . . . ζ n / D p Gondret et al.Dry contactWet contact . . . . . t/t ref . . . . . . ζ n / D p Gondret et al.Dry contactWet contact a) b)
Figure 8:
Effect of including (“wet”) or excluding (“dry”) fluid forces during contacton rebound trajectories. a)
St = 152 and b)
St = 27 . hydrodynamic stress to the buoyant weight of the particle. This characteristicnumber is classically known as the Shields parameter Sh = τ w / (( ρ p − ρ f ) gD p ),where τ w is the wall shear stress (Shields, 1936). It is, therefore, very desirableto include gravitational and buoyant forces in (3) during the contact phase. Inour experience, including the gravitational force during contact has a negligibleeffect in changing the desired T c and e dry . In fact, the results presented so farhave all been generated by including buoyant weight during contact.On the other hand, including the fluid forces during contact can lead tosignificant drag on the particle throughout the collision. Figure 8 shows howexcluding fluid forces during contact gives us excellent agreement with the ex-perimental results, while including fluid forces during contact leads to excessivedamping because the fluid surrounding the particle has not been able to adaptto the change of the kinematics of the particle. Indeed, the simulations of Sime-onov and Calantoni (2012) show coefficients of restitution below experimentalvalues for moderate Stokes numbers (20 < St < (cid:29)
1, redefining (3) as follows: m p d u p d t = F h,p + F g,p + F c,p ζ n > F g,p + F c,p ζ n ≤ . (21)26osta et al. (2015) implemented a similar method for particle-wall collisions,but they turned off fluid forces when the collision overlap exceeded the expectedoverlap due to the particle’s weight, i.e. ζ n < − (1 − ρ f /ρ p ) gm p /k n . For thecases shown in Figure 8, the timescale of the contact phase is much smallerthan the timescale of the general fluid flow, i.e. the timescale of the particlerebound. Thus, while neglecting fluid forces has an important effect on realizingthe correct e wet , it has a minimal effect on the general flow.However, neglecting fluid forces can lead to unphysical situations for endur-ing contact, which we define to be when the timescale of contact matches orexceeds that of the general flow. Consider, for example, a single particle at restand in contact with a wall. If we then impose a shear flow over the particle, itshould be swept up into the flow, or at the very least be carried downstream.However, in a simulation using (21), because the particle is in contact with thewall, it does not experience the hydrodynamic forces. It will therefore continueto sit on the wall, oblivious to the flow around it, until another particle collideswith it. This was addressed in Kempe et al. (2014) by switching on the hydro-dynamic forces for all collisions regardless of the Stokes number, even though itwas not explicitly mentioned in this reference (Kempe & Fr¨ohlich, 2016, privatecommunication). We address this problem in detail in the subsequent Section5 to introduce a suitable threshold for the inclusion of the hydrodynamic forcesin (3) and (4).
5. Enduring contact model
As shown in the results from Section 4, neglecting fluid forces acting onthe particle during contact produces a good match with the experimental dataof Gondret et al. (2002), which involve collisions of finite duration. However,problems can arise in the limit of enduring contact. We therefore propose toinclude fluid forces during contact below some threshold Stokes number St crit .For collisions above St crit , the contact duration should be finite ( T c = 10∆ t ) so27hat no major loss of physicality is encountered. For collisions below St crit , theparticle is not going to experience an appreciable rebound so that the particlemotion is not governed by collision forces during contact, but by hydrodynamicforces. Neglecting hydrodynamic forces in the low-Stokes number regime intro-duces artifacts in particle mobility. Indeed, this was observed in Vowinckel et al.(2016) for the situation of a horizontal turbulent open-channel flow laden withparticles heavier than their critical threshold of motion. Using the same methodfor collisions, these particles formed a closed bed of resting particles. In thisreference, it was shown that a collision with a fast moving particle was necessaryfor almost all of the erosion events recorded to dislodge a particle out of thesediment packing. However, it has not been possible to clarify to what extentthis triggering collision is merely a consequence of the collision procedure.To investigate what the critical value for the Stokes number may be, wecompared particle-wall collisions that include hydrodynamic forces during con-tact (“wet” contact) to those that exclude hydrodynamic forces during contact(“dry” contact), as illustrated in Figure 9. For this scenario, we used the sameparameters as those summarized in Table 1 (St = 27) and repeated the simula-tions for ever-decreasing St. The Stokes number was controlled by prescribingthe particle’s velocity until it made direct contact with the wall. Unlike theprevious simulations, we did not allow the lubrication layer to slow the particlebefore contact.For the cases with higher Stokes numbers, we can clearly see how includinghydrodynamic forces during contact leads to significant undershooting of therebound trajectory. As the Stokes number decreases, however, the significanceof this undershooting also decreases. For St < crit = 5.This value is consistent with the work of other researchers (Gondret et al., 2002;Joseph et al., 2001), who experimentally observed no rebounds for St <
10. Notethat the Stokes numbers reported in Figure 9 and our resulting St crit are basedon the particle velocity at contact, i.e. when ζ n = 0, whereas most other authors28 . . . . . t/t ref − . . . . . . . ζ n / D p Dry contactWet contact . . . . . t/t ref − . . . . . . . ζ n / D p Dry contactWet contact . . . . . t/t ref − . . . . . . . ζ n / D p Dry contactWet contact . . . . . t/t ref − . . . . . . . ζ n / D p Dry contactWet contact a) b)c) d)
Figure 9:
Comparison of trajectories of particle-wall collisions without (“dry”) and with(“wet”) hydrodynamic forces during contact for various Stokes numbers a) St = 20 ,b) St = 15 , c) St = 10 , d) St = 5 , where u in is measured at ζ n = 0 . report Stokes numbers at some distance from the wall, before the lubricationlayer has fully slowed the particle. With this enduring contact model, we canexpand the particle equation of motion (21) to m p d u p d t = F h,p + F g,p + F c,p , ζ n > F g,p + F c,p , ζ n ≤ ∧ max { St } > St crit F h,p + F g,p + F c,p , ζ n ≤ ∧ max { St } ≤ St crit , (22)where the max { St } function represents the maximum Stokes number among allactive collisions for particle p , and fluid forces acting on the particle are onlyincluded from non-disabled markers. The same consideration applies for theangular momentum (4). Using this scheme now allows us to include the fullmomentum balance for particles in enduring contact, i.e. the hydrodynamicstresses as well as the buoyant weight of the particle, so that the considerations29f the Shields parameter become applicable. In the case of St (cid:28)
1, the impact velocity u in approaches zero. This meansin turn that the computed stiffness in (B.1d) would approach infinity. Thisproblem is addressed by Kempe and Fr¨ohlich (2012a) who have introduced acritical Stokes number St crit , which establishes a minimum impact velocity tolimit k n for enduring contact: u in,crit = 9St crit ρ f ν f ρ p D p . (23)In the present study, this critical impact velocity was used in (B.1d) and (B.1c)to compute k n and d n , respectively, for such collisions. This implementationdiffers slightly from that of Kempe and Fr¨ohlich (2012a), who do not apply anydamping for collisions with St < St crit , i.e. they have set d n = 0. We includedthis damping for enduring contact in order to help reach steady-state conditions.Our implementation also differs in that we use St crit = 5 whereas Kempe andFr¨ohlich (2012a) used St crit = 1.Furthermore, we retain the buoyant weight forces in the equation of motionduring contact as outlined in Section 4.5. This means that, for particle packingsseveral diameters thick, the weight of a single sphere resting on another layerof particles is passed along to deeper layers. This effect enhances the physicalrealism because frictional contact forces increase with depth, but it also resultsin increasing particle surface overlap with depth and ultimately in a changeof porosity of the sediment bed, which has been acknowledged as a crucialparameter to define the hydraulic resistance of a sediment to the flow (Vowinckelet al., 2014). However, a flow with a lower Reynolds number would resultin collisions with lower Stokes numbers such that u in,crit could become largerelative to the particle size and relevant time scales. A large u in,crit would resultin a low k n and hence a large overlap between particles, which is undesirable. Toprevent this large overlap, we enforce a maximum overlap distance (cid:15)R p throughthe following procedure: for a collision with St < St crit , the stiffness is given by30 n = m eff √ u in t ∗ u in > u in,crit max ( k n,crit , k n,grav ) u in < u in,crit (24)where m eff and t ∗ are defined in (A.2) and (B.1b), respectively, k n,crit = m eff (cid:112) u in,crit t ∗ (25)is the stiffness limited by the critical impact velocity, and k n,grav = max (cid:104) m p g ( (cid:15)R p ) − / , m q g ( (cid:15)R q ) − / (cid:105) (26)is the stiffness required for particle p (or q ) to have a steady-state overlap of (cid:15)R p (or (cid:15)R q ) with a wall due to gravity. To have a minimal constant overlap weset (cid:15) = 10 − . Thus, we ensure that a bed of particles contains a uniform set ofcollision stiffnesses that minimize particle overlap. The coefficient of friction for a material can depend on whether the contactis rolling or sliding (Fishbane et al., 1996). The rolling condition implies zeroslip at the contact point, i.e. (cid:107) g t,cp (cid:107) = 0 (cf. Appendix A). As a consequence,particle surfaces are in sticking contact for rolling motion until a critical thresh-old of static friction F s = µ s (cid:107) F n (cid:107) is exceeded, where µ s is the coefficient ofstatic friction. As soon as this condition is met, significant slip occurs and thecontact condition changes from sticking to sliding, so that the threshold for ki-netic friction F k = µ k (cid:107) F n (cid:107) must be used, where µ k is the coefficient of kineticfriction, with µ s always greater than µ k . Apart from the physical reasoningpresented above, limiting the frictional forces also becomes important from anumerical point of view whenever two or more collision partners are involved.Otherwise the multiple contact points competing for no-slip conditions can leadto instabilities in the calculation of the frictional forces.In the present study, the distinction between rolling/sticking and sliding ismade by the following scheme, which is comparable to that of Luding (2008):31ase Dry oblique collision Rolling in shear flow R p ρ p /ρ f e dry ν µ k µ s g ν f t = 2e − Table 2:
Simulation setup for oblique and rolling sphere simulations. • While the particle is sticking, i.e. || F t,LS || < || µ F n || , we set µ = µ s totest for the onset of slipping. • Once slipping occurs, i.e. || F t,LS || > || µ F n || , we set µ = µ k until thefriction force falls below the Coulomb friction force.The aim of the present study is to simulate natural sediment. Hence we parametrizedthe coefficients of friction with typical values of silicate materials, yielding µ k = 0 .
15 based on the work of Joseph and Hunt (2004), who worked withglass spheres, and µ s = 0 . ψ out = u t,out u n,in , (27)32 . . . . . . . . ψ in − . − . . . . . . . . ψ o u t Foerster et al.Present
Figure 10:
Rebound vs. impact angles for a particle-wall oblique collision. which depends on the impact angle ψ in = u t,in u n,in . (28)Here, u n,in is the impact velocity normal to the wall, while u t,in and u t,out arethe impact and rebound velocities, respectively, of the particle’s contact pointtangential to the wall ( u t = u p + R p ω p,z for a particle obliquely colliding inthe x -direction). For a particle with no initial rotation, ψ in is the tangent ofthe angle the particle makes with the wall from the normal ( ψ in = 0 meansno relative tangent motion). The rebound angle is zero when the contact issticking perfectly at the time of release. However, the rebound angle is negativewhen, at the point of release, | u p | < | ω p,z | (since ω p,z < ψ out at low impact angles.To test both situations, rolling and sliding, we simulated a particle in aCouette flow. We placed a sphere of radius R p /H = 0 . H . We initialized the particle at rest at a distance ζ n /R p = 1 . × − above the bottom wall. We subsequently exposed the sphereto a linear shear flow, holding it fixed for a short time ( tU/H = 0 .
01) to allowthe flow to develop around it before releasing it. The numerical parameters33 . . . . . . . . . t U w /H . . . . . . u / U w u p − ω p,z R p,cp . . . . . . . . . t U w /H . . . . . . u / U w u p − ω p,z R p,cp a) b) Figure 11:
Translational and rotational velocities of a particle exposed to a linear shearflow. a) sliding motion for Re H = 10 and b) rolling motion for Re H = 50 . are summarized in Table 2. We found that slipping motion occurs for a lowerReynolds number of Re H = U H/ν f = 10, where U is the lid velocity. On theother hand, perfect rolling motion occurs at Re H = U H/ν f = 50.Figure 11 shows how the particle accelerates until it reaches a steady-statetranslational velocity. As soon as the particle makes contact with the wall,gravity holds it there with a slight overlap according to the conditions definedin Section 5.2. Accounting for fluid forces during contact allows the particle toaccelerate to a steady-state speed while in contact with the wall. As expectedthe particle achieves perfect rolling without slip (Figure 11b), marked by thematch between the translational velocity u p and the rotational velocity relativeto the particle center − ω p,z R p,cp . Accordingly, the particle shows significantslip for the lower Reynolds number (Figure 11a), where the increased viscosityleads to increased drag on the particle, which in turn overpowers the frictionfrom the particle’s weight.
6. Flow over dense sediment
We presented a detailed validation of binary particle-wall collisions in Sec-tions 4 and 5. To address the bulk behavior of a dense granular bed shearedby a laminar Poiseuille flow, we carried out numerical simulations to reproduce34 igure 12:
Left plot: Instantaneous snapshot of a slice through the xy -plane for caseA10M. Contours show the streamwise component of the fluid velocity. Particles arecolored grayscale according to their velocity. Right plot: Streamwise and spanwiseaverages of fluid and particle velocities. Arrows correspond to the length scales for theclear fluid, h f , the particle bed, h p , the mobile bed layer, h m , and the motionless bedlayer, h c . the experimental results of Aussillous et al. (2013), who studied pressure-drivenflows over glass spheres with a mean diameter D p = 1 . σ ( D p ) = 0 . h f /D p and Reynolds numbers inthe laminar regime, where h f is the height of the clear-water layer above thesediment bed illustrated in Figure 12. We define h f to be the height abovewhich the average particle volume fraction φ < .
05, which is the threshold fornegligible impact of particle-particle interaction on the flow (Capart and Frac-carollo, 2011). We define the mobile bed height h m to be the portion of theparticle bed above which the mean particle velocity is higher than 1% of thevalue at the fluid/particle interface. 35a 0.397 R p (m) 0.0444 ρ p /ρ f ν f (m /s) 0.219 g (m/s ) 9.81 e dry ζ n,min (m) 3 . × − R p µ k µ s L x × L y × L z ) 1 . × . × . × ns × pGrid cells in x -direction 256 D p /h . Table 3:
Simulation parameters to match the experiments of Aussillous et al. (2013).Boundary conditions are periodic (p) and no-slip (ns). The Galileo number Ga isdefined in (29) . In their experiments, Aussillous et al. (2013) filled a long chamber with par-ticles and then applied a constant pressure gradient, which eroded the particlesfrom the chamber. Initially, the fluid height h f was small and the pressuregradient drove a large number of particles so that the height of the mobile bedlayer, h m in Figure 12, was large. Since no new particles were added to thechamber, h f increased as the particles eroded away until, at long periods oftime, the experiment reached a steady-state configuration where the influx ofparticles into the observation window remained in equilibrium with the outflux.Due to our use of periodic boundary conditions, we will only try to replicatethe long-term steady-state flow conditions, of which there are only a few datafrom Aussillous et al. (2013) 36ase A1 A2 A9 A10 A10MRe b h f /D p (Exp.) 7 . ± . . ± . . ± . . ± . . ± . h f /D p (Sim.) 7 .
15 8 .
31 10 .
33 11 .
29 11.05Sh (Exp.) 0 . ± .
03 0 . ± .
03 0 . ± .
04 0 . ± .
03 0 . ± . .
224 0 .
222 0 .
358 0 .
343 0.357 q f (m /s) 0.0659 0.0880 0.220 0.251 0.251 N p,m ∗ N p,f
132 132 132 132 132 T avg (s) 139.5 137.9 126.2 127.3 111.0 Table 4:
Parameters that vary between the different cases. The bulk Reynolds number Re b is defined in (30) and the Shields number Sh is defined in (31) . The fluid height h f (and hence Shields number) do not exactly match between the experiments (Exp.) andsimulations (Sim.). ∗ Polydisperse particle diameters follow a Gaussian distributionwith a standard deviation of σ ( D p ) = 0 . D p . We executed several simulations in an attempt to match four of the exper-iments of Aussillous et al. (2013) at different flow rates and fluid heights. Thephysical and numerical parameters associated with these simulations are listedin Table 3, and the differences between the four cases are listed in Table 4.These experiments can be characterized by the Galileo numberGa = (cid:113) ( ρ p /ρ f − gD p ν f , (29)the bulk Reynolds number Re b = q f ν f , (30)where q f is the fluid flow rate, and the Shields numberSh = 6Re b Ga (cid:18) D p h f (cid:19) , (31)which represents the ratio of the shear stress acting on the particle bed to thebuoyant weight of a particle. Aussillous et al. (2013) reported an uncertainty37or the determination of the bed height as h f ± R p , which we have included inTable 4 as the deviations in h f and Sh, which depends on h f .We required a low CFL = 0 . D p /h = 22 . D p /h = 17 . N p,m particles in a computational domain with periodic x - and z -boundaries above a layer of N p,f fixed particles, which we arranged ina hexagonal packing with random heights varying from 0 < y < D p . Thesefixed particles were used to avoid over-idealized smooth conditions at the lowerwall. We subsequently allowed the non-fixed particles to settle under “dry”conditions, i.e. without considering hydrodynamic forces. We then applied alarge pressure gradient to produce a fluid flow rate 8 times that of the finaldesired flow rate, mobilizing the entire bed. This mobilization also caused thebed to dilate, or have the average local volume fraction decrease, which in turndecreased h f . Once h f dropped to about 0 . D p below the desired value, weimmediately decreased the flow rate to the final flow rate reported in Table 4,which is defined as q f = 1 L x L z (cid:90) L z (cid:90) L y (cid:90) L x (1 − φ ) u d x d y d z , (32)where φ is the particle volume fraction. We adopted this procedure because wenoticed a hysteresis in the particle flux between an increased flow rate and adecreased flow rate, which has also been observed by Clark et al. (2015). Notethat this procedure more closely resembles the experiments, where the particlebed is largely mobilized and then settles into a lower particle flux.However, one problem with this procedure is that we cannot determine thefinal bed height a priori . The dilation and contraction accompanying the two38 .
00 0 .
05 0 .
10 0 .
15 0 . u/ ( q f /h f ) − . − . − . − . − . . ( y − h p ) / h f FluidParticlesAussillous et al. .
00 0 .
05 0 .
10 0 .
15 0 . u/ ( q f /h f ) − . − . − . − . − . . ( y − h p ) / h f FluidParticlesAussillous et al. .
00 0 .
05 0 .
10 0 .
15 0 . u/ ( q f /h f ) − . − . − . − . − . . ( y − h p ) / h f FluidParticlesAussillous et al. .
00 0 .
05 0 .
10 0 .
15 0 . u/ ( q f /h f ) − . − . − . − . − . . ( y − h p ) / h f FluidParticlesAussillous et al. a) b)c) d)
Figure 13:
Wall-normal profiles of average fluid and particle velocities near the parti-cle/fluid interface compared to the wall-normal particle velocity profile from Aussillouset al. (2013). a) Case A1, b) Case A2 c) Case A9 d) Case A10. flow rates is difficult to predict without executing an iterative procedure of run-ning simulations with varying numbers of particles. Due to the computationalcosts of the simulations, we did not iterate on this method, but instead ac-cepted the values we obtained for h f , which, with the exception of case A10M,are larger than those in the experiments, as seen in Table 4. In Figure 13 we compare the particle velocity profiles of the simulation tothe experimental results of Aussillous et al. (2013). We calculated the particlevelocity profile u p ( y ) from our simulations by averaging the velocities of all theparticles in the streamwise and spanwise directions whose center fell within agiven range of heights. We used bins of width R p arranged such that the topmostbin extended from y = h p − R p to y = h p . We evaluated the fluid velocity profile39y averaging the u -velocity field in the streamwise and spanwise directions foreach grid cell of the y -coordinate. For this calculation, we used the particle cellvolume fractions φ to exclude fluid velocities existing within the particles: (cid:104) u (cid:105) xz = (cid:82) L z (cid:82) L x (1 − φ ) u d x d z (cid:82) L z (cid:82) L x (1 − φ ) d x d z (33)The fluid velocity profiles exhibit a parabolic shape in the clear fluid above thebed, as shown in Figure 12. At the interface between the clear fluid and particlebed, we observe some slip between the fluid and the particles, but within thebed the two velocity profiles are very similar, with only a slight difference dueto flow between the particles. The particle velocity profiles from the simulationscompare very well with the experiments for cases A9 and A10, and reasonablywell for cases A1 and A2.Part of the discrepancy between our experiments and the simulations is dueto the differences in bed heights and Shields numbers as seen in Table 4. Inthis table, we can see that cases A1 and A2 exhibit the largest differences in thefluid height between the simulations and experiments, which may have resultedin the larger deviations in the velocity profiles seen in Figure 13. Likewise, forthese two cases we can also see larger differences in the Shields number, whichcan be sensitive to the fluid height h f . We ran the simulation until it reached a constant particle velocity flux q v ,defined as q v = (cid:90) L y u p ( y ) d y , (34)where u p ( y ) is the particle velocity profile as defined in the previous section.Unlike q f , which had no variability in time, q v did vary as particles occasionallylocked in place or rolled over one another. We therefore evaluated a time-averaged value of the particle velocity flux (cid:104) q v (cid:105) T = 1 T avg (cid:90) t f t s q v d t , (35)40ase Experimental value Simulation value (cid:104) q v (cid:105) − T /q f (cid:104) q v (cid:105) T /q f (cid:104) q v (cid:105) + T /q f (cid:104) q v (cid:105) T /q f A1 5.e-3 1.2e-2 1.9e-2 6.71e-3A2 5.e-3 1.1e-2 1.6e-2 5.17e-3A9 4.87e-3 8.20e-3 1.15e-2 7.61e-3A10 5.17e-3 7.04e-3 8.91e-3 5.84e-3A10M 5.17e-3 7.04e-3 8.91e-3 7.56e-3
Table 5:
Comparison of the velocity flux (cid:104) q v (cid:105) T between our simulations and the experi-ments of Aussillous et al. (2013). (cid:104) q v (cid:105) + T /q f and (cid:104) q v (cid:105) − T /q f represent the mean (cid:104) q v (cid:105) T /q f plus and minus the standard deviation over the averaging time, respectively. where t f is the time at the end of the simulation, t s is the time at which theparticle flux reached steady-state, and T avg = t f − t s is the time interval overwhich the data was averaged. The values of T avg are given in Table 4.In Table 5, we can see a good agreement between our numerical results andthe experimental values of the velocity flux q v . Because these quantities arederived from the particle velocity profiles, we expect to see the similar trends,namely that we underestimate the mean values from the experiments and obtainbetter matches for cases A9 and A10. However, our results still fall within themargin of error of the experiments. Furthermore, we conducted another simulation to show the effect of poly-dispersity. In experiments, it is impossible to have a perfectly monodisperse setof particles. In their article, Aussillous et al. (2013) reported having a set ofspheres with diameters following a Gaussian distribution of mean D p = 1 . σ ( D p ) = 0 . . . . . . . . . φ − . − . − . − . − . . . ( y − h p ) / h f A10A10M
Figure 14:
Wall-normal profiles of average particle volume fractions.
In Figure 14, we do not see any appreciable changes in the particle bedvolume fractions between cases A10 and A10M. The average volume fractionwithin the bed is φ = 0 . h f compared to that of A10, which results in a higher Shieldsnumber, as shown in Table 4. Therefore, we also obtain a velocity flux thatoverpredicts the mean experimental value, as shown in Table 5. However, theparticle velocity profile and velocity flux still agree very well with the exper-imental results, and the results suggest that using monodisperse spheres is avalid approximation to polydisperse spheres for this experimental setup.
7. Conclusions
In the present study, we presented and validated a contact model for thepurpose of phase-resolved Direct Numerical Simulations, in which the dispersephase is represented by the Immersed Boundary Method. The present modelingapproach allows for actual particle contact and takes all relevant contact forcesinto account without introducing parameters that require arbitrary calibration.42 .
00 0 .
05 0 .
10 0 .
15 0 . u/ ( q f /h f ) − . − . − . − . − . . ( y − h p ) / h f FluidParticlesAussillous et al.
Figure 15:
Wall-normal profiles of average fluid and particle velocities for the simulationwith polydisperse particles (case A10M).
These forces include lubrication forces for small inter-particle gaps, normal re-pulsive forces to resolve inelastic collisions, and tangential forces to representparticle friction. We demonstrated that an improved integration scheme is nec-essary to obtain consistent results for particle-wall collisions. Subsequently, wepresented enhancements that extend the model to deal with simulations of flowsover dense granular sediments. It turns out that these enhancements are cru-cial in order to deal with thick sediment packings. The measures taken allowus to generate sediment packings several diameters thick that are numericallystable as the packing reaches a steady-state condition. The simulations areperformed by retaining the full momentum balance of a particle in enduringcontact, which includes the hydrodynamic forces and the buoyant weight of aparticle. Including these forces is crucial to represent phenomena like erosionand resuspension of particles. Moreover, the enhanced model allows for rollingand sliding contact, distinguishing between sticking and sliding conditions.Altogether, the present approach yielded satisfactory agreement with thebenchmark test cases for binary collisions as well as the collective motion ofparticles for a horizontal flow over a dense granular packing. In addition, afirst test case involving polydisperse sediment was presented. The high degree43f accuracy achieved is going to enable us to analyze phase-resolved numericalsimulation data in great detail. Although not explicitly stated, we believe thatthe present approach is also applicable for the situation of vertical channel flowsas well as neutrally buoyant particles in laminar and turbulent conditions. Itcan therefore provide a valuable tool to generate high-fidelity data even on thegrain scale of any kind of multiphase flows involving rigid spheres.
8. Acknowledgements
This research is supported in part by the Department of Energy Office ofScience Graduate Fellowship Program (DOE SCGF), made possible in part bythe American Recovery and Reinvestment Act of 2009, administered by ORISE-ORAU under contract no. DE-AC05-06OR23100. It is also supported by thePetroleum Research Fund, administered by the American Chemical Society,grant number 54948-ND9. BV gratefully acknowledges the Feodor-Lynen schol-arship provided by the Alexander von Humboldt foundation, Germany, andEM thanks Petrobras for partial support. The authors thank ´E. Guazelli and P.Aussillous for stimulating discussions on their experimental work. P. Gondret,N. Mordant, and A. Ten Cate are acknowledged for kindly providing their data.Computational resources for this work used the Extreme Science and Engineer-ing Discovery Environment (XSEDE), which was supported by the NationalScience Foundation, USA, Grant No. TG-CTS150053.
ReferencesReferences
Apte, S.V., Martin, M., Patankar, N.A., 2009. A numerical method for fullyresolved simulation (FRS) of rigid particleflow interactions in complex flows.Journal of Computational Physics 228, 2712–2738.Aussillous, P., Chauchat, J., Pailha, M., M´edale, M.and Guazzelli, E., 2013.Investigation of the mobile granular layer in bedload transport by laminarshearing flows. Journal of Fluid Mechanics 736, 594–615.44alachandar, S., Eaton, J.K., 2010. Turbulent dispersed multiphase flow. An-nual Review of Fluid Mechanics 42, 111–133.Breugem, W.P., 2012. A second-order accurate immersed boundary method forfully resolved simulations of particle-laden flows. Journal of ComputationalPhysics 231, 4469–4498.Capart, H., Fraccarollo, L., 2011. Transport layer structure in intense bed-load.Geophysical Research Letters 38, 2–7.Charru, F., Andreotti, B., Claudin, P., 2013. Sand ripples and dunes. Annu.Rev. Fluid Mech. 45, 469–493.Clark, A.H., Shattuck, M.D., Ouellette, N.T., O’Hern, C.S., 2015. Onset andcessation of motion in hydrodynamically sheared granular beds. PhysicalReview E - Statistical, Nonlinear, and Soft Matter Physics 92, 1–7.Costa, P., Boersma, B.J., Westerweel, J., Breugem, W.P., 2015. Collision modelfor fully resolved simulations of flows laden with finite-size particles. PhysicalReview E 92, 053012.Cox, R., Brenner, H., 1967. The slow motion of a sphere through a viscous fluidtowards a plane surface. Small gap widths, including inertial effects. ChemicalEngineering Science 22, 1753–1777.Derksen, J.J., 2011. Simulations of granular bed erosion due to laminar shearflow near the critical shields number. Physics of Fluids 23, 113303.Derksen, J.J., 2015. Simulations of granular bed erosion due to a mildly turbu-lent shear flow. Journal of Hydraulic Research 53, 622–632.Dieterich, J.H., 1972. Time-dependent friction in rocks. Journal of GeophysicalResearch 77, 3690–3697.Fadlun, E.A., Verzicco, R., Orlandi, P., Mohd-Yusof, J., 2000. Combinedimmersed-boundary finite-difference methods for three-dimensional complexflow simulations. Journal of Computational Physics 161, 35–60.45enton, J.D., Abbott, J.E., 1977. Initial movement of grains on a stream bed- effect of relative protrusion. Proceedings of the Royal Society of LondonSeries A-mathematical Physical and Engineering Sciences 352, 523–537.Fishbane, P.M., Gasiorowicz, S., Thornton, S.T., 1996. Physics for scientistsand engineers. Pearson College Division.Foerster, S.F., Louge, M.Y., Chang, H., Allia, K., 1994. Measurements of thecollision properties of small spheres. Physics of Fluids 6, 1108.Fukuoka, S., Fukuda, T., Uchida, T., 2014. Effects of sizes and shapes of gravelparticles on sediment transports and bed variations in a numerical movable-bed channel. Advances in Water Resources 72, 84–96.Garcia, M., 2008. Sediment transport and morphodynamics, in: Garcia, M.(Ed.), American Society of Civil Engineers, Manuals and Reports on Engi-neering Practice 110, pp. 21–168.Glowinski, R., Pan, T.W., Hesla, T.I., Joseph, D.D., 1999. A distributed La-grange multiplier/fictitious domain method for particulate flows. Interna-tional Journal of Multiphase Flow 25, 755–794.Glowinski, R., Pan, T.W., Hesla, T.I., Joseph, D.D., Priaux, J., 2001. A ficti-tious domain approach to the direct numerical simulation of incompressibleviscous flow past moving rigid bodies: Application to particulate flow. JournalComputational Physics 169, 363–426.Gondret, P., Lance, M., Petit, L., 2002. Bouncing motion of spherical particlesin fluids. Physics of Fluids 14, 643–652.Hertz, H., 1882. ¨Uber die Ber¨uhrung fester elastischer K¨orper. J. f. reine u.angewandte Math. 92, 156–171.Izard, E., Bonometti, T., Lacaze, L., 2014. Modelling the dynamics of a sphereapproaching and bouncing on a wall in a viscous fluid. Journal of FluidMechanics 747, 422–446. 46oseph, G.G., Hunt, M.L., 2004. Oblique particle-wall collisions in a liquid.Journal of Fluid Mechanics 510, 71–93.Joseph, G.G., Zenit, R., Hunt, M.L., Rosenwinkel, A.M., 2001. Particle-wallcollisions in a viscous fluid. Journal of Fluid Mechanics 433, 329–346.Kajishima, T., Takiguchi, S., Hamasaki, H., Miyake, Y., 2001. Turbulencestructure of particle-laden flow in a vertical plane channel due to vortex shed-ding. JSME International Journal Series B-Fluids and Thermal Engineering44, 526–535.Kempe, T., Fr¨ohlich, J., 2012a. Collision modelling for the interface-resolvedsimulation of spherical particles in viscous fluids. Journal of Fluid Mechanics709, 445–489.Kempe, T., Fr¨ohlich, J., 2012b. An improved immersed boundary method withdirect forcing for the simulation of particle laden flows. Journal of Computa-tional Physics 231, 3663–3684.Kempe, T., Vowinckel, B., Fr¨ohlich, J., 2014. On the relevance of collisionmodeling for interface-resolving simulations of sediment transport in openchannel flow. International Journal of Multiphase Flow 58, 214–235.Kidanemariam, A.G., Chan-Braun, C., Doychev, T., Uhlmann, M., 2013. Directnumerical simulation of horizontal open channel flow with finite-size, heavyparticles at low solid volume fraction. New Journal of Physics 15, 025031.Kidanemariam, A.G., Uhlmann, M., 2014. Interface-resolved direct numericalsimulation of the erosion of a sediment bed sheared by laminar channel flow.International Journal of Multiphase Flow 67, 174–188.Lajeunesse, E., Malverti, L., Charru, F., 2010. Bed load transport in turbulentflow at the grain scale: Experiments and modeling. J. Geophys. Res.: EarthSurfarce 115. 47ucci, F., Ferrante, A., Elghobashi, S., 2010. Modulation of isotropic turbulenceby particles of taylor length-scale size. Journal of Fluid Mechanics 650, 5–55.Luding, S., 2008. Cohesive, frictional powders: contact models for tension.Granular Matter 10, 235–246.Meiburg, E., Kneller, B., 2010. Turbidity currents and their deposits. Annu.Rev. Fluid Mech. 42, 135–156.Mittal, R., Iaccarino, G., 2005. Immersed boundary methods. Annu. Rev. FluidMech. 37, 239–261.Mordant, N., Pinton, J.F., 2000. Velocity measurement of a settling sphere.The European Physical Journal B-Condensed Matter and Complex Systems18, 343–352.Ni˜no, Y., Garc´ıa, M., 1998. Experiments on saltation of sand in water. Journalof Hydraulic Engineering 124, 1014–1025.Patankar, N.A., Singh, P., Joseph, D.D., Glowinski, R., Pan, T.W., 2000. A newformulation of the distributed Lagrange multiplier/fictitious domain methodfor particulate flows. International Journal of Multiphase Flow 26, 1509–1524.Picano, F., Breugem, W.P., Brandt, L., 2015. Turbulent channel flow of densesuspensions of neutrally buoyant spheres. Journal of Fluid Mechanics 764,463–487.Rai, M.M., Moin, P., 1991. Direct simulations of turbulent flow using finite-difference schemes. Journal of computational physics 96, 15–53.Ray, S., Kempe, T., Fr¨ohlich, 2015. Efficient modelling of particle collisions usinga non-linear viscoelastic contact force. International Journal of MultiphaseFlow , accepted.Roma, A., Peskin, C., Berger, M., 1999. An Adaptive Version of the ImmersedBoundary Method. Journal of Computational Physics 153, 509–534.48antarelli, C., Fr¨ohlich, J., 2015. Direct numerical simulations of sphericalbubbles in vertical turbulent channel flow. International Journal of MultiphaseFlow 75, 174–193.Seminara, G., 2010. Fluvial sedimentary patterns. Annu. Rev. Fluid Mech. 42,43–66.Shao, X.M., Wu, T.H., Yu, Z.S., 2012. Fully resolved numerical simulationof particle-laden turbulent flow in a horizontal channel at a low Reynoldsnumber. Journal of Fluid Mechanics 693, 319–344.Shields, A., 1936. Anwendung der ¨Ahnlichkeitsmechanik und der Turbulenz-forschung auf die Geschiebebewegung. Ph.D. thesis. Mitteilungen der Preußis-chen Versuchsanstalt f¨ur Wasserbau und Schiffbau, Berlin (in German).Sierakowski, A.J., Prosperetti, A., 2016. Resolved-particle simulation by thephysalis method: Enhancements and new capabilities. Journal of Computa-tional Physics 309, 164–184.Simeonov, J.A., Calantoni, J., 2012. Modeling mechanical contact and lubri-cation in Direct Numerical Simulations of colliding particles. InternationalJournal of Multiphase Flow 46, 38–53.Ten Cate, A., Nieuwstad, C.H., Derksen, J.J., Van den Akker, H.E.A., 2002.Particle imaging velocimetry experiments and lattice-boltzmann simulationson a single sphere settling under gravity. Physics of Fluids 14, 4012–4025.Thornton, C., Cummins, S.J., Cleary, P.W., 2011. An investigation of the com-parative behaviour of alternative contact force models during elastic collisions.Powder Technology 210, 189–197.Thornton, C., Cummins, S.J., Cleary, P.W., 2013. An investigation of thecomparative behaviour of alternative contact force models during inelasticcollisions. Powder Technology 233, 30–46.49hlmann, M., 2005. An immersed boundary method with direct forcing forthe simulation of particulate flows. Journal of Computational Physics 209,448–476.Uhlmann, M., 2008. Interface-resolved direct numerical simulation of verticalparticulate channel flow in the turbulent regime. Physics of Fluids 20, 053305.Vowinckel, B., Jain, R., Kempe, T., Fr¨ohlich, J., 2016. Erosion of single particlesin a turbulent open-channel flow: a numerical study. J. Hydraul. Res. 54, 158–171.Vowinckel, B., Kempe, T., Fr¨ohlich, J., 2014. Fluid-particle interaction in tur-bulent open channel flow with fully-resolved mobile beds. Advances in WaterResources 72, 32–44.Zhang, Z., Prosperetti, A., 2005. A second-order method for three-dimensionalparticle simulation. Journal of Computational Physics 210, 292–324.Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, a.B., 2008. Discrete particle simulationof particulate systems: A review of major applications and findings. ChemicalEngineering Science 63, 5728–5770. 50 ppendixAppendix A. Definitions for particle-particle and particle-wall colli-sions
In order to discuss collisions in a general manner, we provide definitions forseveral variables that describe the contact. Some definitions will depend onwhether the interaction is between particle p and a wall (particle-wall interac-tion, or P-W) or between particle p and particle q (particle-particle interaction,or P-P). For most of the definitions, collisions between a fixed particle anda mobile particle are handled identically to particle-particle collisions, unlessindicated otherwise (particle-fixed, or P-F). • R eff – effective radius R eff = R p R q R p + R q (P-P) (A.1a) R eff = R p (P-W) (A.1b) • m eff – effective mass m eff = m p m q m p + m q (P-P) (A.2a) m eff = m p (P-W, P-F) (A.2b) • x w – point on wall closest to particle • n – unit vector normal to the surface of contact, points from x p to x q (P-P) or directly towards the wall (P-W) n = x q − x p || x q − x p || (P-P) (A.3a) n = x w − x p || x w − x p || (P-W) (A.3b) • ζ n – distance between surfaces of the two bodies (negative value indicatesoverlap) ζ n = || x q − x p || − R p − R q (P-P) (A.4a) ζ n = || x w − x p || − R p (P-W) (A.4b)51 x cp – location of contact point between surfaces, halfway between surfaceoverlap (P-P) x cp = x p + (cid:18) R p + ζ n (cid:19) n (P-P) (A.5a) x cp = x w (P-W) (A.5b) • R p,cp – radius of particle p with respect to the contact point R p,cp = || x cp − x p || (A.6) • g – relative velocity between particle centers of mass g = u p − u q (P-P) (A.7a) g = u p (P-W) (A.7b) • g cp – relative velocity of surface contact point g cp = g + R p,cp ( ω p × n ) + R q,cp ( ω q × n ) (P-P) (A.8a) g cp = g + R p,cp ( ω p × n ) (P-W) (A.8b) • g n,cp – component of g cp normal to surface g n,cp = ( g cp · n ) n (A.9) • g t,cp – component of g cp tangent to surface g t,cp = g cp − g n,cp (A.10) Appendix B. Calculating the normal contact model coefficients
In order to obtain the stiffness and damping coefficients k n and d n , Ray et al.(2015) use nonlinear transformations and a series expansion of (11) to yield thefollowing algebraic expressions: λ = 1 α τ c, (cid:32) − Cη + (cid:114) C η + α τ c, η (cid:33) , (B.1a)52 ∗ = T c τ c, (cid:112) − Aλ − Bλ , (B.1b) d n = 2 λm eff t ∗ , (B.1c)and k n = m eff (cid:112) u in t ∗ , (B.1d)where A = 0 . , B = 0 . C = 0 . α = 1 . τ c, = 3 .
218 areconstants. The parameter η = (ln e dry ) accounts for the restitution coefficient,and we measure the impact velocity to be u in = g n,cp · n at the first occurrenceof ζ n ≤ Appendix C. The tangential displacement vector
Tangential models based on spring systems require a displacement as definedby (13), which represents the accumulated relative motion between two surfacesWe calculate ζ t in a discrete sense as follows: (cid:101) ζ t = ζ k − t − (cid:16) ζ k − t · n (cid:17) n (C.1a) (cid:98) ζ t = || ζ k − t |||| (cid:101) ζ t || (cid:101) ζ t (C.1b) ζ kt = (cid:98) ζ t + 2 α k ∆ t g t,cp . (C.1c)Equations (C.1a) and (C.1b) rotate the displacement from the previous timesteponto a plane tangent to the two surfaces. Luding (2008) implemented thisrotation to account for the change in reference frame that can take place betweentwo timesteps. Without this rotation, the linear spring could contribute to thenormal force acting between two particles.Furthermore, when the two surfaces slip according to the Coulomb frictioncriteria, the displacement vector should not grow as the two surfaces continueto slide past one another. Instead, we reset the displacement to that whichachieves the Coulomb friction force: ζ t = − || µ F n || t + d t g t,cp k t if || F t,LS || > || µ F n || ..