Numerical study of a sphere descending along an inclined slope in a liquid

The descending process of a sphere rolling and/or sliding along an inclined slope in a liquid involves interactions between the hydrodynamic forces on the sphere and the contact forces between the sphere and the plane. In this study, the descending process of sphere in a liquid was examined using coupled LBM–DEM technique. The effects of slope angle, viscosity and friction coefficient on the movement of a sphere were investigated. Two distinct descending patterns were observed: (a) a stable rolling/sliding movement along the slope, and (b) a fluctuating pattern along the slope. Five dimensionless coefficients (Reynolds number (Re), drag coefficient, lift coefficient, moment coefficient and rolling coefficient) were used to analyze the observed processes. The vortex structure in the wake of the sphere gives a lift force to the sphere, which in turn controls the different descending patterns. It is found that the generation of a vortex is not only governed by Re, but also by particle rotation. Relationships between the forces/moments and the dimensionless coefficients are established.


Introduction
This study examines the descending process of a sphere on an inclined slope submerged in a liquid using a micro-scale fluid-solid interaction numerical model that couples the Lattice Boltzmann method (LBM) for the fluid and the discrete element method (DEM) for the sphere. Unlike a sphere settling down in a fluid, the descending process observed in this problem is complex due to the interactions between (a) the hydrodynamic forces between the sphere and the liquid and (b) the contact forces between the sphere and the slope. The numerical model is used to find the underlying fluid-solid interaction mechanisms of a sphere descending on an inclined slope under different conditions (slope angle, fluid viscosity, etc.). This process is fundamental in understanding more complex problems, such as submarine debris flow and saltation of sands in civil and hydraulic engineering. Understanding this simple system may reveal underlying physical mechanism, which could help solve related engineering problems.
Carty [1] was the first to study this problem using experiments. The inclined slope was smooth and thus the friction effect was negligible. By varying the particle Reynold number from 0.02 to 9000, it was found the drag coefficient versus Re curve had a trend similar to a sphere settling in an infinite body of fluid (i.e. without a slope, see below) but the drag coefficient value for a given Re was much larger. Smart et al. [2] measured the translational and rotational velocities of a sphere rolling down an inclined plane and developed a theoretical model including the effect of particle surface roughness. The model predictions were compared with the experimental results, which were conducted in very low par-ticle Reynolds numbers (Re < 2.2 × 10 −4 ). Unlike the previous two studies where the surrounding fluid was static, Prokunin et al. [3] studied a sphere in an upward fluid flow parallel to the inclined slope to induce a shear flow condition. Similar to Carty's results, for a given particle Reynolds number (30 < Re < 3000), they found that the measured drag coefficients were larger than those of a sphere settling in an infinite body of fluid.
More recently, Wardhaugh and Williams [4] conducted experiments, in which spheres of different materials and friction coefficients were placed on an inclined slope with different slope angles immersed in stationary fluids with different viscosities. The experiments produced results with a large range of particle Reynolds number from 15 to 50,000 and the effects of slope angle and friction as well as rotation were also studied. When a sphere rolled down an inclined wall, the relationship between drag coefficient and particle Reynolds number was similar to that obtained by Carty [1]. When particle slipping (at large slope angle) was observed, the measured drag coefficients is smaller than those cases which experience only rolling. This is an important finding, which will be discussed in this paper.
Earlier studies on the sphere-wall interaction problem considered a sphere moving parallel to a wall under gravity but without any contact with a surface. In this case, the contact forces do not come into play. Miyamura et al. [5] performed experiments to examine the wall effects on the steady-state settling velocities of a single solid sphere placed inside triangular and square cylinders and between parallel plates. They introduced a wall correction factor to describe the wall effect on the settling velocity. Takemura and Magnaudet [6] examined the lateral migration of a small spherical buoyant particle in a wall-bounded linear shear flow with Re < 2. An empirical fit capable of predicting the migration velocity for varying distances between the particle and the wall was developed.
Analytical work of a sphere moving parallel to a plane was first conducted by Dean and O'Neill [7] and O'Neill et al. [8]. Goldman et al. [9] derived asymptotic solutions of the Stokes equations for both the translational and rotational motions of a sphere moving parallel to a plane wall in the limit where the gap width tends to become zero. Cox and Hsu [10] and Vasseur and Cox [11] proposed an analytical expression of the lateral inertial motion of a sphere near a wall at low particle Reynolds numbers (Re < 0.3). Zeng et al. [12,13] performed numerical simulations of a rigid sphere translating parallel to a flat wall using the spectral element method. Both lift and drag forces were studied and the relationship between the drag and lift coefficients with the Reynolds number (from 2 to 300) at different distances from the wall were obtained. These studies show the importance of examining the interaction of drag and lift forces to find the underlying mechanism of the sphere-wall problem, which is useful in investigating the problem addressed in this study.
For the descending process of a sphere placed along an inclined slope, most previous studies focused on the kinetic parameters (e.g. velocity, Reynolds number, drag force and drag coefficient). Limited work has been done in terms of the mechanical variables, such as the contact forces and the hydrodynamic forces. In this paper, we, therefore focus on finding the mechanism of interactions between the hydrodynamic forces and the contact forces and their role in determining the pattern of movement of the sphere. Various factors such as the properties of the sphere (density, size), slope inclination (slope angle, friction) and fluid (density, viscosity) can play a role in determining the pattern of movement of the sphere. This study pays special attention to the viscosity of the surrounding fluid, slope angles and the friction coefficient between the sphere and the slope, while the density and the diameter of the sphere are kept constant for brevity.
In this study, a three-dimensional coupled lattice Boltzmann method (LBM) and discrete element method (DEM) code was developed to study the problem of a sphere falling down an inclined slope. The particle scale modelling allows examination of the underlying mechanisms such as velocity, forces and vorticity. The other advantage of LBM-DEM simulation is that the properties of fluid, sphere and wall can be modified independently. The results from the LBM-DEM simulations can provide quantitative information so that the underlying physical mechanism can be captured. An in-depth understanding of this simple system may reveal the underlying physical mechanism, which can be used to solve related engineering problems such as hydroplaning observed in submarine landslides.
This paper is organized as follows. Section 2 describes the coupled LBM-DEM algorithm adopted in this study, and Sect. 3 validates the coupled method with two benchmark problems (sedimentation of a particle in a square cylinder and a sphere rolling down an inclined slope). After validation, the developed method is further applied to study the movement of a sphere descending through a liquid along an inclined plane. The results are summarized and discussed in Sect. 4. The conclusions are presented in Sect. 5.

