Granular flow around a cylindrical obstacle in an inclined chute

Granular flows are among the most commonly encountered processes in industry and nature, but understanding the granular media behavior still remains a very challenging task.^1–3 It is understood that a process involved with granular media can be characterized by solid, liquid, or gas regimes;⁴ however, the transition between different regimes often depends on a complex set of conditions associated with both microstructural and macroscopic behaviors of granular materials.^5,6 While a friction criterion is valid for describing the transition between solid and liquid regimes in a sand pile, it is also subject to conditions arising from the initiation and deformation of the flow^7,8 to the size of the system in a rotating drum.^3,9 The transition between liquid and gas regimes, on the other hand, may depend on a combination of macroscopic parameters (e.g., volume fraction) with microscopic parameters such as the coefficient of restitution.^1,10–13

Following the continuum modeling of granular materials by Savage and Hutter,¹⁴ simulations based on depth-averaged granular models (e.g., Ref. 15) allow us to gain a deep insight into the bulk flow behavior of granular materials at fast-moving regions.^16–19 With the further understanding of the constitutive and rheological laws of granular materials,^5,20–25 continuum simulations of granular flows continue to provide reliable solutions to more complex problems. In recent work by Tregaskis et al.,²⁶ for example, they successfully simulated the granular flow at both subcritical and supercritical speeds on a rough slope using a continuum model. In the study of a granular flow around multiple triangular wedges by Khan et al.,²⁷ however, the “dynamic instability” induced in their shock–shock interactions suggests that even for a fast-moving granular avalanche, a single set of parametric conditions may not be sufficient for modeling this more complex phenomenon.

With the introduction of the discrete element method (DEM) by Cundall and Strack,²⁸ DEM-based simulations have become increasingly important for its direct understanding of the granular media behavior incorporating microstructural details.^2,23,29–34 Moreover, the coupling of the computational fluid dynamics (CFD) calculation of the fluid phase with the DEM computation of the granular phase has offered an enhanced platform to understand mixed solid–liquid–gas behavior under more realistic conditions such as granular materials immersed in an ambient fluid.^13,35,36 For the collapse of a granular column, for example, its run-out may behave very differently, depending on whether the use of coarse or fine grains is immersed partially in water.^37,38 On the other hand, the coupled CFD-DEM simulation has made it possible to adopt a much larger time step in computation through the “unresolved” or “resolved” approach,³⁹ where a resolved time step is based on a length-scale smaller than the particle diameter, while an unresolved step can use a much bigger value based on a “parcel,” i.e., a cloud of particles sharing similar properties.

In this paper, we will investigate both the microscopic behavior of particle's using DEM simulation and the bulk flow behavior of granular material using a depth-averaged granular model (hence called the “continuum simulation”). The particle flow will be assessed during its propagation around a cylindrical obstacle down an incline and will be compared with a chute flow experiment of 100's and 1000's granular sprinkles (or nonpareils, of a spherical shape and a nominal diameter of 1 mm) of the same conditions. Figure 1 shows an example of such granular material impacting a cylindrical obstacle, where three different regimes of the flow, solid (quasi-static), liquid (dense flow), and gas (dilute flow),^1,2 are captured during the flowing process. In Sec. II, we will explain the setup of the chute flow experiment, and the determination of the key parameters and conditions used for the DEM and continuum simulations. In Sec. III, we will introduce a method for obtaining the coefficient of friction for particle–boundary and inter-particle interactions by matching the analytical solution of the continuum model with the flow velocity measured by particle image velocimetry (PIV). In Sec. IV, we will investigate the details of the granular flow past the obstacle from onset to run-up and then to the steady state, where the solutions of shock wave and granular vacuum will be compared and analyzed between experimental data, DEM and continuum simulations, and the granular shock theory. A further comparison is made in Sec. VI for the run-out of the flow, of which time-dependent patterns of the flow will be discussed. To show a full development of the granular flow impacting the obstacle, relevant video and animations have been attached as integral multimedia in Figs. 7, 29, and 30.

The set-up of the experiment is shown in Fig. 2, where the chute has a downslope length of 1000 mm and a cross-slope width of 300 mm, inclined at an angle ζ to the horizontal. The sidewall of the chute has a normal height of 100 mm formed from a transparent plexiglass to enable a high-speed camera to capture the motion of granular particles along the sidewall. The bed of the chute is a flat plywood surface, and a hopper is placed at the top of the chute to allow the granular material, the “100's and 1000's sprinkles,” to flow down the chute with the release of a gate, controlled to a defined height.

Also shown in this figure is the notation of the Cartesian coordinate system, O–xyz, with the origin O being at the top left-hand corner of the chute, the x-axis following the downslope direction, the y-axis along the cross-slope direction, and the z-axis pointing upward in the normal direction. We can, thus, define the flow velocity as $u=(u,v,w)$⁠, where $u,v, and w$ are the components in the x, y, and z directions, respectively. A cylindrical obstacle, which is made of transparent acrylic tube with an outer diameter of 28 mm and a height of 100 mm, can be placed on the bed. For this paper's work, the chute is inclined at an angle $ζ=40°$ with a gate height of 40 mm, and the center of the cylinder is placed at x = 300 mm in the middle of the chute, where the flow captured in a downslope range of $x∈[150,480]$ mm will be the focus of the study.

