1. Introduction
Granular flows are known to exhibit a wide variety of features that are often challenging to study and predict. Much of the difficulties arise from inter-particle friction and inelasticity during collisions. At the particle level, each grain is a macroscopic solid, but a large agglomeration of these particles can easily exhibit solid-, liquid- and gas-like states (Jaeger, Nagel & Behringer Reference Jaeger, Nagel and Behringer1996; Goldhirsch Reference Goldhirsch2003; Andreotti, Forterre & Pouliquen Reference Andreotti, Forterre and Pouliquen2013a). A large number of granular flows involve interaction with obstacles such as in mixers, grinders, conveyors, landslides, avalanches, etc. (Wassgren et al. Reference Wassgren, Cordova, Zenit and Karion2003; Delannay et al. Reference Delannay, Valance, Mangeney, Roche and Richard2017). Protective structures erected against landslides and avalanches is just one example where a good understanding of the dynamics of granular flows around obstacles is needed (Faug et al. Reference Faug, Childs, Wyburn and Einav2015). Interaction of granular flows with solid obstacles results in the formation of shock waves across which there exist large gradients of velocity and volume fraction (Amarouchene, Boudet & Kellay Reference Amarouchene, Boudet and Kellay2001; Rericha et al. Reference Rericha, Bizon, Shattuck and Swinney2001). Analogous to gas-dynamic shocks (Anderson Reference Anderson2004), shock waves in granular flows can be attached or detached depending on the shape of the obstacle (Gray, Tai & Noelle Reference Gray, Tai and Noelle2003; Hákonardóttir & Hogg Reference Hákonardóttir and Hogg2005; Cui, Gray & Johannesson Reference Cui, Gray and Johannesson2007; Gray & Cui Reference Gray and Cui2007; Khan et al. Reference Khan, Hankare, Kumar, Kumar, Verma and Prakash2019, Reference Khan, Hankare, Verma, Jaiswal, Kumar and Kumar2022).
Shock waves form at supersonic speeds, which essentially means that the relative velocity of the obstacle with the fluid exceeds the velocity at which finite perturbations can travel. In the case of molecular fluids, disturbances travel via molecular vibrations and the interaction of conservative force fields associated with the molecules because of which the propagation is efficient and fast. Disturbances in granular medium, however, propagate through direct collisions between macroscopic grains, which is relatively a less efficient mechanism. Therefore, sonic velocity in granular flows is much smaller than molecular fluids, such as air. Earlier studies have predicted sonic speed in granular flows to be less than 
$1\,{\rm m}\,{\rm s}^{-1}$, and therefore shock waves in granular flows are easily formed at speeds that are commonly encountered in day-to-day applications (Savage Reference Savage1988; Rericha et al. Reference Rericha, Bizon, Shattuck and Swinney2001; Heil et al. Reference Heil, Rericha, Goldman and Swinney2004; Amarouchene & Kellay Reference Amarouchene and Kellay2006; Khan et al. Reference Khan, Verma, Hankare, Kumar and Kumar2020).
A typical detached granular shock wave, as demonstrated by Amarouchene et al. (Reference Amarouchene, Boudet and Kellay2001), is parabolic in shape and consists of an inner static zone surrounded by an envelope of moving grains such that the velocity increases nonlinearly in the radial direction. Using particle velocities obtained from experiments, Boudet, Amarouchene & Kellay (Reference Boudet, Amarouchene and Kellay2008) demonstrated that the transition from high velocity in the free stream to very low velocity in the heap (downstream of the shock) occurs via multiple granular collisions within the shock wave region. In more recent investigations, it was demonstrated that the particles around the shock wave have two prominent sub-populations; one corresponding to the particles upstream of the shock front having supersonic velocities and the other corresponding to the subsonic particles in the downstream (Boudet et al. Reference Boudet, Amarouchene and Kellay2008; Vilquin, Boudet & Kellay Reference Vilquin, Boudet and Kellay2016; Vilquin, Kellay & Boudet Reference Vilquin, Kellay and Boudet2018; Khan et al. Reference Khan, Hankare, Verma, Jaiswal, Kumar and Kumar2022). They extended the bimodal distribution proposed by (Mott & Harold Reference Mott and Harold1951) and later investigated by other researchers experimentally (Holtz & Muntz Reference Holtz and Muntz1983; Pham-Van-Diep, Erwin & Muntz Reference Pham-Van-Diep, Erwin and Muntz1989; Mazouffre et al. Reference Mazouffre, Vankan, Engeln and Schram2001) and numerically (Bird Reference Bird1970, Reference Bird1978) for molecular gases to account for dissipative effects in granular shock waves. Both continuum and molecular-based approaches can be used to model granular shock waves. A continuum-based approach is particularly suited for large-scale applications such as landslides and avalanches but cannot predict the inner structure of the non-equilibrium shock wave (Gray et al. Reference Gray, Tai and Noelle2003; Cui & Gray Reference Cui and Gray2013). The discrete element method (DEM), however, can expose the inner structure of the shock wave at the microscopic level but is restricted to small domains owing to the computational cost (Padgett, Mazzoleni & Faw Reference Padgett, Mazzoleni and Faw2015). Rericha et al. (Reference Rericha, Bizon, Shattuck and Swinney2001) successfully used continuum and kinetic theory simulations to capture the bulk properties of shock waves around a wedge. Results from both the techniques compared well with the experiments.
Shock waves in granular flows remain a challenging subject owing to the limited measurement techniques and difficulties in theoretical modelling. Although earlier investigations have successfully revealed many of the interesting features of the shock waves in granular medium (Amarouchene et al. Reference Amarouchene, Boudet and Kellay2001; Padgett et al. Reference Padgett, Mazzoleni and Faw2015; Vilquin et al. Reference Vilquin, Boudet and Kellay2016; Karim & Corwin Reference Karim and Corwin2017; Garai, Verma & Kumar Reference Garai, Verma and Kumar2019), there persists a dearth of understanding on some of the very fundamental aspects. The objective of the present work is to probe deeper into the fluid-dynamic aspects of the dynamic zone comprising the shock wave and the dense heap close to the cylinder. While the shock wave is a gas-like zone with strong agitations in particle motion, the grains within the heap are slow-moving owing to granular collapse. The dynamics result in substantial shear and pressure variations within the system that are not given due considerations in the earlier works on the subject. The present work elucidates salient features of shock waves that are formed when a granular stream passes around a circular cylinder using the open-source DEM solver – LAMMPS Improved for General Granular and Granular Heat Transfer Simulations (LIGGGHTS) (Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012). The simulation model is validated by experimental results performed in a rectangular channel confined by two parallel glass plates. A comprehensive analysis is then carried out for a more generic case where a vertically falling granular stream interacts with a stationary solid cylinder, with no influence from the side walls. The complicated fluid-dynamic aspects of the shock wave and the dynamic heap near the cylindrical obstacle are elucidated and are discussed in detail. The origin and the role of the streaming pressure and the shear stress give important insight into the dynamic nature of granular shocks. It is interesting to see how a gas-like stream of particles could transform into a fluid- and then into a solid-like state in such a simple flow situation.
The remainder of the article is structured in the following manner. First, the numerical method is explained in brief in § 2. Then a detailed discussion on the validation is presented where the experimental set-up and preliminary results are compared with the results from DEM simulations in § 3. This is followed by the results and discussion in § 4, where a detailed discussion on the flow field is performed purely on the basis of numerical simulations. Finally, concluding remarks are presented in § 6.
2. Numerical method
In the present simulations, grains are modelled as soft spheres, where inter-particle forces during granular collisions can be modelled using various contact theories. Here, the Hertz-contact theory is used to model elastoplastic behaviour of the forces using springs and dashpots in the normal and the tangential directions (Hertz Reference Hertz1896; Johnson Reference Johnson1987; Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012; Andreotti, Forterre & Pouliquen Reference Andreotti, Forterre and Pouliquen2013b). The normal force, 
$F_n$, is formulated using (2.1) and the tangential force, 
$F_t$, using (2.2) (Di Renzo & Di Maio Reference Di Renzo and Di Maio2004; Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012). Note that these forces are nonlinear because the coefficients vary with the normal displacement according to the Hertz contact theory (Hertz Reference Hertz1896; Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012). The forces act only when the distance 
$r$ between the centres of two particles of radii 
$R_i$ and 
$R_j$ is less than their contact distance 
$d' = R_i + R_j$,
$$\begin{gather} F_n = K_n \,\delta_n - \gamma_n\,{v_n}^{rel}, \end{gather}$$
$$\begin{gather}F_t = \min\displaystyle\left(\left|K_t\delta_t - \gamma_t\,{v_t}^{rel}\right|, \mu_f\,|F_n|\right). \end{gather}$$
Here, 
$\delta _n = d' - r$ is the normal overlap, 
$v_n^{{rel}}$ is the relative velocity between two particles, 
$\delta _t$ is the tangential displacement, which is measured as the relative tangential movement at the surface from the time of contact, and 
${v_t}^{rel}$ is the tangential displacement velocity (Di Renzo & Di Maio Reference Di Renzo and Di Maio2004; Ai et al. Reference Ai, Chen, Rotter and Ooi2011). The normal force has two terms, an elastic-repulsive force with a spring coefficient (
$K_n$) and a damping force with a damping coefficient (
$\gamma _n$). The tangential force also has two components: a varying shear force with a spring coefficient (
$K_t$) and a damping force with a damping coefficient (
$\gamma _t$). The tangential displacement gives rise to the shear force between the particles for the duration of the time they are in contact. Here, 
$\mu _f$ is the coefficient of Coulomb friction, while subscripts 
$n$ and 
$t$ denote the normal and tangential components, respectively. The asterisk (
$*$) denotes equivalent quantities derived in the case of contact between two different spherical particles (Johnson Reference Johnson1987). The parameters 
$K_n$, 
$K_t$, 
$\gamma _n$ and 
$\gamma _t$ are calculated from Young's modulus (
$Y^*$), shear modulus (
$G^*$), Poisson's ratio (
$\nu$) and the coefficient of restitution (COR) (Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012). The maximum tangential force is limited to 
$\mu _f\,|F_n|$, where 
$\mu _f$ is the coefficient of friction (COF) to satisfy the Coulomb friction criterion, as shown in (2.2). Note that the particles studied here are assumed to be perfectly spherical without any rolling resistance. The numerical integrations are carried out using the velocity Verlet scheme of Kruggel-Emden et al. (Reference Kruggel-Emden, Sturm, Wirtz and Scherer2008).
3. Experimental set-up and validation
A new set of experiments are performed in a very controlled environment for the purpose of validating the numerical model. The experimental set-up used to generate a dilute granular stream is shown in figure 1. It consists of an inclined chute made of a smooth glass panel that is 300 mm wide and 900 mm long. The cylindrical disk of diameter 39.8 mm and thickness of 4.7 mm is placed at a distance of 500 mm from the top that acts as an obstruction to the incoming flow. Another glass panel is placed on the top (and parallel to the bottom panel of the chute) such that it touches the top surface of the cylindrical obstacle, and thus maintains a constant channel spacing of 4.7 mm throughout. The granular material used for the experiments is glass beads with a nominal diameter of 1.9 mm and specific gravity 
$\gamma = 1.6$. The hopper at the top is used for feeding grains such that its opening matches perfectly with the channel inlet. The velocity of the incoming stream is changed by adjusting the inclination of the channel, 
$\alpha$. Once the hopper gate is removed, grains slide down the inclined channel such that they are mostly in contact with the bottom surface. A light-emitting diode (LED)-based monochrome light panel illuminates the chute from the bottom. This additional illumination is needed because images are acquired with an exposure of the order of nanoseconds. Monochrome images of the flow field are acquired with an i-Speed 713R high-speed camera at a frame rate of 4000 frames per second with a resolution of 
$2048\times 1536$.
Figure 1. Experimental set-up. The inset shows the coordinate system and the flow direction with respect to the obstacle.
The numerical simulations for validating the flow field are carried out by simulating the glass channel section of the experimental set-up. An insertion surface is placed 
$10$ cm upstream of the cylinder from which particles are inserted at a velocity and mass flow rate matching the experimental values. Table 1 summarises the parameters used in the simulations. A lower Young's modulus than the actual value of the material helps in using a larger time step for the simulation without much sacrifice on the accuracy (Yan et al. Reference Yan, Wilkinson, Stitt and Marigo2015). The Poisson's ratio was assumed to be 
$0.45$. A first guess for the values of coefficient of friction and the coefficient of restitution for glass beads is obtained from the literature (Tang et al. Reference Tang, Song, Dong and Song2019). Fine-tuning of these parameters was then done by running several simulations around the base value. An excellent match between experimental and numerical results is obtained and the final values based on these simulations were then fixed, as shown in table 1.
Table 1. Parameters used for validation simulations.
Validation is based on the comparison of the velocity field around the cylinder obtained from experimental and numerical simulations. A snapshot of the flow field is shown in figure 2(a) from experiments and figure 2(b) from simulations for the channel inclination of 
$70^\circ$. The mass flow rate of the oncoming stream is 0.4 kg per unit channel width per second, and the velocity just ahead of the shock front is 
$2.1\,{\rm m}\,{\rm s}^{-1}$. A sudden clustering of grains near the cylinder nose, wake and the grain-free vacuum region bears a good resemblance in the two cases.
Figure 2. Snapshot of the flow field: (a) from experiments and (b) from DEM simulations. (c) Velocity across the shock wave along the dashed line shown in panel (d). (e,f) Probability distribution function of velocity for a small window around the cylinder shown in panel (d). Mass flow rate and free stream velocity are 
$m = 0.55\,{\rm kg}\,{\rm s}^{-1}$ and 
$V_\infty =1.92\,{\rm m}\,{\rm s}^{-1}$ for 
$\alpha = 40^\circ$ (e), and 
$m = 0.4\,{\rm kg}\,{\rm s}^{-1}$ and 
$V_\infty =2.51\,{\rm m}\,{\rm s}^{-1}$ for 
$\alpha = 70^\circ$ (f).
Velocity from the experimental data is obtained by performing particle image velocimetry (PIV) on a set of 100 instantaneous images acquired through a high-speed camera. An averaged flow field is obtained by averaging the velocity data from all the instantaneous frames. Figure 2(c) compares the velocity across the shock wave just ahead of the cylinder nose (along the dashed line shown in figure 2(d) for reference) for channel inclination, 
$\alpha = 40^\circ$ and 
$70^\circ$ (a zero value of z corresponds to the cylinder surface facing the free stream). A sharp decrease is observed initially inside the shock wave, which is followed by a more gradual decrease downstream of the shock up to the cylinder surface. A sharp decrease in velocity is a typical character of a shock wave. The minimum velocity near the cylinder surface approaches an almost zero value. A good match in the overall velocity profile from the cylinder surface to the free stream is observed for both the channel inclinations.
Figure 3(a) compares the velocity (
$V_z$) distribution of particles within a small square region of size approximately equal to 1.3 times the cylinder diameter, as illustrated by the rectangular box in figure 2(d). A reasonably good match in the velocity distribution indicates the robustness and accuracy of the numerical model in resolving the flow structure at the molecular (or granular) level. This is particularly important because one of the important objectives of the present study is to elucidate the physics of the structure of the shock wave and the dense flow between the shock wave and the cylinder wall. The probability distribution function represents a typical bi-modal structure that was earlier observed for gas-dynamic shocks (Mott & Harold Reference Mott and Harold1951; Pham-Van-Diep et al. Reference Pham-Van-Diep, Erwin and Muntz1989) and granular shocks (Vilquin et al. Reference Vilquin, Boudet and Kellay2016, Reference Vilquin, Kellay and Boudet2018). The velocity on the 
$x-$axis is normalized by the free stream value, so that the peak at 
$V/V_\infty =$ 1 represents the undisturbed supersonic flow. Another peak at 
$V/V_\infty = 0$ corresponds to densely clustered grains around the cylinder surface. The spread in the central region of the distribution corresponds to the particles within the shock wave, or in other words, those particles from the free stream that have just entered into the shock wave region and are undergoing collisions. A detailed discussion on these features is given in the work of Vilquin et al. (Reference Vilquin, Kellay and Boudet2018) and the references therein. A small portion of particles have negative velocity indicating particle motion towards the upstream direction; this is owing to the particles getting reflected from the cylinder surface and or from the dense heap. This also explains our reason for choosing a particular window for making probability distribution functions (p.d.f.s). While the choice of window size and location is arbitrary, it should actually give a good representation of the flow field around the shock wave. For example, if the window is too large in the upstream direction, it will result in overpopulation of free stream particles, and therefore, the distribution of the particles in the shock wave and the inner heap will not be well resolved. Further, a choice of taking grains with a diameter of 1.9 mm is made so that the size of the grain and the channel is comparable. This is because the PIV algorithm assumes the flow field to be two-dimensional. Therefore, if the grain size is too small, there will be many layers of particles inside the heap, which will remain hidden in the still images. In such a case, the p.d.f.s will not give an accurate representation of the dense region of the shock and the heap.
Figure 3. Simulation domain.
The final domain used for the main case study, as discussed in the next section, is different from that used in the validation study. This is because the authors wanted to perform a detailed analysis for a very fundamental and generic case – flow past a circular cylinder (without any complication of sidewalls). However, setting up such a flow in reality is difficult. So, we performed experiments using a closed channel where experiments and diagnostics are relatively easy. Therefore, in § 4, details of the simulation domain are first given before explaining the results.
4. Results and discussion
Once the simulation model is validated, a relatively simple domain (shown in figure 3) is used for detailed analysis. Spherical grains with a diameter 
$d = 800\,\mathrm {\mu }{\rm m}$ and density of 
$2500\,{\rm kg}\,{\rm m}^{-3}$ are released continuously from an insertion surface into the channel. The periodic boundary condition is used in the 
$x$- and 
$y$-directions, whereas grains leaving the 
$z$-direction are deleted. For post-processing, an in-house code is used, where the properties are averaged over bins of 
$1 \times 1\,{\rm mm}^2$ area in the 
$x$–
$z$ plane. The values of Young's modulus and Poisson's ratio used in the simulations are 
$5 \times 10^6$ N m
$^{-2}$ and 
$0.45$, respectively. The values of the cylinder diameter (
$D$), coefficient of restitution and the coefficient of friction are 40 mm, 0.4 and 0.5, respectively, unless stated otherwise.
Particles are inserted at an 
$x$–
$y$ plane placed 0.045 m upstream of the cylinder, as shown in figure 3. The number of particles inserted is estimated in accordance with the predefined mass flow rate. Particles are allowed to fall freely around the cylindrical obstacle resulting in a detached bow shock.
Figure 4(a) shows a sudden change in the volume fraction around the bow shock and a wake in the shadow of the obstacle completely devoid of grains. The volume fraction (
$\phi$) increases from a value of 0.2–0.5 immediately across the shock front but varies slightly thereafter in the dense regime downstream of the shock. A volume fraction of 0.55–0.62 is established close to the cylinder surface, which is close to the random close-packing fraction of spherical particles. A maximum value of 
$\phi = 0.62$ is achieved near the cylinder surface. The flow expands as the grains move downstream of the compression heap, resulting in a gradual decrease in the volume fraction in the wake. The expansion and compression effects are quantified by calculating the divergence of velocity around the cylinder (figure 4b). High magnitudes of divergence near the shock wave points to the compressible nature of the shock wave, whereas, inner shear layer in the wake exhibits high magnitude owing to the expansion effects. The phenomenon is analogous to expansion waves that are observed at the trailing regions of obstacles in gas dynamics (Rericha et al. Reference Rericha, Bizon, Shattuck and Swinney2001; Anderson Reference Anderson2004).
Figure 4. Contours of steady-state flow quantities obtained from DEM simulations for free stream velocity 
$= 1\,{\rm m}\,{\rm s}^{-1}$ and volume fraction, 
$\phi _\infty = 0.22$: (a) volume fraction, 
$\phi$; (b) divergence of velocity, 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {V}^*$; (c) speed, 
$V$ 
$({\rm m}\,{\rm s}^{-1})$; (d) granular temperature, 
$T^*$; (e) streaming pressure, 
$P_s^*$, and (f) total pressure, 
$P_t^*$. Here, 
$X^*$ = X/d, 
$Z^*$ = Z/d, 
$T^*$ = T/
$\sqrt {gd}$ and 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {V}^*$ = 
$\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {V}/(V_\infty d^{-1}$). Note: the value of 
$V$ and 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {V}^*$ is undefined for the grain-free region in the lee side of the cylinder, and therefore, it is artificially set to zero during post-processing for simplicity.
Velocity contours in figure 4(c) reveal a sharp jump across the shock front and a relatively gradual change from the shock to the cylinder surface. After multiple collisions within the shock region, grains lose energy owing to friction and inelasticity that eventually leads to granular collapse near the cylinder walls. The resulting flow structure resembles a dynamic heap with a solid-like dense inner zone surrounded by an envelope of energetic granular fluid. These flow features resemble those reported earlier by Amarouchene et al. (Reference Amarouchene, Boudet and Kellay2001) and Buchholtz & Pöschel (Reference Buchholtz and Pöschel1998). The presence of a dynamic heap is typical to granular shocks with no such correspondence to the molecular-based fluid.
The granular temperature, 
$T$, which is defined as one-third of the variance of the velocity field, is an important ingredient in the application of gas kinetic theory for granular flows. The speed of sound, 
$C_s$, is calculated using the following relation (Savage Reference Savage1988; Vilquin et al. Reference Vilquin, Kellay and Boudet2018):
\begin{gather} C_s= \sqrt{T\beta \left(1+\frac{1}{3}\beta + \frac{\phi}{\beta}\frac{\mathrm{d} \beta}{\mathrm{d} \phi}\right)}, \end{gather}
\begin{gather} \beta = 1+\frac{2(1+e)\phi} {\left(1-\left(\dfrac{\phi}{\phi_{max}}\right)^{(({4}/{3})\phi_{max})}\right)}, \end{gather}
where 
$e$ is the coefficient of restitution and 
$\phi _{max}=0.65$ is the maximum random close packing fraction. The Mach number, 
$M$, is defined as the ratio of the average particle velocity to the local speed of sound. The sonic speed and the Mach number across the shock wave are shown in figure 5. Here, 
$C_s$ is small near the surface of the cylinder owing to the absence of granular motion in the dense heap. The value increases rapidly near the shock wave region owing to high fluctuations in the flow velocity. Upstream of the shock wave, the value of 
$C_s$ again reduces to zero owing to the absence of fluctuation in the free stream velocity. Following the sonic velocity, the Mach number increases from zero in the static heap region to high values away from the cylinder. The free stream value of M approaches exceedingly high values owing to a very low granular temperature. The sonic values obtained in the present study agree well with the estimates reported in the previous studies (Rericha et al. Reference Rericha, Bizon, Shattuck and Swinney2001; Amarouchene & Kellay Reference Amarouchene and Kellay2006). The sonic line that demarcates the subsonic and the supersonic regime is shown in the inset, and as can be anticipated, it depends on the free stream flow conditions.
Figure 5. Quantities calculated from the gas granular theory. Here, 
$R$ is the radius of the cylinder and 
$\delta$ is the heap height (see figure 4f). (a) Speed of sound. (b) Mach number. The contour of the sonic lines is shown in the inset of panel (b).
Granular particles generate stresses as they interact with each other and with obstacles (Savage & Jeffrey Reference Savage and Jeffrey1981). The total Cauchy's stress tensor (
$\sigma$) is generally symmetric for granular flows and is composed of the streaming stress component 
$\sigma _s$ and the collisional stress component 
$\sigma _c$ (Savage & Jeffrey Reference Savage and Jeffrey1981; Savage Reference Savage1988). The streaming component originates owing to the transport of momentum of the particles as they move through the bulk of the material, and the collisional component appears owing to the forces between colliding particles.
For many-body interactions, the total stress tensor can be computed as
\begin{equation} \sigma = \sigma_s + \sigma_c . \end{equation}The streaming stress component is given by (Campbell Reference Campbell1989)
\begin{equation} \sigma_s = \frac{1}{V_b} \sum_{i=1}^{N_b}m~C_i\otimes C_i, \end{equation}
where 
$V_b$ is the volume of the bin, 
$N_b$ is the number of particles inside the bin, 
$m$ is the mass of the particle, 
$C_i$ is the instantaneous velocity of the particle 
$i$ and 
$\otimes$ is the dyadic product. The collisional stress component for the monodisperse medium is given by (Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012).
\begin{equation} \sigma_c = \frac{1}{V_b} \sum_{i=1}^{N_b}\sum_{j=1}^{N_i} \frac{1}{2} (F_{ij}\otimes r_{ij}), \end{equation}
where 
$N_i$ is the number of neighbouring particles that are in contact with particle 
$i$, 
$F_{ij}$ is the force exerted by particle 
$j$ on 
$i$ and 
$r_{ij}=r_i-r_j$ is the relative position of the particle 
$i$ with respect to particle 
$j$. The streaming pressure (
$P_s$) is calculated as the average of the trace of the streaming stress tensor. The contours of normalized streaming pressure 
$P_s^* = P_s/(mgd^{-2})$ are shown in figure 4(e). Because the streaming pressure is generated owing to granular fluctuations, it highlights the regions of high interlayer mixing of momentum. The value of streaming pressure is high near the shock front and the cylinder surface (away from the stagnation region) where granular fluctuations are relatively high.
The total pressure (
$P_t$) is computed as the average of the trace of the total stress tensor. The contours of the normalized total pressure 
$P_t^* = P_t/(mgd^{-2})$ are shown in figure 4(f). It can be observed that the total pressure is significantly higher than the streaming pressure, especially in dense regions. This indicates that the collisional stresses contribute the most to the total pressure in the dense regions.
The acquired stress values have been compared with the stress values obtained through the conservation of momentum equations in the vertical direction (at 
$\theta = 0^{\circ }$ figure 4b) given by
\begin{equation} \rho\left(\frac{\partial V_z}{\partial t}+V_x\frac{\partial V_z}{\partial x}+V_z\frac{\partial V_z}{\partial z}\right)+\rho g = \frac{\partial\sigma_{xz}}{\partial x} +\frac{\partial\sigma_{zz}}{\partial z}. \end{equation}The right-hand side of the above equation has been calculated directly from the stress field obtained using (4.3) and compared with the results obtained by evaluating the left-hand side of the above equation (see figure 6). A good match, except in the shock wave region owing to the non-equilibrium effects, is achieved. The stress calculations made by the in-house post-processing code show good agreement with the continuum fields obtained from theoretical equations.
Figure 6. Comparison of the sum of differential stresses 
$G={\partial \sigma _{xz}}/{\partial x} +{\partial \sigma _{zz}}/{\partial z}$ obtained directly from stress calculations and from the momentum conservation equation.
5. Anatomy of granular heap
The detailed physics of the inner structure of the heap is studied through the variation of flow properties across the heap in different radial directions normal to the surface of the cylinder (see figure 4b for the schematic of radial lines). Radial lines 
$r^*=\sqrt {x^2+z^2}/d$ are drawn at angular positions 
$\theta = 0^{\circ }, 15^{\circ }, 30^{\circ }, 45^{\circ }, 60^{\circ }$ and 
$65^{\circ }$. The dashed line at 
$r^* = 25$ denotes the location of the cylinder wall and the right-dashed line at 
$r^* = 48$ marks the location of the compression front (for the 
$\theta = 0^{\circ }$ case).
Reading figure 7(a) from right to left, the volume fraction undergoes a sharp rise through the shock wave, attains a plateau and then decreases suddenly near the cylinder wall. The presence of a similar plateau region for all radial lines indicates that the volume fraction remains nearly the same inside the heap.
Figure 7. Variation of properties with radial distance from the centre of the cylinder. The dashed lines encompass the dense region for the 
$\theta = 0^{\circ }$ case: (a) volume fraction 
$\phi$; (b) normal (
$V_n$) and tangential (
$V_t$) velocities; (c) granular temperature; (d) shear stress; (e) streaming pressure and (f) total pressure (absolute value).
The velocities transformed in the normal/radial (
$V_{n}$) and the tangential (
$V_{t}$) directions are shown in figure 7(b). Normal velocity shows a piecewise linear profile with two sections: one corresponding to the shock wave region with a higher slope and the other to the inner heap with a smaller slope for all the radial lines. There is a sharp decrease in the magnitude of 
$V_{n}$ around 
$r^* = 48$ owing to the presence of the shock wave, but unlike the volume fraction, the velocity continues to decrease within the heap (the region between the shock wave and the cylinder wall). It is interesting to note that there is a slight overshoot of 
$V_{n}$ for 
$\theta = 60^\circ$ and 
$65^\circ$ near 
$r^* = 25$, which suggests the presence of flow separation at the cylinder walls. The tangential velocity 
$V_{t}$ (see figure 7b) remains zero along the 
$\theta = 0^\circ$ line, which is anticipated owing to the lateral symmetry of the flow field. For other values of 
$\theta$, 
$V_{t}$ decreases slightly within the shock wave and the downstream region, and finally drops to zero near the cylinder wall.
Figure 7(c) shows that the value of the granular temperature in the free stream is almost zero and rises rapidly in the shock wave region owing to an increase in the fluctuations in velocity. The energy is dissipated within the shock wave owing to the high rate of collisions, and consequently, cooling takes place inside the heap. A slight increase in the temperature is again observed near the wall region, prominently for high values of 
$\theta$, owing to the fluctuating components generated in the velocity as the grains slide down the cylinder surface.
Figure 7(d) shows the variation of the normalized shear stress 
$\tau ^*_{x'z'}=\tau _{x'z'}/(mgd^{-2})$ along radial lines. It remains zero along the 
$\theta$ = 
$0$ line owing to the lateral symmetry of the flow field. A sharp increase in the shear stress is observed near to the cylinder for theta values (
$\theta = 15^\circ - 45^\circ$) that is attributed to the strong shearing action provided by the grains sliding atop the cylinder. The shearing is evident from the high gradients in the tangential velocity, as discussed earlier (see figure 7b). Moving further along the radial direction, the value of shear stress drops to a minimum in the central region and then again rises owing to the shearing action of the free stream near the shock front. For higher values of 
$\theta$ (
$\theta = 60^\circ$ and 
$65^\circ$), the spike in the shear stress at the surface is absent. Interestingly, the value of shear stress near the surface becomes zero somewhere between 
$\theta = 60^\circ$ and 
$65^\circ$ locations. At low angles, the pressure near the surface is higher (see figure 4f), which means that the grains exert larger forces on each other and on to the surface. This results in more sustained contacts, which eventually give rise to high Coulombic friction and hence high shear stresses. However, for high values of 
$\theta$, grains tend to detach from the cylinder walls resulting in weak and short-lived inter-particle contacts. The position of zero shear stress appears to be the location where the flow starts to expand. Analogous to fluid dynamics, the point of zero shear stress in the present case may be crudely regarded as the point of flow separation. As we move away from the cylinder walls, the shear stress for 
$\theta = 60^\circ$ and 
$65^\circ$ increases gradually. This gradual increase is attributed to the shearing action provided by the direct interaction of the free stream near the shock front.
As shown in figure 7(e), the behaviour of streaming pressure correlates well with the granular temperature (figure 7c). They both reach maxima inside the shock wave, where the normal velocity gradient is the highest, and reduce to low values in the central region. The maximum reached is significantly higher for lower angles. A slight increase is observed near the surface owing to high tangential velocity gradients, as can be deduced from figure 7(b).
Figure 7(f) shows the variation of total pressure inside the heap and the shock wave region. The total pressure demonstrates a piecewise linear variation with a steep rise in the shock wave followed by either a rise or a fall in the downstream region of the heap. The initial rise in the pressure arises from the sudden deceleration of particles owing to multiple collisions within a thin region of the shock wave. The second branch of the pressure has a contribution from the deceleration as well as the weight of the particles acting on top of each other analogous to the hydrodynamic pressure. For the 
$\theta = 0$ line, the pressure rise in the second branch is steepest because the free stream is oriented normal to the cylinder surface that results in maximum deceleration of grains. As the value of 
$\theta$ is increased (
$\theta = 30^\circ$ for example), the component of the velocity of the incoming grains normal to the cylinder decreases, which results in the net increase in total pressure owing to granular deceleration being relatively less. Additionally, for higher theta values, the grains pick up tangential momentum, which causes them to tend to lose contact with the surface resulting in lower stresses near the surface while exhibiting an increase in stresses towards the shock wave. This particular state of the heap, where grains are densely packed but are in motion, exhibits a viscoplastic behaviour. In such a scenario, the volume fraction remains constant and the pressure distribution varies across the flow cross-section. A particularly interesting case is observed for 
$\theta = 45^\circ$ where the pressure in the second branch is almost constant, thus exhibiting a linear plug-like flow profile. For high 
$\theta$ lines (
$\theta = 60^\circ$ and 
$65^\circ$), the grains leaving the shock wave region do not get compressed into a heap, rather flow down into the wake, and thus the pressure decreases.
5.1. Heap height
Because higher pressures give rise to increased frictional forces, as per the Coulomb friction model, the heap can sustain a large number of particles as long as the pressure within the heap is high. The pressure from the impact of the grains near the shock front is most likely transmitted through the network of force chains established through the contact point of neighbouring grains (Radjai et al. Reference Radjai, Jean, Moreau and Roux1996). Therefore, the size of the heap correlates with the total pressure rise around the cylinder surface. In the absence of a continuous free stream flow, the pressure at the shock front reduces to zero and the heap breaks down layer by layer.
The heap height varies almost exponentially with the volume fraction, 
$\phi$, for constant free stream speeds, as shown in figure 8(a). A denser flow will increase the pressure near the cylinder surface, giving rise to more frictional forces. The increase in frictional forces aids in the formation of a larger heap. The speed of the flow has relatively less impact on the heap height; though an increase in speed of the flow with constant volume fraction would effectively increase the mass flow rate, the particles are more energized within the heap and the heap formed is less stable. The inset image on the bottom right of figure 8(a) shows the variation of heap height with different coefficient of friction values. The higher values of COF allow the particles to transmit more frictional forces from the surface of the cylinder, thus stabilizing the heap further; this results in an increased heap height. The graph shows an asymptotic variation because fewer particles reach the maximum Coulomb friction criteria as the COF value is increased. Figure 8(b) shows the variation of the heap height with 
$D/d$ for two values of 
$\phi$. For dotted curves, 
$D/d$ is varied by changing the cylinder diameter 
$D$ from 0.02 m to 0.16 m, whereas for solid curves, 
$D/d$ is varied by changing the grain diameter 
$d$ from 0.8 mm to 1.6 mm. For higher 
$D/d$ ratios, the non-dimensionalized heap height almost remains a constant, signifying a linear relationship between 
$S$ and 
$D$. Additionally, a change in grain diameter at high 
$D/d$ ratios does not change the non-dimensionalized heap height, signifying that it does not have any effect on the heap height.
Figure 8. Variation of heap height with (a) volume fraction, 
$\phi$, and (b) cylinder to grain diameter ratio, 
$D/d$. The ratio 
$D/d$ is varied either by changing 
$D$ or 
$d$ at a time, where the constant value is denoted by the arrow in panel (b) for that particular branch of the plot.
6. Conclusions
In this work, we studied the formation of shock waves and the inner structure of dynamic granular heaps using the discrete element method for the flow of a dry granular stream past a circular cylinder. A two-dimensional granular stream is generated by allowing grains to fall vertically over a cylinder. The shock wave consists of a thin non-equilibrium zone followed by a dense region near the cylinder surface where grains are almost stationary. The region within the shock wave is dominated by granular collisions resulting in a sharp increase in the granular temperature and the streaming stress component. Strong dissipation in the kinetic energy owing to granular collisions results in a sharp transit from the supersonic to subsonic regime across the shock wave, and the formation of a granular heap near the cylinder surface. The heap height increases nearly exponentially with volume fraction but remains almost insensitive to the free stream flow velocity. Additionally, the height of the heap is found to be correlated with the total pressure values, and is sustained as long as the grains within the heap continue to receive energy from the incoming stream. The results indicate that the formation and the stability of the heap is a consequence of frictional and inelastic dissipation during granular collisions.
Acknowledgements
Authors acknowledge the use of the computing facility at the High Performance Computing Facility of the Computer Center, IIT Kanpur, for carrying out the computations associated with this work. Authors thank Dr N. Tiwari for providing the high-speed camera for experiments.
Funding
This research was supported by Science and Engineering Research Board (SERB) through the Grant Numbers CRG/2019/003989 and CRG/2020/000504, and by the Indian Space Research organization through the Grant Number STC/AE/2019438.
Declaration of interests
The authors report no conflict of interest.