Numerical model
In this section, a brief description of the lattice Boltzmann method (LBM), the discrete element method (DEM) and the coupling strategy is presented.

D3Q15 LBM
Originated from the lattice gas automata (LGA), LBM solves the continuous Boltzmann method with discrete lattices and converges to the Navier-Stokes equations at low Mach numbers of less than one. LBM has great advantages over traditional numerical methods, resulting in wide applications in hydraulic, civil, chemical and even medical engineering [14]. LBM can simulate various hydrodynamics [19] and offers intrinsic parallelism due to its local nature. The lattice Boltzmann approach has the advantage of accommodating large particle sizes and the interaction between the fluid and the moving particles can be modelled through relatively simple fluid-particle interface treatments [35], which allows for modelling complex boundary conditions.
Because of its simplicity, lattice Boltzmann Bhatnagar-Gross-Krook (LBGK) [15] model is adopted in this study. The evolution equation can be expressed as where f i (x, t) is the density distribution function of a particle moving with velocity c i at position x and time t, τ is the relaxation time and f eq i (x, t) is the equilibrium distribution function of f i (x, t) given by where ρ is the fluid density, u is the fluid velocity, c s is the speed of sound in the lattice, and ω i represents the weights related to the lattice speeds. For the D3Q15 lattice model [16] used in this study, the weighting factors ω i are given as while the 15 discrete velocities are The macroscopic variables can be obtained as where p is the pressure and can be obtained directly from density. The kinematic viscosity υ of the fluid is determined by where c s = c/ √ 3 for D3Q15 model and c = x/ t is the lattice speed. x is the lattice spacing and t is the discrete time step. It is noted that the relaxation time τ is the most important parameter of the LBGK model, which controls the relaxation speed of density distribution functions from nonequilibrium state to equilibrium state and is related to the viscosity of fluid. More details on LBM can be found in [14].
Since turbulence is involved in this problem, large eddy simulation (LES) is incorporated in LBM to account for turbulent effects. Turbulent dynamics is intrinsically nonlinear, the influence of the small removed scales on the large resolved ones must be considered and have been widely adopted within the Navier-Stokes framework. The LES strategy consists of removing the smallest turbulent scales of the flow, whose contribution to global features of the flow is assumed to be small, therefore allowing for the use of a much coarser grid resolution and a significant reduction in the number of degrees of freedom of the computational model. Since nonlinearities in LBM equations differ from those of the Navier-Stokes equations, the use of an eddy-viscosity model can be interpreted as a convenient trick to close the LBM-LES equations [17]. The widely used single-parameter Smagorinsky subgrid model [18] is adopted in this study for numerical stability [19][20][21]. The turbulence is modelled by adding an eddy viscosity term (which models the turbulence) to the true fluid viscosity, which is related to the relaxation time according to Eq. (6). Thus, both relaxation time and viscosity can be split into two parts: fluid part υ, τ and turbulent part υ t , τ t .
The turbulent viscosity is related to turbulent relaxation time by τ t = 3υ t while υ t can be calculated explicitly from the strain rate tensor, which is obtained directly from the second order moments,Q, of the non-equilibrium distribution functions Q i j : where S c is the Smagorinsky constant; S is the characteristic value of the strain rate tensor and τ * = τ + τ t is the total relaxation time. In this way, the turbulent viscosity and relaxation time can be obtained directly from the non-equilibrium distribution functions and are used in the evolution equation at every time step. It can be observed that the lattice Boltzmann equation is not changed and its simplicity is maintained. Further details can be found in [22,23].
In LBM, all the evolutions and calculations are expressed in non-dimensional lattice units. The basic principle of the conversion between LB and SI units is dimensional analysis using fundamental parameters like lattice spacing x, time step t and density ρ. For example, when sphere velocity in lattice unit v l is obtained, then sphere velocity in the physical unit is represented as v p = v l x/ t. More details on the conversion between the LB and SI units can be found in [22,23].