The DEM simulation of chute flow is conducted through a coupled CFD-DEM method (e.g., Refs. 40–42) on a STAR-CCM+ platform, where the discrete element method is used to model the granular particle system, and the volume-averaged Navier–Stokes equation is used to solve the fluid flow through a multiphase-Lagrangian framework; the modeling details are described in  Appendixes A and  B.

In our DEM simulation, the computational domain is chosen to cover the same downslope distance of 330 mm as used in the experiment, but to save computational time, we adopt a lateral width of 200 mm for the cross-slope dimension. This is because the cylinder only has a diameter of 28 mm, such a cross-slope dimension is sufficiently wide to treat the sidewall as a far field boundary without affecting the flow around the cylinder. Correspondingly, the inlet condition is defined at x = 150 mm, which is not the same as the gate condition obtained at x = 0 mm. This treatment has the advantage of simplifying the simulation while avoiding the complexity of the gate flow that can often be subcritical. For convenience, the key physical and control parameters used for the DEM simulation are summarized in Table I.

Since the simulation is operated through a coupling between the CFD and DEM calculation, the control of the time step in the simulation is based on an interactive assessment to both fluid phase and solid phase, in which the time step of the Lagrangian solid phase is based on the size of parcels in an unresolved approach. As shown in Table I, a time step, $Δt=0.001$ s, is used in the implicit unsteady solver in our simulation. For the Lagrangian multiphase calculation, a local time step $δtp$ is further calculated by allowing the Courant number C_o to vary, say, from 0.05 to 0.35 here, which suggests a local time step to be in a range of

where $Δx$ represents a characteristic length scale of the cell containing the parcel. For example, if assuming that $Δx$ is five times the particle diameter, i.e., $Δx=5$ mm, and the maximum speed of the flow field is $max(|V|,|Vp|)=2$ m/s, we may get $δtp$ in a range of $1.25×10−4$–$8.75×10−4$ s. To further improve the accuracy of the DEM calculation, a sub stepping, say with a maximum of $2×104$ steps, is used to the Lagrangian phase, suggesting that a local time step can be as small as around $10−8$ s when resolving the particle details at the microstructural level.

The governing equations used for the continuum simulation are based on the depth-averaged model for free-surface granular flows,^18,19,26 where the flow velocity has been averaged over the flow thickness in the normal z direction by assuming the shallowness of the flow. Therefore, the normal component of velocity w is not considered, and to reflect the depth-averaged treatment, we use notation, $u¯=(u¯,v¯)$⁠, to denote the velocity and its components in the x and y directions, respectively. The continuity equation and the momentum equation can be then given as

where $·$ and $⊗$ are the dot product and dyadic product, respectively, of the gradient operator $∇=(∂/∂x,∂/∂y)$⁠. The source term $S=(Sx,Sy)$ on the right-hand side takes into account of the effect of the gravitational force and the frictional resistance exerted on the basal surface and is given by

where μ is the coefficient of friction, $|u¯|=(u¯2+v¯2)12$⁠, and $i$ is the unit vector in the $x−$ direction.

With the use of the conservative variables h, $m=hu¯$⁠, and $n=hv¯$⁠, we can re-write the system of Eqs. (2) and (3) into a non-strict hyperbolic form

where $U=(h,m,n)T$ with the superscript T denoting the transpose to a row vector, and the respective fluxes and source term vector are

For the finite difference computation of the flow around an obstacle, a body fitted coordinate system,⁴³ namely, $ξ=ξ(x,y)$ and $η=η(x,y)$ as shown in Fig. 3, is required to represent the physical generation of a computational grid around the obstacle. Accordingly, the governing equation (5) can be converted to

with $τ=t$⁠. The transformed variable, fluxes, and source term are

where the unbracketed subscripts denote differentiation with respect to that subscript and the Jacobian coefficient $J=ξxηy−ξyηx=(xξyη−xηyξ)−1$⁠.

The governing equation (7) is solved numerically in an O-grid, where the orthogonality of transformation from the physical domain O – xy to the computational domain $O′−ξη$ is maintained. As shown in Fig. 3, since the study of the flow is focused on a rectangular region of 330 (downslope) by 200 (cross-slope) mm for both experiment and simulation, we set the radius of the outer boundary of the O-grid to R = 245 mm around a circular obstacle of radius r = 28 mm. In the continuum simulation, a shock capturing method based on the non-oscillatory central (NOC) scheme⁴⁴ is used for improving the shock resolution with the use of a “min-mod” limiter.⁴⁵ 

The set-up of the boundary conditions is as follows. A uniform inlet condition of $h0=15$ mm, $u¯0=1.152 1$ m/s, and $v¯0=0$ m/s is given to the upstream region for $x≤x1=150$ mm when t > 0 s. For the part of the outer boundary where $150<x<300$ mm, at the same time, a one-dimensional flow without obstacle is calculated when t > 0 to provide a time-dependent solution of h(x, t) and $u¯(x,t)$⁠, following the propagation of the flow. By doing this, the flow upstream of the shock can be kept as one-dimensional with the use of an O-type grid. For the downstream part of the outer boundary, an outflow condition is applied. To the surface wall of the circular cylinder, the slip condition, $n·u¯=0$⁠, is used, where $n$ represents a unit vector normal to the wall. Moreover, a symmetrical condition can be used to the central streamline if a half of the O-grid is used for simulation. As the numerical method is explicit, the time step is automatically given according to the CFL (Courant–Friedrichs–Lewy) stability condition.^43,45,46

