2 Simulation implementation

Simulations were implemented using the discrete element method for soft particles in two dimensions, which allows long-lasting multiple contacts between the particles. The original code was developed by Wassgren (Reference Wassgren1996) and was adapted to include periodic boundary conditions in the vertical dimension. Our numerical code considers only the contact forces and gravity. For the normal component of the contact forces, we used a model proposed by Walton & Braun (Reference Walton and Braun1986). In this model, the normal contact force during a contact on a particle $i$ due to another particle $j$ , $\boldsymbol{F}_{N}^{ij}$ , is modelled using a linear spring during loading, and another linear, but stiffer, linear spring during unloading:

(2.1) $$\begin{eqnarray}\boldsymbol{F}_{N}^{ij}=\left\{\begin{array}{@{}ll@{}}-k_{l}{\it\delta}\hat{\boldsymbol{n}} & \text{loading}\\ -k_{u}({\it\delta}-a_{0})\hat{\boldsymbol{n}} & \text{unloading},\end{array}\right.\end{eqnarray}$$

where $k_{l}$ ( $k_{u}$ ) is the loading (unloading) spring constant, ${\it\delta}$ is the overlap between particles, $\hat{\boldsymbol{n}}$ is the unit vector normal to the contact and $a_{0}$ is the overlap at which the unloading force is zero due to the plastic deformation of the particles. For each pair of particles in a specific contact, $a_{0}={\it\delta}_{max}(k_{u}-k_{l})/k_{u}$ , where ${\it\delta}_{max}$ is the maximum overlap for that contact: whenever ${\it\delta}\leqslant a_{0}$ the contact finishes. With these considerations, the coefficient of restitution is

(2.2) $$\begin{eqnarray}{\it\varepsilon}=\sqrt{\frac{k_{l}}{k_{u}}}.\end{eqnarray}$$

The tangential force, $\boldsymbol{F}_{S}^{ij}$ , is modelled as a Coulombic sliding friction element in series with a linear tangential spring (Cundall & Strack Reference Cundall and Strack1979):

(2.3) $$\begin{eqnarray}\boldsymbol{F}_{S}^{ij}=-\min ({\it\mu}|\boldsymbol{F}_{N}^{ij}|,k_{t}|{\bf\gamma}|)\frac{{\bf\gamma}}{|{\bf\gamma}|},\end{eqnarray}$$

where ${\it\mu}$ is the coefficient of sliding friction, $k_{t}$ is the tangential spring constant, ${\bf\gamma}={\it\gamma}{\hat{s}}$ is the total tangential displacement and ${\hat{s}}=-n_{y}\hat{\imath }+n_{x}\hat{\jmath }$ . The value of the tangential displacement is restricted to be smaller than the maximum value ${\it\gamma}_{\text{max}}={\it\mu}|\boldsymbol{F}_{N}^{ij}|/k_{t}$ . The stiffness of the tangential spring is determined as a fraction of the normal loading spring. In our case, $k_{t}/k_{l}=0.5$ .

Table 1. Coefficients of friction ${\it\mu}$ and restitution ${\it\varepsilon}$ , and stiffness of the springs used in the simulations.

The parameters of our simulations are shown in table 1 and they are the same as those used in a previous work by our group (Solano-Altamirano et al. Reference Solano-Altamirano, Caballero-Robledo, Pacheco-Vázquez, Kamphorst and Ruiz-Suárez2013). In choosing these values we assumed that the macroscopic behaviour of the system was not sensitive to the precise value of the stiffness constants because that is what Wassgren (Reference Wassgren1996) found while simulating glass beads. However, we are now working with extremely light particles and such an assumption may no longer hold. It has been shown that when the normal contact stiffness is sufficiently large compared to a characteristic pressure of the system, the simulation is in the asymptotic limit of rigid grains and the macroscopic behaviour is insensitive to the stiffness value (Kneib et al. Reference Kneib, Faug, Dufour and Naaim2016). The simulation is in the limit of rigid grains whenever the dimensionless number $N_{0}=k_{n}/P_{0}d\geqslant 10^{4}$ , where $P_{0}$ is a characteristic pressure of the system. In our case $N_{0}\sim 6\times 10^{3}$ , which means that we are not in the limit of rigid grains, but we are not far from it. Although it will be interesting to study the effect of varying $k_{n}$ on the macroscopic behaviour of our system, we believe that the current values of our parameters are reasonable since in our previous work we were able to quantitatively match our simulations with experimental results (Solano-Altamirano et al. Reference Solano-Altamirano, Caballero-Robledo, Pacheco-Vázquez, Kamphorst and Ruiz-Suárez2013).

Figure 1. (a) Schematic representation of the simulation cell: a dense granular flow confined by lateral walls with periodic boundary conditions in the vertical direction. Small walls at the bottom of the cell control flow velocity and erase memory. (b) Intruders are static disks of diameter $D=2.5~\text{cm}$ , separated by a distance $X_{sep}$ . The medium is composed of grains of diameter $d=0.52~\text{cm}$ . (c) Snapshot of the region around the intruders during flow. Positive lift and drag forces are defined, as are the ‘in’ and ‘out’ measurement regions.

3 Simulation set-up

The system consists of a channel, delimited by a pair of vertical walls separated by distance $L$ , and two static disks of diameter $D=2.5~\text{cm}$ , located side-by-side in the middle of the channel and separated from each other by a distance $X_{sep}$ (see figure 1). We placed two staggered rows of small walls at the bottom of the channel, whose sizes can be varied, to allow for the control of the flow velocity $U_{\infty }$ and to erase the flow memory. The second row of the small walls is useful to further break the memory of the granular medium. From these control walls and the periodic boundary conditions, long-lasting simulations and a homogenous velocity field along the channel are possible (see figure 1 in the supplementary material available at http://dx.doi.org/10.1017/jfm.2016.384).

In each simulation, we created 8000 grains with random initial radii, velocities and positions, which are left to settle under the effect of gravity. Initially, the grains are only a fraction of their final size to avoid overlap, and they grow with time until they reach an average diameter of $d=0.52~\text{cm}$ , with a size dispersion of 6 %. We simulated the grain material as polystyrene, with a mass density of ${\it\rho}_{g}=0.014~\text{g}~\text{cm}^{-3}$ (for the sake of simplicity, we determined the density and mass of the grains as if the grains were spheres of diameter $d$ , although the simulation is performed in two dimensions). The angle of repose of the simulated granular material was measured to be ${\it\theta}=24.5^{\circ }$ , which corresponds to an effective friction coefficient of ${\it\mu}_{eff}=\tan ({\it\theta})=0.45$ .

A wall, located at the bottom of the channel to contain the bed, is removed at time $t=0~\text{s}$ to begin the flow. The height of the column above the intruders stabilises at approximately $h\sim 33~\text{cm}$ (see supplementary movies 1–6). Although $h$ may vary, mainly because of the size of the voids behind the intruders, the variations are small and do not account for the change in drag and lift forces, which were found to depend linearly on $h$ . We measured the lift force $F_{L}$ and drag $F_{D}$ on the intruders, whose positive directions are depicted in figure 1(c) for the right intruder. The mirror image of these forces correspond to the values for the left intruder. We systematically varied the flow velocity, $0.28\leqslant U_{\infty }\leqslant 1.06~\text{m}~\text{s}^{-1}$ , and the separation between intruders, $1.25\leqslant X_{sep}\leqslant 11.25~\text{cm}$ , and for each configuration the simulation lasted approximately 7 s, with a simulation time step of $2.7\times 10^{-6}~\text{s}$ . We registered force and position data at a rate of 21 717 samples per second which, for simulations of 7 s, makes the standard error of the mean of any averaged value to be extremely small. That is why in most of the plots the error is smaller than the symbol size.

Figure 2. (a) and (b) Lift force on the right intruder as a function of time for a single simulation corresponding to, (a) $U_{\infty }=1.1~\text{m}~\text{s}^{-1}$ and $X_{sep}=2.5~\text{cm}$ and (b) $U_{\infty }=0.29~\text{m}~\text{s}^{-1}$ and $X_{sep}=5~\text{cm}$ . The thick black lines are sliding averages (with a time window of 0.05 s) showing attraction and repulsion. The force is registered at a rate of 21 717 samples per second. (c) Average lift force $\langle F_{L}\rangle$ as a function of intruders’ separation $X_{sep}$ for different flow velocities $U_{\infty }$ . Positive lift means repulsion. Inset: drag force as a function of $U_{\infty }$ . Drag is reduced when there are two intruders, but it does not depend on $X_{sep}$ . (d) $\langle F_{L}\rangle$ as a function of $U_{\infty }$ for intruders at different separation distances. The point with $X_{sep}/d=3.6$ corresponds to a single simulation doubling the size of the grains, which suggests that the quantity $X_{sep}/d$ is the relevant parameter concerning the lift force, although more simulations are needed to be conclusive on this point. The error bars (standard error of the mean) are smaller than the size of the symbols.

4 Lift and drag

When the flow begins, the lift and drag forces exhibit a transient behaviour, which lasts for less than half a second, after which the forces stabilise and fluctuate around a well-defined average value (see figure 2 a,b). Interestingly, the average lift force can be either attractive or repulsive, depending on the gap between the intruders and the flow velocity (figure 2 c,d): a repulsive lift is associated with a small gap and high velocities. At intermediate separation distances, the lift becomes attractive, and reaches a minimum (maximum attraction), and tends to zero at larger $X_{sep}$ values (we could not test greater separation distances to determine the exact range of intruder–intruder interactions because we did not want the distance from the intruders to the walls to be less than the size of the gap between the intruders). This behaviour reveals that there are two competing mechanisms involved in determining lift. Later, we demonstrate that repulsion results from the jamming of grains flowing between the intruders. In contrast, the origins of the attraction are not yet clear, but we can demonstrate that a Bernoulli-like effect is unlikely to explain either the attraction or repulsion, as has been suggested in other studies (Zuriguel et al. Reference Zuriguel, Boudet, Amarouchene and Kellay2005; Pacheco-Vazquez & Ruiz-Suarez Reference Pacheco-Vazquez and Ruiz-Suarez2010; Solano-Altamirano et al. Reference Solano-Altamirano, Caballero-Robledo, Pacheco-Vázquez, Kamphorst and Ruiz-Suárez2013).

Figure 3. Average drag and lift forces versus $L$ , $h$ and $D$ (defined in figure 1). The data in (b) and (e) are the only ones that follow a clear linear relation and we performed a linear fit. The quantity $h_{eff}$ refers to the effective height of the column above the intruders due to the rain of grains impacting the surface of the granular column (see § 6). Simulations with smaller grain size $d$ are very demanding computationally so we have only performed one simulation doubling the size of the grains for $X_{sep}=3.75~\text{cm}$ (see figure 2 d).

The behaviour of the average drag force $\langle F_{D}\rangle$ is also interesting – compared to the drag that occurs with a single intruder, drag is reduced when there is another intruder placed alongside (see inset in figure 2 c). Surprisingly, this drag reduction does not depend on the gap size between the intruders, at least not for the separation distances that we explored. This reflects the fact that the range of interaction of the two intruders is much larger than the lateral size of our container and therefore lateral confinement plays an important role. Indeed, lift and drag forces are very sensitive to variations in the distance $L$ between the lateral walls (see figure 3 a,d).

The inset in figure 2(c) shows that the drag force, on either one or two intruders, is almost constant despite variations in the flow velocity. This indicates that our system, even for the fastest flows, is a quasistatic regime dominated by friction and gravity, in contrast with the inertial regime that is characterized by a dependence of the force on the square of the flow velocity (Faug Reference Faug2015). However, although it is clear that our system is not in the inertial regime, the small variation of the drag force as a function of $U_{\infty }$ indicates that it is not in a purely quasistatic regime either, where stresses are rate independent and can be described with Mohr–Coulomb friction models. Instead, this behaviour, which has been observed in other granular systems (Geng & Behringer Reference Geng and Behringer2005; Caballero-Robledo & Clément Reference Caballero-Robledo and Clément2009; Kneib et al. Reference Kneib, Faug, Dufour and Naaim2016), seems to be consistent with the inertial rheology for sheared dense granular materials, where the bulk friction of the material is controlled by the non-dimensional inertia number $I=t_{micro}/t_{macro}$ , where $t_{micro}$ is a microscopic time scale related to rearrangements and $t_{macro}$ is a macroscopic time scale related to shear rate (Andreotti, Forterre & Pouliquen Reference Andreotti, Forterre and Pouliquen2013; Kneib et al. Reference Kneib, Faug, Dufour and Naaim2016). It is worth noting that the flow behind the obstacles may be in the inertial or collisional regime, but it is well established that the force on circular objects immersed in granular flows depends solely on the grains that are upstream from the centre of the object (Ding et al. Reference Ding, Gravish and Goldman2011).

Figure 4. (a–d) Dependence on flow velocity of the difference in packing fraction times the square of the flow velocity on both sides of the intruders, calculated in the regions adjacent to the upper quadrants of the intruders (see inset in a), and divided by a power of the flow velocity with different exponents. (e–h) Dispersion of data around the best linear fit in a plot of the average lift force versus ${\rm\Delta}({\it\phi}v^{2})$ divided by the corresponding power of $U_{\infty }$ . It is clear that the expression that best represents $\langle F_{L}\rangle$ , among those tested here, is ${\rm\Delta}({\it\phi}v^{2})/U_{\infty }^{3/2}$ .

5 Local flow velocity

What are the roles played by the local flow velocity and local velocity fluctuations around the intruders in the mechanisms responsible for attraction and repulsion? To tackle this question we define square regions where field quantities such as velocity, velocity fluctuations and mass density can be measured. Determining the positions and size that these regions must have is important for making relevant measurements. We define ‘in’ regions as those that are on the side of the gap between the intruders, and ‘out’ regions as those that are in between the intruders and the walls (inset in figure 4 a). This is reasonable, since we were hypothesising that the lift force is related to the difference between the field quantities in the ‘in’ and ‘out’ regions. In these regions we measure the mean grain velocity $\langle v(x,y)\rangle =(1/n(x,y,{\it\tau}))\sum _{t=t_{1}}^{{\it\tau}}(\sum _{i}v_{i}(x,y))_{t}$ , where $v_{i}(x,y)$ is the magnitude of the velocity of the $i$ th grain that is inside the region centred at $(x,y)$ at time $t$ , $t_{1}$ is the time at which the flow can be considered stationary, ${\it\tau}$ is the end of the simulation and $n(x,y,{\it\tau})$ is the total number of grains counted on both sums; and we also measure the average mass density $\langle {\it\rho}(x,y)\rangle \approx {\it\rho}_{g}\langle {\it\phi}(x,y)\rangle _{t}$ , where ${\it\phi}(x,y)$ is the area fraction occupied by the grains in the region of measurement at a given time. The difference in the field quantities in the ‘in’ and ‘out’ regions that we found useful were ${\rm\Delta}\langle v\rangle ^{2}=\langle v(\text{out})\rangle ^{2}-\langle v(\text{in})\rangle ^{2}$ and ${\rm\Delta}({\it\phi}v^{2})=\langle {\it\phi}(\text{out})\rangle \langle v(\text{out})\rangle ^{2}-\langle {\it\phi}(\text{in})\rangle \langle v(\text{in})\rangle ^{2}$ . Interestingly, we found that the flows around or below the intruder equator are not related to the forces on the intruder. This finding is consistent with the lift on a single intruder moving horizontally in a granular medium (Ding et al. Reference Ding, Gravish and Goldman2011).

Zuriguel et al. (Reference Zuriguel, Boudet, Amarouchene and Kellay2005) suggested that the origin of repulsive lift forces between intruders is a Bernoulli-like effect, where a difference in flow velocity on both sides of the intruders results in a pressure difference with the form ${\rm\Delta}p={\it\rho}\langle v(\text{out})\rangle ^{2}-{\it\rho}\langle v(\text{in})\rangle ^{2}$ . In our system, the quantity ${\rm\Delta}\langle {\it\phi}v^{2}\rangle$ , calculated in the regions adjacent to the upper quadrants of the intruders and plotted versus flow velocity, reproduces the lift force behaviour reasonably well, especially for fast flows and for both repulsion and attraction (figure 4 a,e) (if the quantity ${\rm\Delta}\langle v^{2}\rangle$ is plotted instead of ${\rm\Delta}\langle {\it\phi}v^{2}\rangle$ versus $U_{\infty }$ , the plot does not change substantially; is like multiplying $v^{2}$ by a factor of the order of $0.7$ (see figure 2 in the supplementary material)). However, the similarity with the lift force is much better when the velocity difference is divided by $U_{\infty }^{3/2}$ (see figure 4 c,g), which is also better than dividing by $U_{\infty }^{1}$ or $U_{\infty }^{2}$ (see figure 4). In addition, if measurements are taken in the regions shown in figure 5(a), the similarity with the lift force is even better (see figure 5). Therefore, although the difference in the square of the flow velocity between the ‘in’ and ‘out’ regions seems to be related to lift force, its relation with pressure does not seem to be solely captured by multiplying by ${\it\rho}$ , as suggested by Zuriguel et al. (Reference Zuriguel, Boudet, Amarouchene and Kellay2005). The fact that ${\rm\Delta}\langle {\it\phi}v\rangle ^{2}$ must be divided by $U_{\infty }^{3/2}$ suggests that the origin of repulsion and attraction is not a Bernoulli effect. As such, we are obliged to conduct a dimensional analysis.

Figure 5. (a) Regions of measurements found empirically by trial and error where the relation between the field quantities and the lift force were the best. The squares have sides of size $D/2$ . (b,c) Field quantity ${\rm\Delta}({\it\rho}_{g}{\it\phi}v^{2})/U_{\infty }^{3/2}$ as a function of $X_{sep}$ and $U_{\infty }$ is very similar to lift force in figure 2. (d) Dispersion of data around the best linear fit in a plot of the average lift force versus ${\rm\Delta}({\it\rho}_{g}{\it\phi}v^{2})/U_{\infty }^{3/2}$ . The fit is better than in figure 4(g), which shows that the small regions in (a) are better than those in figure 4(a).

Figure 6. Average lift force and local pressure difference ${\rm\Delta}({\it\rho}_{g}{\it\phi}v^{2})$ arranged following (6.1) to study the dependence of ${\it\xi}$ on (a) $L$ , (b) $h_{eff}$ and (c) $D$ . The lines are only qualitative guides to the eye to compare data with a linear behaviour. The value of $F_{0}=-3.5\times 10^{-4}~\text{N}$ was taken from the fit in figure 7(a).

6 Dimensional analysis

Based on the data presented in the previous sections, we propose an equation for the average lift force as follows:

(6.1) $$\begin{eqnarray}\langle F_{L}\rangle =\frac{f({\it\mu}_{eff}){\it\xi}}{Fr^{3/2}}S{\rm\Delta}({\it\rho}_{g}{\it\phi}v^{2})+F_{0},\end{eqnarray}$$

where $f$ is a dimensionless function of the effective friction coefficient ${\it\mu}_{eff}$ , $S=dD$ is the relevant obstacle surface, $Fr=U_{\infty }/\sqrt{gd}$ is the Froude number, ${\it\xi}={\it\xi}(h,d,D,L)$ is a dimensionless function that depends on the geometry of the system and $F_{0}$ is a function of all the possible variables of the system but varies weakly with $X_{sep}$ and $U_{\infty }$ . By performing variations on $h$ , $D$ and $L$ (see figure 6), we find that, at least as a first-order approximation, the function ${\it\xi}$ has the form ${\it\xi}={\it\xi}(L,h_{eff},D^{-3/2})$ , where $h_{eff}$ is the effective height of the column above the intruders. Indeed, because of the periodic boundary conditions, the hydrostatic pressure in the granular medium is augmented with respect to the ${\it\rho}gh$ value by the pressure exerted by the ‘rain’ of grains impacting the surface of the granular column after passing the obstacles at the bottom of the cell. This extra pressure can be expressed in the form of a virtual column of height $h_{r}=13.8~\text{cm}$ , so that the pressure would be $p_{eff}={\it\rho}g(h+h_{r})={\it\rho}gh_{eff}$ with $h_{eff}\sim 48~\text{cm}$ (see figure 3 b). With these dimensional constraints, the expression for ${\it\xi}$ in (6.1) can be written as

(6.2) $$\begin{eqnarray}{\it\xi}=\left(\frac{h_{eff}}{D}\right)\left(\frac{L}{D}\right)\left(\frac{D}{d}\right)^{1/2}.\end{eqnarray}$$

When the average lift force is plotted against the term ${\it\xi}Fr^{-3/2}S{\rm\Delta}({\it\rho}_{g}{\it\phi}v^{2})$ , with ${\it\xi}$ given by (6.2), the data follow a better linear pattern than those presented in the previous sections (see figure 7 a). From the linear fit in figure 7(a) we obtain the missing parameters in (6.1): $f({\it\mu}_{eff})=4.6\times 10^{-2}$ and $F_{0}=-3.5\times 10^{-4}~\text{N}$ . We note that most of the data plotted in figure 7(a) correspond to simulations in which systematic variations are made on $X_{sep}$ and $U_{\infty }$ , but the plot also includes simulations where $D$ , $L$ , $h$ and $d$ vary, which demonstrates the robustness of (6.1) and (6.2) as a first-order approximation of those variables.

Figure 7. (a) Dispersion of data around the best linear fit in a plot of the average lift force versus ${\it\xi}Fr^{-3/2}S{\rm\Delta}({\it\rho}_{g}{\it\phi}v^{2})$ . Most of the data correspond to systematic variations of $X_{sep}$ and $U_{\infty }$ , but the plot also includes data corresponding to variations of $D$ , $L$ , $h$ and $d$ (see figure 1 for their definitions). From the linear fit we get the parameters $f({\it\mu}_{eff})=4.6\times 10^{-2}$ and $F_{0}=-3.5\times 10^{-4}~\text{N}$ of (6.1). (b) Comparison between the average lift force calculated from the force on intruders and the lift force obtained from (6.1) for three different velocities.

Figure 7(b) shows the average lift force calculated, either directly from the force on intruders, or using (6.1), as a function of $X_{sep}$ for three different velocities. The reproduction of the lift force through the difference in the flow velocity in the ‘in’ and ‘out’ regions works remarkably well, although it is not perfect.

There may be a number of reasons for the dispersion of the data observed in figure 7. First, the term $F_{0}$ in (6.1) may be a function of all the system variables, although its contribution to the lift force is small. The fact that $F_{0}$ varies weakly with $X_{sep}$ and $U_{\infty }$ indicates that it might stem primarily from intrinsic properties of the granular medium. In addition, negative $F_{0}$ means that there exists an attractive lift force even when the local pressure difference given by ${\rm\Delta}({\it\rho}_{g}{\it\phi}v^{2})$ is zero. The origin of this attraction in a symmetric flow remains for the moment as an open question, but we could speculate that fluctuations could play a role here. Second, the expression for the prefactor ${\it\xi}$ in (6.2) is an empirical formula that is at best a first-order approximation of the real behaviour. This follows at once from the fact that it is not compatible with the two limiting cases in which the variables become very large or very small. For example, it is clear that the lift force cannot be infinite for a container of infinite width $L$ ; instead, it should become saturated at some point, probably following the pressure saturation, such as that described by Janssen in a silo (Pacheco-Vázquez et al. Reference Pacheco-Vázquez, Caballero-Robledo, Solano-Altamirano, Altshuler, Batista-Leyva and Ruiz-Suárez2011). An interesting result in our system is that the lift force changes sign when we double the size of the grains in the medium (see figure 2 d). This shows a strong dependence of the local flow velocity, the term ${\rm\Delta}({\it\rho}_{g}{\it\phi}v^{2})$ , on the size of the grains $d$ . Therefore, the scaling of the lift force as $d^{5/4}$ (almost linear) predicted by (6.1) and (6.2) is probably not the only contribution by $d$ to the lift force. Larger forces for larger grains have been reported on objects dragged horizontally (Albert et al. Reference Albert, Pfeifer, Barbasi and Schiffer1999; Ding et al. Reference Ding, Gravish and Goldman2011), the result of which is at odds with the observation reported in other work (Chehata et al. Reference Chehata, Zenit and Wassgren2003; Seguin et al. Reference Seguin, Bertho, Gondret and Crassous2011), in which the drag on an object immersed in a vertical granular flow increases when the size of the grains decreases. Clearly, the role of grain size on drag and lift forces is not well understood and will require further detailed study. Finally, the linear dispersion of data in figure 7 may also be related to the choice of the regions where the flow is measured, which in our case was determined by trial and error. There is no guarantee that they represent the best choices until an analytical derivation can be made of the relation between lift and local flow.

Figure 8. Average flow velocity at the middle point between intruders, $v_{mp}$ , plotted as a function of flow velocity $U_{\infty }$ . The dashed line is the average velocity expected from continuity at the middle area of the intruders due to the reduction in the cross-section: $v_{eq}=LU_{\infty }/(L-2D)$ . A comparison with the data in figure 2(d) shows that repulsive lift is related to the saturation of $v_{mp}$ due to jamming.

7 Jamming between intruders

The variation of the average lift force as a function of $X_{sep}$ reveals the competition of two mechanisms, one of which causes attraction and the other repulsion between the intruders. These two regimes are evident in figure 2(c), which also shows that the repulsive mechanism is favoured by a small intruder separation and large flow velocities. It has been proposed that the origin of repulsion may be the jamming of particles trying to flow through the gap between intruders (Pacheco-Vazquez & Ruiz-Suarez Reference Pacheco-Vazquez and Ruiz-Suarez2010; Solano-Altamirano et al. Reference Solano-Altamirano, Caballero-Robledo, Pacheco-Vázquez, Kamphorst and Ruiz-Suárez2013), but as yet there has been no direct proof. Now, however, with our current set-up, we can support this idea by measuring the average flow velocity at the middle point between the intruders, $v_{mp}$ , and plotting it as a function of $U_{\infty }$ (figure 8). When we compare this plot with the lift force as a function of $U_{\infty }$ (figure 2 d), we observe that whenever $v_{mp}$ varies linearly with $U_{\infty }$ , there is attraction between the intruders, while repulsion is correlated with $v_{mp}$ remaining constant. When the separation between the intruders is 2.5 cm, as in figure 8, the results are especially illustrative in showing how $v_{mp}$ transits from linearity to a constant value when $U_{\infty }$ increases. This saturation of $v_{mp}$ to a constant value is favoured by small gaps between intruders and by fast flows, as is the case for the repulsive lift force. Therefore, it is a signature of the jamming of grains trying to flow through the gap between intruders, imposing the maximum possible flow through the gap and pushing apart the intruders, which causes repulsion when $U_{\infty }$ increases.

What, then, causes attraction? Figure 8 shows that the proximity of the intruders, when jamming has not occurred, accelerates the flow between the intruders and favours attraction. Moreover, the lower the velocity, the stronger the attraction. Therefore, the presence of a neighbour intruder at an intermediate distance facilitates flow in the gap between the intruders, an effect that is enhanced at slow flow velocities. This scenario supports that proposed by Seguin et al. (Reference Seguin, Bertho, Gondret and Crassous2011), where the drag force on a cylinder penetrating a granular medium can be understood as analogous to a hot cylinder penetrating a fluid whose viscosity decreases when the temperature increases, i.e. when the flow is slow, the temperature between the intruders increases further and the viscosity is at its lowest. This means the extension of the kinetic theory to dense granular flows is a good candidate for describing drag and lift in systems with multiple intruders.