Discrete element method
DEM solves movements of particles, like collisions between particles and wall, congestion and migration, which has been successfully applied in rock mechanics, geotechnical and chemical engineering [24].
The motion of every sphere follows Newton's second law: where m i , I i are the mass and moment of inertia of particle i, V i , ω i are the translational and rotational velocity, respectively while a superposed dot indicates a time derivative, F c i j indicates the contact force of particle i due to particle j and R i is the vector connecting the center of particle i to the location of the contact point. F e i denotes the external force (gravity), while F f i and T f i are the hydrodynamic force and moment exerted by fluid, respectively, and their expressions are presented in the next section. For this study, the soft contact model was applied to calculate the contact forces which are composed of normal force and tangential force, i.e., F c i j = F n n + F t t, where n, t are the unit normal and tangential vectors while F n , F t are the normal force and shear force of particle i due to particle j, which are calculated as follows: in which k n , k t are the normal and tangential contact stiffness while c n , c t are the normal and tangential viscous damping coefficient; δ n , δ t are the normal overlap and tangential displacement while v n , v t are the normal and shear velocity at the contact point, respectively.

Coupling LBM and DEM
To solve this problem, both the contact forces between the sphere and slope, and the hydrodynamic forces exerted by the fluid need to be coupled. In this study, the immersed boundary technique proposed by Noble and Torczynski [25] is used for LBM-DEM coupling because it establishes an accurate and smooth lattice representation of solid particles to reduce the fluctuation of the computed hydrodynamic forces [22]. A local ratio ε s is defined as the area fraction of the nodal cell covered by a solid particle. An additional collision term, s i , which accounts for the effect of solid particle on fluid, is added in the evolution equation of those lattice nodes (fully or partially) covered by a solid particle, where B s is a weighting function for the additional collision operator s i based on ε s .B s , s i are written as follows, where i donate the opposite direction of i, v p is the rigid body velocity of particle node. The hydrodynamic force exerted on one solid particle by liquid can be obtained by summing the momentum change of all the boundary nodes: where n is the number of boundary nodes and Q = 15 for D3Q15 model. x n is the coordinates of the lattice node n while x p is the coordinates of the particle center. The interactions between fluid and particles are bi-directional. The discrete time step also needs to be coupled and the coupling strategy adopted is described in [22,23]. The time step for LBM t is determined by Eq. (6). The time step for DEM t D is calculated as per Eqs. (44) and (45) in [22] and is related to the normal stiffness and mass of the particles. Thus, there are two time steps used in the coupled LBM-DEM, t for the fluid flow and t D for the particles. Since t D is generally smaller than t, t D is reduced to a new value t S so that the ratio between t and t S is an integer n S , n S = t/ t D + 1 and t S = t/n S .
The immersed moving boundary method is used to solve the interactions between fluid and particles and calculate the hydrodynamic forces and torques of the sphere exerted by the surrounding fluid. More accurate results can be obtained when more lattices are used to represent the sphere, leading to an increased computational cost. The error is an artefact of this approach and is not from DEM. It has been proven that the immersed moving boundary method has quadratic convergence with grid refinement and satisfying results can be obtained with at least 5 lattices representing a sphere along its diameter [25]. The specific lattice resolution depends on the required accuracy and computational cost.