The velocity field of the granular flow along the chute is obtained by the following method. In the experiment, the 100's and 1000's sprinkles were released from the hopper and flowed down the chute, and a series of high-speed video clips, at a time interval of 1 ms (i.e., 1000 fps), were taken from different viewing angles to capture as much detail as possible for the post-experiment examination of the flow. The high-speed footage was then processed by a PIVLab package,^47,48 where the experimental images taken from both overhead and sideview angles were analyzed to achieve a more correlated solution.

As shown in Fig. 4, a sideview image of the flow, over a downslope distance from 70 mm to 400 mm, has been processed and velocities calculated, where a mean velocity field is obtained over 100 frames of video to minimize any randomness from the experiment. We then used a similar method to process the overhead footage taken for the same flowing condition, where velocity shows a robust consistency and repeatability. For the flow upstream of the shock, the lateral velocity component, v, is at least two orders of magnitude smaller than the downslope velocity u and, thus, can be neglected. The measurement of the overhead velocity provides a free-surface distribution of u with the downslope distance, which can also be used to calibrate the velocity field obtained from the sideview measurement.

Shown in Fig. 5 are such velocities obtained accordingly, where the overhead measurement is for a downslope range of 0 to 250 mm (circular symbols), and the sideview measurement covers a range of 70 to 400 mm along the downslope. To obtain a free-surface velocity from the sideview measurement, e.g., in Fig. 4, a series of markers were first made to denote the thickness of the free-surface of the flow at various downslope locations (top panel), which were then used to extract the velocity distribution accordingly from the PIVLab data (bottom panel). By comparing this velocity distribution (square symbols) with the overhead result, we can reach a matched distribution over a total range of 0–400 mm in the downslope direction, where the values for $x∈[70,250]$ mm agree precisely between the two measurements. With such a free-surface velocity as a benchmark, we can further obtain the velocity distribution at the basal surface (diamond symbols), which can be used to determine the coefficient of friction for the particle–wall contact. Also shown in Fig. 5 is an averaged velocity distribution based on the free-surface and basal surface results, denoted in delta symbols. It is equivalent to the depth-averaged velocity used in the continuum simulation and can also be used to calculate the height of a granular shock (see Sec. IV C).

With the velocity profile obtained along the sidewall, it is possible to determine the coefficient of friction based on the continuity equation (2) and the momentum equation (3) by assuming that the flow is steady-state and one-dimensional. Therefore, we can let the cross-slope velocity $v¯=0$ and simplify Eqs. (2) and (3) into¹⁸ 

where S₁ is given by

with the angle, $β=ζ−δ$⁠, measuring the difference between the chute inclination and the basal angle of friction. Let $u¯0$ be the flow velocity, h₀ is the flow thickness at the inlet boundary, and the continuity equation (9) can be integrated directly to give a relation that

The momentum balance equation (10) can be simplified by expanding out the derivatives, using Eq. (9) and then dividing through by h to give

Using (12) to substitute for the flow thickness, Eq. (13) can be integrated to give a cubic equation for the flow velocity,

which can be numerically solved for $u¯=u¯(x)$⁠, of which only the positive root of the solution is valid.

For given inlet velocity $u¯0$⁠, flow thickness h₀, and inclination angle ζ, the matching of the calculated velocity using (14) with the measured velocity can be conducted by adjusting the angle of friction δ until the $u¯∼x$ distributions agree completely. Since the cylindrical obstacle is placed at x = 300 mm down the slope from the gate of the chute here, we only performed the fitting for a downslope range of $150≤x≤400$ mm for the requirements of our DEM and continuum simulations.

With such an approach, we can obtain the corresponding values of δ for the velocities at the basal surface, the free-surface, and the averaged depth of the flow, as shown in Fig. 6. For the velocity at the basal surface, an angle of $δ=21.5°$ is obtained and shall be used to model the friction between the particle and basal surface in our DEM simulation, and, hence, $μ=tan 21.5°=0.393 9$⁠, as previously shown in Table I. With the fitting for the averaged velocity, an angle of $δ=18°$ is obtained and can be used in the continuum simulation of Eq. (7). Similarly, an angle of $δ=13°$ is achieved for matching with the free-surface velocity. If assuming that the flow is of a “laminar” form within its thickness, we may say that the dynamic friction angle of inter-particles can vary from $13°$ from the free-surface layer to $21.5°$ at the bottom layer. From our static experiment, however, we found that the static inter-particle friction angle for the 100's and 1000's is equal to $26°$⁠, and we shall discuss the effect of this angle later in the simulation.

Finally, since the flow is initiated by the inlet condition defined at a downslope distance of x = 150 mm, we get the inlet thickness of the flow $h0=15$ mm, and the depth-averaged inlet velocity $u¯0=1.152 1$ m/s according to the profile given in Fig. 6. On the other hand, the mass flow rate used in the DEM simulation is calculated as follows:

where W denotes the cross-slope width of the flow and ρ_s is the particle density, as shown in Table I. In our case, W = 200 mm and $ρs=811.54$ kg/m³, we, hence, have an inlet mass flow rate $ṁ0=2.804 9$ kg/s.

As the experiment and simulation were conducted independently, it is necessary to calibrate both operations to a comparable timescale. To do this, we first matched the high-speed footage with the simulation result to a temporal instant at which the front of the granular flow is about to impact the leading edge of the cylinder. As shown in the top panel of Fig. 7 (Multimedia view), this temporal node corresponds to a numerical time of 0.088 s. For the result in the second panel for t = 0.188 s, therefore, we then use the same time interval, 0.1 s, to extract the experiment snapshot or simulation data. Since the time interval of the high-speed footage is 1 ms, we shall use the next $100th$ snapshot in the experiment to compare with the simulation solution at t = 0.188 s.

For convenience of discussion, we introduce an average particle velocity,

where $Vp,i$ is the resultant velocity of particle i and N_p is the total number of the particles when a DEM simulation approaches a steady-state, for example, $Np=458 300$ at 1.0 s in the current case. Accordingly, a time-dependent history of $Vp,a$⁠, as shown in Fig. 8, can be used to follow the development of the flow from run-up to approaching a steady state and then to the run-out upon the stop of the injection of the inlet flow. Also shown in the figure is the history developed through the continuum simulation, where its $Vp,a$ is averaged by the total number of gridpoints. It is seen that the run-up stage of the continuum simulation is much shorter than that of the DEM simulation and the experiment. A one-dimension solution for the flow without cylinder is also shown in the figure, where a uniform rectangular mesh, with 166 × 101 gridpoints, is applied to the continuum simulation for a region of 330 × 200 mm under the same h₀ and $u¯0$⁠.

It is seen from Fig. 8 that the primary stage of the run-up is up to a time of around 0.3 s for the experiment and DEM simulation. For the results at 0.188 and 0.274 s shown in Fig. 7, there are a few distinctive features worth noting. First, upon impacting the cylinder, the flow thickness becomes much higher than that of a later developed shock wave, for example, at 0.388 s. This could be of practical importance, for example, to the design for a snow avalanche defense as it may provide a more accurate approximation to the height of an effective defensive infrastructure.^49,51,52 Second, a particle-free region, i.e., granular vacuum, starts to develop at the downstream of the cylinder. Following the shock wave, a significant amount of granular material flies freely alongside the shock tail, indicating a strong mixing of liquid (dense) and gas (dilute) behaviors in the flow.^1,2 The result at 0.388 s is close to a steady-state formation where both the bow shock and the granular vacuum have developed sufficiently. Such a formation can be clearly observed in the experiment, but the corresponding simulation solution at this point seems to take a longer time to settle to a similar steadily state (e.g., as shown in the Multimedia view and in Fig. 28 in the  Appendix).

The shock waves developed in both DEM and continuum simulations are shown in the sideview solutions in Fig. 9. For the DEM solutions on the left-hand side, it is seen that an early impact of the flow against the obstacle, namely, at t = 0.188 s, can produce a much higher jump of the flow thickness up from 80 to 95 mm. Figure 10 shows the corresponding development of the flow thickness for these temporal nodes, where the flow upstream of the shock (i.e., jump) has fully developed at 0.274 s, but the downstream flow undergoes particle-level changes throughout the later times. For the continuum solutions on the right-hand side in Fig. 9, however, the shock wave is shown to be sufficiently developed at t = 0.188 s, which remains nearly identical to other later-time solutions, resulting in a faster run-up rate in the simulation (e.g., as shown in Fig. 8). Another important difference between the DEM and continuum simulations is that the DEM solution has a more extended shock region to the downstream, while the shock region is largely confined near the obstacle in the continuum simulation. This leads to a marked difference for the formation of the granular vacuum.

Corresponding to the t = 0.188 s solution, as shown in Fig. 11, the granular vacuum region is shown widely open as the front of the flow sweeps to the downstream. In the DEM solution, a “low speed region” emerges on the lee-side of the granular vacuum and continues to extend downstream with the increase in time, for example, from t = 0.274 to 0.388 s. Compared with the solutions at 0.688 and 1.0 s, we can see that the shape of the granular vacuum undergoes interesting changes including contraction from 0.188 to 0.388 s and expansion from 0.388 to 1 s, where the particle velocity on the vacuum's lee side accelerates (further details can be seen in the attached animation). This phenomenon may be interpreted as a “hysteresis” between modeling the microscopic dynamics of the particles and exhibiting the bulk flow behavior at the macroscopic level, which becomes more obvious when simulating a shock wave, across which the velocity of the granular particles has been greatly reduced. However, such a hysteresis effect is neither difficult to observe in the experiment nor identified in the continuum simulation where the vacuum develops continuously with its lee-side having the fastest velocity.