Validation of the numerical model
The performance of the proposed coupled LBM-DEM method was validated with two benchmark tests with different interactions and boundary conditions; they are (a) sedimentation of a particle in a square cylinder and (b) a sphere rolling down an inclined slope. Results are compared with the analytical solutions or experimental results from the literature. In general, a good agreement is obtained. In this section, only the case of particle sedimentation in a square cylinder is described while the problem of a sphere rolling down an inclined plane (without fluid) is presented in "Appendix".
In order to validate the interactions between the fluid, the particle, and walls, simulations of a sphere settling in a square cylinder were conducted and compared with the experimental results as shown in Fig. 1. The size of the square cylinder is 0.01 m (L) × 0.01 m (L) × 0.06 m and is kept constant while the diameter of the sphere (D) is varied from 0.0012 to 0.0092 m, similar to the experiments [5]. Non-slip boundary condition was applied to all six walls. The kinematic viscosity of water is υ = 1.0 × 10 −6 m 2 /s and the density is ρ = 1000 kg/m 3 while the density of the sphere is slightly larger, 1000.001 kg/m 3 , to achieve creeping flow. The lattice spacing is x = 0.0002 m. A relaxation time of τ = 0.6 and a time step of t = 0.0013 s was adopted. These three parameters were chosen to achieve a satisfactory result with acceptable computational cost. Lattice spacing determines the number of lattices or node number, whereas the relaxation time is closely related to the numerical stability. Once lattice spacing and relaxation time are chosen, the time step can be calculated with Eq. (6). More details on the parameter selection procedure can be found in [22].
In this study, the sphere which was initially fixed at the center of the cylinder was released to settle under gravity. After some time, the sphere reaches a terminal velocity u t and is compared with the free settling velocity u s calculated from Stokes' law. The Reynolds numbers of free settlings (Re = u s D/υ) were between 0.00094 and 0.42, whereas those of hindered settlings (Re = u t D/υ) by the wall effect were between 0.00067 and 0.014.
The computed dimensionless velocities u t /u s are plotted against different diameter/box size ratios of D/L as shown in Fig. 1. It can be seen the numerical results match well with the experimental data of Miyamura et al. [5]. Due to the drag effect of the surrounding walls, the terminal velocity is Fig. 1 Simulation set-up and results of the settlement of one particle in a square cylinder smaller than that of free settling when the presence of the walls influences the movement of the sphere. As expected, the wall effect increases with an increase in D/L, leading to a smaller dimensionless velocity.
It is noted that there is a small deviation between the two results. One reason is that, when the diameter of the sphere is small, only a few LB lattices interact with the particle, resulting in an inaccurate calculation of the hydrodynamic forces exerted on the particle. The other reason is the existence of the top and bottom walls that are not sufficiently far from the sphere, thus influencing the terminal velocity. The effect of the presence of the top and bottom walls on the deviation is more evident in the case where the ratio of D/L is large. The effect of lattice number of the sphere on the deviation is more evident in the case where D/L is small. These two effects both contribute to the final deviation.
Further simulations were conducted with D = 0.002 m to examine the effect of these two issues. It was found the error between the numerical result and experimental data decreased from 7.7 to 2.9% when the lattice spacing was halved to 0.0001 m and the error decreased from 7.7 to 2.4% when the cylinder length was doubled. Therefore, it can be deduced the errors between the numerical predictions and the experimental results will become smaller when a smaller lattice spacing or a longer cylinder is used.