For the DEM solutions of 0.188–1.0 s in Fig. 11, a region of “upstream influence” is observed at the front of the shock wave, within which the velocity of the flow becomes substantially lower than the rest of the upstream flow. By extracting the DEM data along the central streamline of y = 100 mm, we can obtain an u–w velocity field following the flow thickness. From Fig. 12, we can see that the velocity starts to decrease around x = 270 mm at the free surface, but the velocity at the basal surface undergoes a more dramatic change toward the leading edge of the obstacle. For the region confined between the free surface, the basal surface, and the forward side of the obstacle, a “channelized” flow is, thus, formed by pushing the particles to a maximum height at the front of the obstacle. The appearance of this upstream influence region is an interesting result captured in the DEM simulation because, usually, a shock wave only imposes a downstream influence through its characteristics. With the increase in the roughness of basal surface, such an influence can become stronger and, hence, extended further to the upstream.²⁶ The time-dependent development of the flow velocity is also given in Fig. 13, where the velocity at the free surface experiences a more dramatic change across the shock, but the velocity at the basal surface takes a longer time to be fully developed.

Figure 14 compares the flow thickness obtained in the experiment, DEM simulation, and continuum simulation. With comparison to the experimental shock height of 72 ± 2 mm, it is seen that the continuum solution overshoots the shock with a height of 82 mm, while the DEM simulation undershoots the shock with a height of only 62 mm. Furthermore, we can calculate the shock condition along the central streamline using the normal shock theory.⁴⁹ By ignoring the effect of the source term in the momentum equation (10), and with the continuity equation (9), we can obtain the jump condition across a normal shock by

where the subscripts “1” and “2” represent the value on the forward and rearward sides of the shock, respectively. By substituting $u¯2$ in (18) through (17), we can have the following jump condition for the flow thicknesses:

where $Fr1=|u¯1|/gh1 cos ζ$ is the local Froude number upstream of the shock. As shown in Table II, the shock height h₂ calculated according to (19) is in good agreement with but smaller than the simulation result, which is mainly caused by the omission of the source term. It is noticed that the flow thickness h₁ obtained in the DEM simulation is significantly smaller than that of the experiment and continuum simulation, hence making h₂ much smaller too. Note also that the upstream velocity $u¯1$ has used an averaged value according to the velocity profile given in Fig. 15, where the DEM's velocity is lower than that obtained in other methods. In this figure, since the velocity distributions obtained from the PIVLab (in square and diamond symbols) are based on a simple one-dimensional flow without the obstacle being considered (see Fig. 6), their trend follows a monotonic increase with the downslope distance, hence showing a clear difference with other velocities obtained from the DEM and continuum simulations for the flow around the cylinder.

Shown in the last row of Table II is the shock standoff distance (SSD) formed between the shock wave and the forward side of the obstacle, which is normalized by the cylinder diameter. We can see that there is an excellent agreement of SSD between the experiment and the continuum simulation, which also agrees extremely well with the analytical solution of SSD.⁵⁰ The SSD solution from the DEM simulation, however, shows a significantly smaller value. If comparing between Figs. 14 and 17 (a solution discussed later in Sec. V A), we can see that the SSD solution can be improved substantially with the enhanced local flow condition upstream of the shock wave in the DEM simulation.

The boundaries of the granular vacuum are compared in Fig. 16 by projecting the numerical solutions of the flow thickness into the x–y plane, where the continuum solution is shown in a flood-type contour, overlapped by the DEM solution in dashed dotted lines. The results are based on the 1 s solutions, where the DEM thickness of the flow has been averaged in the z-direction to form a 2D map in the x–y plane. The boundary of the granular vacuum obtained in the experiment is marked with a solid line with error bars, showing a better agreement with the DEM solution.

In summary, compared with the experimental result, the DEM simulation produces a smaller shock height h₂ due to an underestimated thickness h₁ upstream of the shock, but forms a better granular vacuum downstream of the shock. On the other hand, the continuum modeling produces a more realistic flow thickness upstream of the shock, giving a better overall agreement for the shock height, but there is small discrepancy for the granular vacuum. Moreover, the normal shock relation (19) shows an excellent agreement with the experiment for the shock height, at an error of about $4%$⁠.

From the shock solutions shown in Figs. 14 and 15, it is seen that under the same inlet condition, the flow thickness and velocity still develop differently in the experiment and in the simulation. Between the experimental and continuum results, although the shock height overshoots in the simulation, the flow thickness and velocity (depth-averaged) upstream of the shock are in good agreement. For the DEM simulation, on the other hand, its flow thickness and velocity show a different but more detailed result due to its modeling complexity. We now select three heights, 15, 21, and 25 mm, for the inlet thickness h₀ and discuss how it affects the shock and granular vacuum solutions in the DEM simulation. To keep the inlet condition comparable, we let the velocities at the free-surface and basal surface remain unchanged by prescribing a linear distribution over the depth of h₀. Therefore, the inlet mass flow rate $ṁ0=2.804 9$ kg/s if $h0=15$ mm (previously used); $ṁ0=3.923 1$ kg/s if $h0=21$ mm; $ṁ0=4.598 5$ kg/s if $h0=25$ mm. Also, the friction angle arising from the particle–wall interaction is equal to $21.5°$⁠, and that from the inter-particle interactions is equal to $26°$⁠.

The effect of h₀ on the flow thickness is shown in Fig. 17, where the shock height based on $h0=25$ mm agrees very well with the experiment. From the thickness distributions between $h0=15$ and $h0=25$⁠, we can see that the thickness difference of the flow changes from 10 mm at the inlet (x = 150 mm) to about 4 mm at the upstream of the shock (x = 270 mm). A “transition layer” may be, thus, marked to denote such a change of the free-surface from a parabolic to a linear form, but such a transition is only observed in the DEM simulation. Moreover, the corresponding velocity at the free surface, as seen in Fig. 18, shows a clear drop around this transition layer. At the basal surface, the velocity becomes smaller too than the experimental result. The change of the flow thickness and velocity affects not only the shock formation but also the granular vacuum. In Fig. 19, we map the DEM solutions of 1 s to the x–y plane to show the flow formation around the cylinder. For the boundaries of granular vacuum, i.e., along the h = 1 mm lines, we can see that the boundaries of $h0=15$ and $h0=21$ agree well since their basal velocities remain almost identical. The vacuum boundary formed for $h0=25$ mm is, however, significantly smaller due to its basal velocity being smaller in the vacuum region. [A further comparison of the flow formation for different h₀ is also given in Fig. 29 (Multimedia view) in the  Appendix].

For simplicity, we continue to follow the form, $μ=tan δ$⁠, to define the coefficient of friction and, thus, use δ_wp to denote the friction angle formed between the basal surface and particles, and δ_pp the friction angle arising from inter-particle interactions. The benchmark solution of DEM is based on the condition that $h0=25$ mm, $δwp=21.5°$⁠, and $δpp=26°$⁠.

The effect of δ_wp on shock height, flow velocity, and granular vacuum is compared in Fig. 20, where δ_wp is used $21.5°, 15°$⁠, and $10°$⁠. For the shock height in Fig. 20(a), the solutions of $15°$ and $10°$ are slightly smaller. This is an interesting result since the basal flow velocities for these two angles, as shown in Fig. 20(b), are greater than that of the $21.5°$ by about 0.14 m/s (or about 14%) at the upstream of the shock. If comparing the flow thickness upstream of the shock, say, around x = 270 mm, we can see that the flow thickness of $15°$ and $10°$ is smaller than that of the $21.5°$ by about 1.5 mm (or about $10%$⁠). This may suggest that the flow thickness upstream of the shock can be more dominant in affecting the shock height for a certain range of wall-particle friction angles if the free-surface velocity remains largely unchanged [Fig. 20(b)], which is consistent with the shallowness assumption for the depth-averaged continuum model of granular flow where the inter-particle interactions are not considered. However, the reduction of the basal velocity does make the granular vacuum smaller at $δwp=21.5°$⁠, as shown in Fig. 20(c).

The effect of the static friction coefficient arising from inter-particle interactions is shown in Fig. 21, where three angles of δ_pp, $45°, 26°$⁠, and $13°$ are compared. It is seen that its effect on both shock height and granular vacuum is not as marked as that of δ_wp but the shock region becomes more extended downstream with the decrease in δ_pp [see Figs. 21(c) and 29 in the  Appendix]. In other words, with the increase in the inter-particle friction, particles tend to become more compact and, hence, behave more “liquid-like” when moving around an obstacle on a slope.

One advantage of simulating the granular flow using the continuum model is its accuracy and simplicity of calculating the flow thickness upstream of the shock, e.g., as shown in Fig. 14 and Table II. It is, however, dependent on the correct setup of the conditions including h₀, $u¯0$⁠, and the basal friction coefficient $μ(I)$⁠, where the inertial number I^5,24 represents the ratio of a macroscopic deformation time $(1/γ̇)$⁠, with $γ̇$ being the shear rate, to an inertia timescale $(d2ρs/P)0.5$⁠, with ρ_s being the particle density, P is the normal hydrostatic pressure, and d is the particle diameter. For free-surface granular flows, a bulk inertial number⁵³ has been given in a simpler form:¹⁹ 

which is directly related to the Froude number, particle diameter, and flow thickness. To improve the calculation for the granular vacuum, we now propose a $μ(Ib)$ rheology in a form

where the parameters are given as $μ1=0.3939, μ2=0$⁠, $Ib1=1.0, and Ib2=1.1$ for the present example.

The relation among I_b, $Fr$⁠, and h is shown in Fig. 22, where two types of μ-relation, one with a simple definition $μ=tan (21.5°)=0.3939$ shown in solid line, and another with $μ1=0.3939,μ2=0, Ib1=1.0$⁠, and $Ib2=1.1$ in dashed dotted line, are compared. It is seen that a bulk inertial number $Ib=1.0$ [Fig. 22(b)] corresponds to a Froude number of 7 [Fig. 22(a)], above which the flow regime becomes largely dilute illustrating strong gas behavior with the near-zero basal frictional effect, and we, thus, let $μ2=0$⁠. It is also found from preliminary tests that the transition from a dense flow regime, i.e., for $Ib≤Ib1$⁠, to a dilute regime for $Ib>Ib2$ can take place within a narrow window, for example, a range of I_b from 1.0 to 1.1 can be used here.