Overview
For a sphere in a liquid placed on an inclined slope, three sets of forces and moments exist: (i) gravity, (ii) contact forces and moments and (iii) hydrodynamic forces and moments, which  Fig. 2. In this study, the gravity direction is changed to represent the change in the slope angle. The effective gravitational acceleration is g = 9.81(1−ρ f /ρ s ) m/s 2 to account for buoyancy effect. For the simulations, the density of the sphere is ρ s = 1200 kg/m 3 , whereas that of the liquid is ρ f = 1000 kg/m 3 . However, the simulation results are normalized by the gravity force to investigate the behavior of the sphere in general.
In terms of the equilibrium state in the x direction, and in the z direction, while in the y direction, where m, D are the mass and diameter of sphere and θ is the slope angle, while F d , F l , F s , F n , M h , M c are the drag force, lift force, shear force, normal force, hydrodynamic moment and contact moment of the sphere, respectively (see Fig. 2). It is noted that the force balance is different from that of Wardhaugh and Williams [4], who did not take the hydrodynamic moment into account. In their analysis, the normal force is not acting through the center of mass but displaced by a distance b from the center, which contributes a moment M = bF n balanced by an opposing moment from shear force M c = F s D/2. The particle Reynolds number is defined as Re = u x D/υ (where υ is the kinematic viscosity of fluid). For comparison purpose, the drag force and shear force are normalized by mg sin θ , whereas the lift force and normal force are normalized by mg cos θ . The physical time t is normalized as The following four dimensionless coefficients are therefore applied to better understand the descending process, i.e., in which C d , C l , C m , ω are the drag coefficient, lift coefficient, moment coefficient and rolling coefficient, respectively. u x , ω y are the translational and rotational velocity of the sphere. The rolling coefficient is an indication of rolling and sliding of the sphere, with ω = 1 meaning only rolling and ω = 0 only sliding [4]. The diameter of the modelled sphere is 0.02 m. The stiffness of the sphere is k n = 5×10 5 N/m and k t = 1×10 5 N/m while the damping ratios are c n = c t = 0.5. The Smagorinsky constant is S c = 0.1 for turbulent flow modelling according to Feng et al. [22]. Although the model parameter values selected are those used widely [22,23], the effect of those on the results is a subject of future study.
In Sect. 3, the error between numerical results and experimental results decreases with increase in the number of lattices representing the sphere diameter. The error can be neglected when about 15 lattices represent the sphere diameter (0.003 m). Since the flow condition in this section is more complex, 20 lattices over the sphere diameter is chosen to obtain better results with acceptable computational cost. Thus, a lattice spacing of x = 0.001 m was used. The time step t = 0.00033 s was adopted for both LBM and DEM iterations, resulting in a lattice speed 3 m/s to ensure a Mach number of less than one.
The fluid domain was bounded by six stationary walls and the boundary effects were examined prior to the main study. Apart from the bottom wall that represents a slope, the other five walls, especially the upper wall, have an effect on the movement of the sphere as discussed by many researchers [6,10,12,26]. In order to estimate the impact of the upper and side walls on the results, several simulations with different domain sizes were conducted. The slope angle was 75 • , whereas the friction coefficient was 0.2679 (tan 15 • ). The fluid viscosity was 1.0 × 10 −4 m 2 /s. The boundary condition for all walls was no-slip. The sphere was first fixed at   Table 1. When the height of the domain is larger than 0.06 m or the width is larger than 0.1 m, the final parameters of the sphere like Re, velocity, and forces remains almost unchanged. That is, when the domain size was larger than 0.6 m × 0.1 m × 0.06 m, the boundary effect became small and the error due to the presence of the upper and side walls was neglected. Thus, 0.6 m × 0.1 m × 0.06 m was chosen as the domain size for this work to take account of both accuracy and computational cost.
The descending process of the sphere is closely related to Re of the sphere, which is determined by gravity, contact forces and hydrodynamic forces. Different Re of the sphere can be obtained by changing the slope angle, viscosity or the friction between the sphere and slope. Table 2 shows a summary of 62 case studies in which the slope angle, viscosity and the friction between the sphere and slope were varied. In the first series, the effects of slope angle and fluid viscosity on the descending pattern were studied. Three slope angles (45 • , 60 • , 75 • ) were used while the fluid viscosity was changed from 1.0 × 10 −6 to 6.0 × 10 −4 m/s 2 (i.e. the relaxation time τ was modified from 0.501 to 1.1). The friction coefficient, μ, between the sphere and slope was fixed to be 0.2679 (tan 15 • ). In the second series, the effect of slope friction on Re and other parameters was further studied. The friction coefficient was varied from 0 to 11.43 (i.e. the friction angle varies from 0 • to 85 • ). The slope angle was fixed to be 75 • , whereas the viscosity was set to 1.0 × 10 −5 m 2 /s. Following the work of Feng et al. [23], we used the simple and convenient D3Q15 model for LBM and incorporated the Smagorinsky subgrid model for numerical stability. D3Q15 model may break the rotational invariance due to truncation errors at high Reynolds number [30][31][32][33]. D3Q27 model with Multiple-Relaxation Time (MRT) is currently the best model to simulate turbulence at high Reynolds flows. Most of the data used to obtain the quantitative relations in this study are for cases in which Re is less than 300. More efforts will be made in the future to study this topic with D3Q27 model with MRT model.