With the use of the $μ(Ib)$ relation (21), it is seen in Fig. 22(c) that the granular vacuum has improved significantly with comparison to the simple use of $μ=0.393 9$ where the vacuum closes at the downstream end. Figure 23 shows a further comparison for a larger basal friction $μ=tan (26°)=0.487 7$⁠, where the reattachment point of the vacuum moves further upstream, but the use of the $μ(Ib)$ rheology shows a clear improvement. In this figure, a more abrupt condition with a step change for μ from 0.4877 to 0 at $Ib=0.6$ is also tested and denoted in dashed line, where we can see that the boundary of the vacuum becomes widely open [Fig. 23(c)]. While the use of different basal friction coefficient μ has affected the boundary of granular vacuum greatly, its effect on the shock height is not so obvious, as shown in Fig. 24, which is similar to that which is observed in the DEM simulation, see Fig. 20(a). This is because the flow thickness and velocity upstream of the shock remain approximately the same with a moderate change of δ, e.g., from $18°$ to $26°$ here. Also, because the $μ(Ib)$ rheology has been defined to improve the basal friction condition in the vicinity of the granular vacuum, its impact is localized and does not affect the shock formation either.

The run-out of the granular flow is important to understand how the flow behaves and the granular material deposits upon the stop of the incoming flux from the inlet. To initiate a run-out in our DEM and continuum simulations, we stopped injecting particles after t > 1.0 s by giving zero mass flow rate and zero velocity at the inlet of the domain. To keep it consistent with the discussion for the run-up and steady-state solutions in Sec. IV A, the conditions used for the DEM and continuum simulations remain unchanged (see Table I). Since the inlet flow thickness h₀ still uses 15 mm in the DEM simulation, its shock height is smaller than the experimental result, but the granular vacuum keeps continue to agree well between experiment and simulation. According to the evolution history of $Vp,a$ in Fig. 25, we can see that the time scales developed in the DEM simulation for different conditions discussed in Sec. V remain nearly identical particularly for the run-up and run-out stages and are comparable to that achieved in the one-dimensional continuum simulation where the cylinder is not considered.

The comparison between experiment and simulation for the run-out stage adopts a similar approach as used in Sec. IV A, where four temporal instants, 1.065, 1.120, 1.182, and 1.251 s, are used to capture typical patterns of run-out [see Fig. 30 (Multimedia view) in the  Appendix]. For the continuum simulation, however, the timescale is shorter for a similar solution. It is seen from Fig. 26 that at the instant when the upstream flow is deposited at the front of the cylinder, the solution of the DEM (and thus the experiment) corresponds to t = 1.120 s, but for the continuum simulation, a similar solution can be achieved at 1.102 s. This is consistent with the evolution of the average flow velocity shown in Fig. 8 where the continuum simulation runs out at a faster speed than the DEM simulation.

The sideview solutions obtained in the DEM and continuum simulations are compared in Fig. 27, where we can see that during the primary stage of the run-out, namely, for the time up to 1.182 s, the shock region and granular vacuum remain mostly “uninfluenced” even after the rest of the flow runs much further to the downstream. This observation suggests that the height of a granular shock and the boundary of a granular vacuum could be “marked” more permanently upon an avalanche impacting an obstacle, for example, in a snow avalanche event, where the field observation of the aftermath could offer useful information to re-calculate the avalanching condition to a greater accuracy.^{49,51,52,54,55} On the other hand, the “shielding” of the shock wave to the granular vacuum could ensure the safety zone formed behind an obstacle to sustain over the entire course of an avalanche. For the continuum solutions shown in Figs. 26 and 27, the shock structure and granular vacuum are also very well maintained during the primary stage of the run-out, but the shock height drops much more quickly.

The appearance of the shock wave during the granular flow around a cylindrical obstacle shows a mixing of solid, liquid, and gas behaviors of the granular material, as early shown in Fig. 1, of which details can be captured by the coupled CFD-DEM simulation and high-speed camera. The solid regime is focused at the front of the obstacle where the shock region becomes totally stationary due to its dramatically increased thickness and near-zero velocity of the flow, while the liquid and gas regimes are mixed within the shock extension region to the downstream with the propagation of the flow. The continuum simulation, on the other hand, provides a simpler and faster solution to the shock solution particularly at the front of the obstacle, but it has a much shorter extension to the downstream. This is because that the depth-averaged granular equations (2) and (3) are yet to take the effect of gaseous phase into account in their present form.

The granular vacuum formed at the downstream of the obstacle shows another mixed liquid and gas behavior near the boundary of the vacuum as the flow only has a thin thickness while moving at high velocity. In this circumstance, the local Froude number can usually reach between 7 to 10, hence yielding a bulk inertial number I_b of above 1.0. It suggests that the flow near the vacuum boundary is of strong gaseous behavior, and the local coefficient of friction between particles and basal surface, thus, becomes extremely small, i.e., can quickly reduce to zero, as modeled in the $μ(Ib)$ relation of (21). Therefore, understanding the behavior of granular material during its transportation process is of ultimate and prerequisite importance for modeling and simulating the resulted flowing event in a realistic and appropriate manner. While the friction coefficient $μ(I)$ may reach to a higher critical value when a granular material deposits at the boundary of granular vacuum on a rough surface,^1,20,26 such frictional effect may become negligibly small when granular particles move almost freely close to the vacuum boundary on a smooth surface, and, hence, the use of $μ(Ib)=0$ provides a better solution to the vacuum formation, as shown in our present simulation of the continuum model.