Descending patterns
For a sphere moving down an inclined slope without fluid, the velocity increases at a constant acceleration while the contact forces remain unchanged [27]. From the fluid coupled simulations conducted in this study, two distinct patterns of particle behavior were observed as shown in Fig. 3: (i) stable and (ii) fluctuating. For the stable case, all the forces do not change with time when the steady state was reached. For the fluctuating case, the forces were changing and oscillating with time.
The drag force and shear force in the x-direction act as a resistance to the movement of the sphere, balancing the driving force caused by the gravity component parallel to the slope. Since the drag force generally increases with Re, the sphere tends to reach a terminal velocity under the actions of gravity, drag force and shear force instead of accelerating as in the no fluid case.
The lift force works together with the normal contact force to balance the gravity component normal to the slope. According to Zeng et al. [12,13] who conducted simulations of a sphere settling down close to a vertical wall, there are three primary contributions to the lift force: (i) shearinduced lift, (ii) rotation-induced lift, and (iii) wall-induced lift. The lift force of the sphere originates from the vorticity Fig. 3 The instantaneous Re and normalized normal force with time for the two descending patterns and pressure differences around the sphere. In this problem, the existence of the bottom wall breaks the symmetry of the pressure and vorticity fields and is the main reason for the lift force. Thus, the wall-induced lift force is dominant and is a result of two fundamental mechanisms: pressure difference effect [4,12] and vorticity effect [12]. Figures 4 and 5 show the instantaneous velocity field and vorticity field around the sphere, respectively. In Figs. 4a and 5a, the results of the plan view at z = 0.01 m (i.e. one particle radius distance from the slope) are presented for the stable case of Re < 100). Both the velocity field and vorticity field are symmetrical on this plane as expected. The side view at y = 0.05 m for the two patterns are shown in Figs. 4b, c and 5b, c. In all patterns, both the velocity field and vorticity field are asymmetrical in this plane. The asymmetry of the vorticity distribution in the wake of the sphere contributes to an asymmetric induced flow. This generates a lift force on the particle directed away from the slope [13].
When Re is smaller than 100, the vortex in the wake behind the sphere is stable, as shown in Figs. 4b and 5b. When Re > 100, one-sided double threaded wake vortices are observed as shown in Figs. 4c and 5c and the effect of vorticity on lifting becomes particularly strong [12]. The vortex swings and even sheds at a certain frequency. The pressure and vorticity distributions on the surface of the sphere are changing and the lift force fluctuates. With an increase in Re, the hydrodynamic forces increase and thus the lift force increases. The normal force then becomes smaller due to the force balance in the z-direction. The motion of the particle is determined by the balance between gravity, hydrodynamic forces, and contact forces. As the particle is always in contact with the wall in the present simulation cases, the normal force is not zero and there is friction.
In Fig. 6, the averaged lift coefficient and lift force of the last 5000 steps of the stable and fluctuating case simula-tions are plotted against Re. Figure 6 shows differences in the lift coefficient with slope angle. The differences are caused by rotation-induced lift, i.e. Magnus lift effect [12,23]. To investigate the magnitude of the Magnus lift coefficient, simulations were conducted by forcing the rotation to be zero (sliding only). The cases with two values of Re (1.07 and 107) were examined.
As shown in Fig. 7, the rolling coefficient increases with decreasing slope angle (i.e. more rolling) but it is zero for the zero rotation cases. The lift coefficients for different slope angles are almost identical at the same Re, when the sphere is not allowed to rotate. Therefore, different rolling coefficients and lift coefficients are obtained for different slope angles with same Re when the sphere can rotate freely.
As observed in the streamline patterns in Fig. 8, the effect of rotation is to drive the vortex away from the sphere, reducing the vorticity driven lift effect. When the slope angle decreases, the contact moment caused by shear force increases due to increase in the normal force and hence the sphere tends to roll. More rotation gives less lift effect (or reduction in the lift coefficient as shown in the right graph in Fig. 7) for a given Re; this is more pronounced in the Re = 107 cases. In the zero rotation cases, the lift coefficient is large, making the particle to fluctuate. In summary, the role of rotation is to reduce the overall lift coefficient. The lift force obtained in our simulations is the combined effect of wall-induced lift and rotation-induced lift.
If there is no rotation of the sphere, the existence of the two regimes is a Reynolds number effect: at higher Re fluid structures becomes unstable leading to the fluctuating behavior. The effects of changing fluid viscosity and slope angle are directly related to the Reynolds number. If the sphere is allowed to rotate, the situation is different. The existence of the two regimes is a combined effect of both Reynolds number and rotation.

Slope angle and fluid viscosity effects
The effect of slope angle and fluid viscosity on the descending behavior of a sphere was investigated. Figure 9 shows the variations of normalized normal force with time for different slope angles and fluid viscosities. As the fluid viscosity decreases, the descending pattern changes from the stable condition to the fluctuating condition (see the cases of slope angle = 75 • ). For the stable cases, lower viscosity means smaller drag force and larger velocity (Re). As a result, the lift force increases while the normal force decreases. For the fluctuating cases, both the velocity and normal force change with time.
To illustrate the effect of fluid viscosity and slope angle more clearly, only results of the 9 cases (3 slope angles, 45 • , 60 • , 75 • , and 3 viscosities, 1.0 × 10 −6 , 1.0 × 10 −5 , 1.0 × 10 −4 m 2 /s) are presented in Fig. 10 for comparison purpose. The temporal variations of Re, rotational velocity, normalized lift force, normalized drag force, normalized nor-  As expected, the sphere velocity (Re) increases with slope angle and decreases with the increase in viscosity. As shown by the Fig. 10e, f, the normal and shear forces decrease (or hydrodynamic and drag forces increase) with time because of the increase in Re with time. The rotational velocity increases to a steady state in the two high viscosity cases (1.0 × 10 −5 and 1.0 × 10 −4 m 2 /s). The rotational velocity decreases with increase in slope angle due to decrease in shear force and moment, which indicates more sliding.

Slope friction effect
In this part, the effect of the friction between the slope and sphere was investigated by changing the friction angle from 0 • to 85 • . The slope angle was fixed to be 75 • , whereas the viscosity was set to 1.0 × 10 −5 m 2 /s.
The instantaneous Re and normalized normal forces of different friction cases are shown in Fig. 11. When the friction angle is less than 60 • , the velocity and normal force become stable with time. When it is greater than 60 • , the normal force starts to fluctuate with time even the Re values are smaller than the cases of friction angle less than 60 • . This is because rotation velocity (or rolling) increases with friction angle, which in turn influences the vortex and the movement of the sphere. Figure 12 shows the velocity and vorticity fields when the friction angle is 0 • , 50 • and 85 • . The vortex structures between the three cases are significantly different. The vortex structures of the friction angle of 0 • and 50 • cases become stable after a normalized time T = 19.8 but the magnitudes of velocity and vorticity are different. The magnitude of the normal force is different between these two cases because of the difference in the vortex structures. The vortex for a friction angle = 85 • case sheds regularly. As a result, the lift force and the normal force dynamically change, as observed in Fig. 11.
The average values of selected variables for the last 5000 steps are plotted against the friction of the slope as shown in Fig. 13. The average Re decreases with increasing friction angle due to an increase in the resisting force, i.e., the shear force. All the variables show a transition in the mechanism when the friction angle is between 30 • and 60 • . As shown by the rolling coefficient plot, when the slope friction angle is greater than 60 • , the sphere movement becomes rolling dominant, whereas it is sliding dominant when the friction angle is less than 30 • .
As discussed earlier, the rotation has a large effect on the vortex structure and hence the lift force. When the sphere behavior becomes rolling dominant, the lift coefficient decreases. The normal and shear forces and moment increase, which reduces Re. But the vortex structure becomes more unstable as shown in the 85 • case in Fig. 12. This, in turn, makes the sphere to fluctuate, as shown by the large friction cases in Fig. 11.
As also discussed in the earlier sections, the simulation results shown in this section highlight the importance of rolling behavior, which is influenced by the friction between the sphere and slope. Enhanced rolling decreases the lift force due to the creation of vortices. This, in turn, reduces the velocity because of increased contact forces. However, rolling also makes the vortices to become unstable, which induces the sphere to exhibit an oscillatory behavior.

Interrelationships of variables
In this section, the averaged results from the last 5000 steps of the stable and fluctuating cases presented in Sect The pattern of the descending process of a sphere is related to the lift coefficient, which is a function of Re as well as the slope angle (or the vortex structure affected by the rolling behavior of the sphere) as shown Fig. 14a. As indicated in [12], for a sphere moving parallel to a wall without contact in a stagnant fluid, the lift coefficient of free rotating condition is the sum of wall-induced lift coefficient without rotation and rotation-induced lift coefficient, that is, C l = C l,W + C l,M , where C l is the lift coefficient for free rotating condition while C l,W , C l,M represent wallinduced lift coefficient and rotation-induced lift coefficient, respectively. When a sphere moving parallel to a vertical wall, the direction of hydrodynamic moment is anticlockwise and thus the sphere rotate anticlockwise. The effect of rotation-induced lift is to increase the total lift coefficient. When the sphere descending along an inclined slope in our simulations, the direction of hydrodynamic moment is anticlockwise but the contact moment due to shear force is clockwise. The contact moment is dominant and make the sphere rotate clockwise. As result, the effect of rotationinduced lift is to decrease the total lift coefficient. Thus, in our cases, C l = C l,W − C l,M . Zeng et al. [13] obtained an empirical relation, 0.313 + 0.812e (−0.125Re 0.77 ) , for the translation-induced lift coefficient when the particle is in contact with a vertical wall but without contact forces. In our case with contact forces, C l,W = 0.41 + 1.2e (−0.17Re 0.64 ) gave a good fit.
For rotation-induced lift coefficient, it is assumed that it is a linear function of rolling coefficient ω [12,28], C l,M = kω. Zeng et al. [12] showed that k is within the range of 0.25 and 0.35, while Bagchi and Balachandar [28] showed that k is around 0.275. It is found that k = 0.3 gives good fit for our simulation data. The total lift coefficient for a sphere descending along an inclined slope in liquid can then be described as follows: The comparison of the simulated data to the values from Eq. (18) and Zeng et al. [13] is shown in Fig. 14a. It is noted that there are errors between the simulated lift coefficient and Eq. (18) when Re > 100 as the vortex structure becomes unstable. Further work is needed to examine unstable vortex structures.
The simulation results show that the drag coefficient depends mainly on Re but is independent of slope angle (or rolling) as shown in Fig. 14b. The following equation is therefore proposed: This is based on the expression of the drag coefficient of a sphere moving down in an infinite fluid domain by Clift et al. [29]. Due to the contact forces, the drag coefficient of a sphere descending along an inclined slope is much larger than that of Clift et al. Unlike the lift coefficient, the relation Fig. 11 The instantaneous Re and normalized normal force of the descending process with different friction angles between drag coefficient and Re was found to be unique for different slope angles. Sphere rotation has little effect on the drag coefficient, which is consistent with the published results [12,28]. The moment coefficient and the rolling coefficient (Fig. 14c, d) are influenced by Re and slope angle.
The results of the averaged hydrodynamic and contact forces and the moments are shown in Fig. 15. As discussed earlier, they are functions of Re as well as slope angle (or rolling). For the normal force, when Re < 100, the best fitting line can be drawn as: while the other forces and moments can be derived by force analysis, which are summarized as follows: in which a(θ ), b(θ )are the coefficients of the fitting function. The curve fitted results are shown in Fig. 15.