The comparison of the granular shock and vacuum solutions shows a good overall agreement between experiment and simulation at both microscopic and macroscopic levels. Further understanding of the dynamic properties of granular particles in a broad parameter space will provide more insight to simulate the bulk flow behavior of granular materials under complex conditions, where examples of study may include the hysteresis effect in the granular vacuum formation, the upstream influence of the flow at the front of the shock wave (Fig. 11), and the transition layer occurred in the flow thickness (Fig. 17).

Above all, the transition among solid, liquid, and gas behaviors may have never been a simple dependence on a single set of parametric conditions, and, thus, a unified approach of combining the modeling and simulation of granular flows at both microscopic and macroscopic levels may offer a more comprehensive and consistent approach to this challenging task at a higher dimension.

The authors would like to acknowledge the support of the Innovate UK under Grant No. BB/S020993, the EPSRC COVID-19 Grant Extension Allocation (Grant No. EP/V521140/1), and the Royal Society under Grant No. SIF/R2/212009 for conducting this research. Also, the authors are grateful to the valuable discussion on rice milling processes with Mr. Alec Anderson from Koolmill Systems Ltd. of UK, and to the support from Siemens and MayaHTT on the exploration of STAR-CCM+ platform including the DEM component. The team is deeply thankful to colleagues from Hallam on the use of computing facilities and high speed camera. Dr. Xinjun Cui is a Royal Society Short Industry Fellow.

The authors have no conflicts to disclose.

Xinjun Cui: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (lead); Methodology (lead); Supervision (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Matthew Harris: Data curation (equal); Investigation (equal); Software (equal); Visualization (equal). Martin Howarth: Funding acquisition (lead); Investigation (supporting); Project administration (supporting); Resources (lead); Software (supporting); Supervision (equal); Writing – review & editing (supporting). Daisy Zealey: Data curation (equal); Investigation (equal); Software (equal). Reegan Brown: Data curation (supporting); Investigation (supporting); Resources (equal); Validation (supporting). Jonny Shepherd: Data curation (supporting); Project administration (supporting); Visualization (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request. A list of experiment video and simulation animations as an integral multimedia: Figures 7 and 30, video showing the experiment of 100's and 1000's sprinkles impacting a cylindrical obstacle in a chute inclined at $40°$ to horizontal. Figures 7 and 30, animation-1 showing the simulation corresponding to the experiment, taken from a sideview angle. Figures 7 and 30, animation-2 showing the simulation corresponding to the experiment, taken from an overhead angle. Figure 29, animation-3 showing the simulation for an inlet thickness $h0=25$ mm (⁠$δwp=21.5°, δpp=26°$⁠), taken from a sideview angle.

The CFD model for the fluid flow adopts the volume-averaged Navier–Stokes equations, where the volume occupied by the fluid within each cell depends of the volume being taken by the solid particles. Let α_f be the volume fraction, $uf$ is the fluid velocity, and ρ_f is the fluid density, the continuity equation of the CFD model can be given as (e.g., Refs. 39, 41, and 42)

and the momentum equation is

where $∇p$ represents the pressure gradient, τ_f is the fluid viscous stress tensor, $Fpf$ represents the term associated with momentum transfer between the fluid and solid phase, $g$ is the vector of gravitational acceleration, and “$·$” is the dot product and $⊗$ the dyadic product.

The DEM model adopted in our work is based on a soft-sphere approach, where the translational and angular motions of each individual particle can be given in the following form (e.g., Refs. 41, 42, and 58):

where $ui$ and $ωi$ are the translational and angular velocities of particle i, respectively, with m_i being its mass and I_i is the moment of inertia tensor. This model allows particle i to interact with a total number of n_i neighboring particles through both normal and tangential contacts, denoted by subscripts “n” and “t,” respectively. Likewise, the normal and tangential damping forces are denoted as $Fnd$ and $Ftd$⁠, correspondingly. $Tt$ is the torque produced by the tangential force, and $Tr$ is the torque produced by the rolling friction. The gravitational force of particle i is modeled by $g=(g cos ζ,0,−g sin ζ)$⁠, with g = 9.80 m/s², to consider the effect of the chute inclination ζ according to the coordinate system used in the DEM simulation, as shown in Fig. 2.

The contact forces in our DEM simulation are based on the Hertz–Mindlin model,^56,57 where the magnitude of the normal forces including the damping force is given by

In the tangential direction, the magnitude of the tangential forces (including the damping force) is given by

if

otherwise

Other notations in the above expressions are explained as follows. E_eq is the equivalent Young's modulus, G_eq is the equivalent shear modulus, and R_eq is the equivalent radius, between particles A and B. d_n and d_t are the overlaps, and v_n and v_t are the relative velocity components in the normal and tangential directions, respectively, at the contact point. C_fs denotes a coefficient of static friction, and $Nn,damp$ and $Nt,damp$ are the normal and tangential coefficients of damping, respectively. In the translational equation of (B1), there is an additional term $Fo$ to account the coupling between the CFD and DEM calculations, which may include the effect of drag, pressure gradient, virtual mass, moving reference frames, user-defined force, etc., and we here primarily consider the effect of drag and pressure gradient.

For the calculation of the angular momentum equation (B2), the torques are given by

where μ_i denotes the coefficient of rolling friction, $Ri$ is the position vector from the particle centroid to the contact point, and $ωp$ is the component of the particle angular velocity parallel to the contact plane.

Further figures are shown here to give a more detailed picture to the relevant sections in the main text.