Conclusions
The descending process of a sphere along an inclined slope in fluids was numerically studied with a coupled 3D LBM-DEM code. The effects of slope angle, fluid viscosity and friction coefficient between sphere and slope on the descending process were investigated. In this paper, results of both hydrodynamic and contact forces and moments are discussed using five dimensionless coefficients, i.e., Reynolds number, lift coefficient, drag coefficient, moment coefficient 1. The movement of a sphere placed along an inclined slope in a fluid is complex because of the interaction between hydrodynamic forces and contact forces. Numerical simulations show that two patterns exist: stable and fluctuating. In the stable pattern, the normal force of the sphere reaches a steady state. In the fluctuating pattern, the normal force fluctuates with time. 2. The existence of an inclined slope breaks the symmetry of vortex and pressure fields around the sphere and thus leads to a net lift force normal to the slope, which is governed by the vortex in the wake of the sphere. The vortex structure determines the lift force of the sphere and thus governs the changing patterns of the normal force. It was found that the vortex structure is not only determined by translation (Re) but also by sphere rotation. 3. Slope angle and slope-sphere friction influence the movement of a sphere. On one hand, they impact the translational motion by influencing the driving force (mg sin θ) and resisting forces (drag force and shear force). On the other hand, they play important roles on the rotational motion, which in turn creates different vortex structures even though the translation motion (Re) is the same. 4. Enhanced rolling decreases the lift force due to the creation of vortices. This, in turn, reduces the translational velocity because of increased contact forces. However, rolling also makes the vortices to become unstable, which causes the sphere to exhibit an oscillatory behavior. For a given friction angle between slope and sphere, both the dimensional forces and moments and the dimensionless coefficients can be related to Re and rotation velocity (or slope angle). Empirical relationships between them are quantitatively proposed, which show that the moving state of a sphere is predominantly affected by the friction coefficient. Thus, the empirical formulas should be modified to incorporate the effect of friction. That is, the forces and dimensionless coefficients are functions of Re, slope angle and friction coefficient, which is planned for future study.
The key point of this paper is to study the different descending patterns of the sphere and the influencing factors like viscosity and slope angle. Following the work of Feng et al. [23], we used the simple and convenient D3Q15 Fig. 15 The averaged hydrodynamic and contact forces and moments and related fitting lines model for LBM and incorporated the Smagorinsky subgrid model for numerical stability. D3Q27 model is currently the best model to simulate turbulent flow, while D3Q15 model and D3Q19 model may break the rotational invariance due to truncation errors [30][31][32][33]. Such improvement on the LBM model will provide an opportunity to examine the studied behavior at very large Re conditions.
In this study, the density and diameter of the sphere, as well as the density of the fluid are all constant but their effects can also be investigated. Furthermore, the deformation of the sphere is also not considered in the current work and the work of Takemura [34] on the movement of a deformable bubble near a vertical wall can be referred to for further research. where θ, φ, φ critical are slope angle, friction angle and critical friction angle between sliding and rolling cases, respectively; m, g are mass of the sphere and gravitational acceleration while N , S represent normal force and shear force of the sphere; u, ω are the translational and rotational velocity, respectively. In the simulations, the slope angle was fixed to 45 • while the friction angle was changed from 0 • to 80 • . The radius of the disk or sphere was 0.1 m and density was 1200 kg/m 3 . The stiffness values of the ball were k n = 5 × 10 5 N/m and k t = 1 × 10 5 N/m while the damp ratios were c n = c t = 0.5. The time step was 0.0001 s. The particle started moving from the stationary state when simulation began and reached the steady state when both translational and rotational acceleration became constant. The numerical results, as well as analytical solutions are both showed in Fig. 16. For easy comparison, the normal force is normalized by mg cos θ , shear force is normalized by mg while translational accelerationu and rotational accelerationω are normalized by g. Figure 14 shows that the numerical values are consistent with the theoretical solutions. In the sliding regime, the shear force and angular acceleration increase with friction angle while the changing trend of the translational acceleration is opposite. When the friction angle becomes greater than the critical friction angle, the sphere rolls and both the contact forces and accelerations are independent of friction angle. In summary, both sliding and rolling occur in the sliding regime, whereas only rolling exists for the rolling regime.