# Mean force and fluctuations on a wall immersed in a sheared granular flow

François Kneib, Thierry Faug, Frédéric Dufour, Mohamed Naaim

► To cite this version:

François Kneib, Thierry Faug, Frédéric Dufour, Mohamed Naaim. Mean force and fluctuations on a wall immersed in a sheared granular flow. Physical Review E, 2019, 99 (5), pp.052901. 10.1103/Phys-RevE.99.052901. hal-02123563

HAL Id: hal-02123563

https://hal.science/hal-02123563

Submitted on 18 Jun 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mean force and fluctuations on a wall immersed in a sheared granular flow

François Kneib* and Thierry Faug†

Univ. Grenoble Alpes, Irstea, UR ETGR, 38000 Grenoble, France

Frédéric Dufour‡

Univ. Grenoble Alpes, CNRS, INP-Institute of Engineering

Univ. Grenoble Alpes, 3SR, 38000 Grenoble, France

Mohamed Naaim§

Univ. Grenoble Alpes, Irstea, UR ETGR, 38000 Grenoble, France

(Dated: March 20, 2019)
Abstract

In a sheared and confined granular flow, the mean force and the force fluctuations on a rigid wall
are studied by means of numerical simulations based on the discrete element method. An original
periodic immersed-wall system is designed to investigate a wide range of confinement pressure and
shearing velocity imposed at the top of the flow, considering different obstacle heights. The mean
pressure on the wall relative to the confinement pressure is found to be a monotonic function of the
boundary macroscopic inertial number which encapsulates the confinement pressure, the shearing
velocity and the thickness of the sheared layer above the wall. The one-to-one relation is slightly
affected by the length of the granular system. The force fluctuations on the wall are quantified
through the analysis of both the distributions of grain-wall contact forces and the autocorrelation of
force time series. The distributions narrow as the boundary macroscopic inertial number decreases,
moving from asymmetric log-normal shape to nearly Gaussian-like shape. That evolution of the
grain-wall force distributions is accompanied at the lowest inertial numbers by the occurrence of
a system memory in terms of the force transmitted to the wall, provided that the system length
is not too large. Moreover, the distributions of grain-wall contact forces are unchanged when the
inertial number is increased above a critical value. All those results allow to clearly identify the
transitions from quasi-static to dense inertial, and from dense inertial to collisional, granular flow
regimes.

<footnote>* francois.kneib@irstea.fr</footnote>

<footnote>† thierry.faug@irstea.fr</footnote>

<footnote>‡ frederic.dufour@3sr-grenoble.fr</footnote>

<footnote>§ mohamed.naaim@irstea.fr</footnote>
## I. INTRODUCTION

Many industrial processes and geophysical problems involve granular matter at rest or in movement. More specifically, the force exerted on solid objects subjected to flowing granular matter is an important topic that has attracted increasing attention over the recent years. The time-averaged force on objects of different sizes and shapes was studied in various configurations (see for instance [1–5]). The drag force on immersed solid objects in dense free-surface granular flows was firstly investigated by the experimental work of Wieghardt [1]. In his pioneering experiments, Wieghardt found that the drag force was only weakly influenced by the flow velocity and observed that a quasi-static pile of grains formed upstream the solid object. Those findings were confirmed later under low to intermediate shear rates in different configurations, including annular shear flows [2], free-fall chute flows [3], immersed rods rotating along an axis perpendicular to its main axis [4], etc. All those studies showed that, whenever the obstacle is in contact with a quasi-static bulk of grains, the drag force linearly depends on both the confinement pressure and the object size rather than on the flow velocity. This picture becomes however no longer valid once the granular flows become fast, more dilute and collisional, as shown for instance in [6, 7]. In this situation, it was found that the mean force acting on the obstacle evolved with the square of the relative grains-object velocity.

The question of the force fluctuations on objects immersed in granular flows is still poorly investigated. Most studies dedicated to force fluctuations in granular materials concerned quasi-static or slowly sheared systems, in Couette-like geometries (see [8–10]). Only a few number of papers addressed the problem of force fluctuations for rapid granular flows. The experiments were based on the discharge of silos equipped with a force sensor, at a boundary wall, whose size was the typical size of the grains [11, 12]. The measurements of contact force distributions in quasi-static or slowly sheared systems were initially computed with the help of carbon paper placed at a system boundary [13], but the use of photo-elastic disks in laboratory experiments, as well as numerical modelling, opened the path to measuring the contact forces inside the granular bulk. Some cross-comparisons were performed in [14–17], showing that—at least in the static case—the grain-grain contact force distributions from the bulk were identical to the grain-wall contact force distributions. A robust finding was established for the shape of the contact forces *f* distributions for values greater than
the mean $\bar{f}$: they decrease exponentially up to the highest values of $f$ (see for instance [11, 18, 19]). This behaviour is sensitive to the shear velocity (or the deviatoric stress tensor in static cases), as the inertia tends to broaden the distributions (see [9, 11, 16]). For the lowest forces ($f < \bar{f}$), different behaviours depending mostly on the sample history are observed from vanishing to very high probabilities near $f = 0$. Some empirical or theoretical models were proposed to fit the entire range of observed distribution shapes for quasi-static systems [13, 20, 21]. A few studies tackled the distributions of forces that were integrated through a spatial dimension for slow to intermediate flow rates. It means that the data sets were systematically recorded on surfaces that were greater than the grain size. The configurations tested were the Couette geometry [9, 20, 22] and a rotating immersed rod [23]. A remarkable trend emerged from all those studies: while increasing the spatial scale considered, the exponential force distributions narrow and finally move towards Gaussian-like distributions. The latter result is especially pronounced for very slow flows, whereas intermediate shear rates result in exponential to log-normal shapes.

Most of the past studies on granular force fluctuations addressed a narrow range of boundary conditions in terms of shearing velocities imparted to the granular sample. As such, the transitions between the different flow regimes (quasi-static, dense inertial, collisional) and their links to force fluctuations still remain poorly investigated. The overarching aim of the present study is to bridge this gap. We propose to investigate the force exerted on a wall-like obstacle immersed in a sheared granular flow by means of numerical modelling. The originality of this work is the wide range of flow regimes investigated, from nearly quasi-static through dense inertial to collisional flows. In order to control the flowing conditions and to focus on the grain-wall interactions, an innovative planar system is designed. While the spherical grains are trapped across the direction perpendicular to the mean flow between two rough walls, a periodic boundary condition is used along the mean flow (shearing) direction. The obstacle is a wall which is fully immersed in the granular bulk and orthogonal to the shearing direction. Note that one initial motivation of the present study was to mimick the problem of gravity-driven free-surface granular flows passing over a rigid wall [24–26]. First, the modelled system (boundary conditions, relevant macroscopic parameters), as well as the numerical simulations based on the discrete element method, are described in Sec. II. Then, the mean (time-averaged) dynamics is analysed by investigating the velocity fields inside the system, the scaling of the mean force experienced by the wall, and the spatial maps of the
92

93

94

95

96

97

98

local rheology (Sec. III). The force fluctuations on the immersed wall are quantified through both the autocorrelation of mesoscopic force signals associated with the system memory (Sec. IV A), and the distributions of force time series at macroscopic and mesoscopic spatial scales (Sec. IV B). The effect of the system length relative to its height is discussed in Sec. V. Sec. VI ends the paper with an extended discussion on the main results and perspectives, including potential applications of this research to the problem of civil engineering structures exposed to natural granular flows.

99

# II. DEM SIMULATION OF A GRANULAR FLOW OVER AN IMMERSED WALL

100

## A. Macroscopic boundary conditions

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

A planar sample of spheres is confined between two walls at the top and the bottom, and by a periodic boundary condition linking the right and the left borders of the system, as depicted by the vertical dashed lines in Fig.1. In the present study, a three-dimensional numerical code [27] is used considering spherical grains for which the kinematics into the $z$-direction (normal to the plane $(x,y)$ in Fig.1) remains at zero. The bottom wall is fixed, rough, and its length is equal to the system length $L$. A constant macroscopic pressure $P$ is applied to the system. It is controlled by the vertical position of the rough top wall, continuously computed by the DEM algorithm and applying a vertical force on the sample. The periodicity of the system also concerns the top wall which has a constant horizontal velocity $U$ shearing the sample. The top wall remains perfectly horizontal by preventing its rotation in the $(x, y)$ plane (see Fig. 1). The roughness of the top and bottom walls is made of aligned spheres. A thin, rigid and smooth vertical wall of finite height $h$, namely the obstacle, is immersed into the spheres assembly. The spheres are initially deposited under gravity to form a rectangular dense packing of height $H$ greater than the obstacle height $h$. Then, the gravity is set off, the confinement pressure applied at the top and the shearing starts for a period of ten seconds over which the data is collected. Only the data arising from the latest nine seconds is exploited because the first second is the typical observed time needed for the system to reach a permanent regime.

119

120

Although the system is gravity-free, the following dimensionless macroscopic parameters are defined for convenience and used to run the numerical simulations over a broad range of
<figure><img src="image_2.png" /><figcaption>FIG. 1. Geometry and boundary conditions of the periodic sheared system simulated by DEM. The granular sample is trapped between two horizontal walls and a periodic boundary condition in the horizontal x-direction (thick dashed lines). The top wall applies the shear displacement U and the confinement pressure P (see text for details), and force measurements are done onto the sidewall of height h, immersed into the granular sample.</figcaption></figure>

both confinement pressure P and shear velocity U:

$$N_P = \frac{1}{\Phi \rho_p g H} P, \qquad (1)$$

$$N_U = \frac{1}{\sqrt{gd}} U. \qquad (2)$$

The particle density was taken equal to $\rho_P = 2500 \text{ kg m}^{-3}$ and the gravity acceleration is $g = 9.81 \text{ m s}^{-2}$. The mean grain diameter $d$ was taken equal to 1 mm and $H/d$ was kept constant $H/d = 30$. A slight polydispersity was introduced by picking randomly the grain diameters between $0.85d$ and $1.15d$, in order to avoid crystallisation effects on the one side (no polydispersity) and migration/segregation processes on the other side (too high polydispersity). An arbitrary constant macroscopic volume fraction $\Phi = 0.6$ is considered here, which roughly corresponds to the random close packing of a three-dimensional assembly of spheres. We recall here that the numerical simulations presented in this study use spherical particles whose centers stay, by construction, on a planar surface [plane (x,y) on Fig. 1]. The choice of defining such dimensionless numbers $N_P$ and $N_U$ based on $g$ is purely arbitrary, has no influence on the results presented in the following, and was essentially motivated by advancing knowledge on the problem of the mean force and force fluctuations exercised by free-surface gravity-driven granular flows on rigid walls, as already discussed in the introduction.
Moreover, it followed the choice already made in a previous study on a different system— namely the granular lid-driven cavity—but for which very similar measurements and analysis were made [28]. Considering those definitions thus allows a direct comparison between the two systems (see discussion in Sec. VI B).

A boundary macroscopic inertial number $I_M$ can be defined from the typical time associated with the top confinement pressure $t_P = d\sqrt{\rho/P}$ and the typical time equal to the inverse of the macroscopic shear rate $t_U = H/U$:

$$I_M = \frac{t_P}{t_U} = \frac{d\sqrt{\rho/P}}{H/U} = \frac{d\sqrt{\Phi\rho_P}}{H} \frac{U}{\sqrt{P}}, \quad (3)$$

where $\rho = \Phi\rho_P$ holds for the bulk density of the granular sample. It is worthy to stress here that the dimensionless number $I_M$ defined above accounts for the variation of both $P$ and $U$ but that $H$ and $\Phi$ are considered as constant and equal to the initial values before shearing the granular sample ($H = 0.03$ m and $\Phi \simeq 0.6$ for a 3D assembly of spheres). Actually, the numerical system studied here is by construction volume-free and therefore some slight variation of the volume fraction is possible. The mean volume fraction of the granular sample which was actually measured in the DEM simulations will be discussed in the concluding section of the paper (see Sec. VI A).

In the present study, the dimensionless numbers $N_P$ and $N_U$ were varied from 0.01 to 100 and from 1 to 20, respectively. This allowed us to investigate a wide range of both $P$ (from about 4.4 to 4.4 $10^4$ N m$^{-2}$) and $U$ (from about $10^{-2}$ to 2 m s$^{-1}$), resulting in inertial numbers $I_M$ that ranged typically from $1.2 \times 10^{-3}$ to 0.6 (considering the constant values discussed above for $H$, $d$, $\rho_P$ and $\Phi$). Simulations with different obstacle heights ($h \in \{5d, 10d, 15d, 20d\}$) were performed in order to analyse the potential influence of this parameter on the time-averaged dynamics. The length of the periodic cell relative to the initial height of the sample, $L/H$, is a parameter that may also influence the system dynamics. Most of the simulations presented in this study were done with $L/H = 5$. The slight sensitivity of the results to doubling $L/H$ (from 5 to 10) will be discussed in Sec. V.

## B. Micromechanical parameters

The discrete elements method (DEM), the contact laws and the associated parameters used in this paper are the same as described in [28, 29]. The general DEM algorithm is
described in [30]. A visco-elastic contact law is used for the normal component of the inter-particle forces (see [31]), while a classical coulomb friction threshold drives the tangential component in the same fashion as in [32]. The following set of equations describes the micromechanical model used in the present study, for the grain-grain and the grain-wall interactions:

$$ \begin{cases} \mathbf{F}_n = \max(0, k_n \delta_n + c_n \dot{\delta}_n) \mathbf{n}, \\ d\mathbf{F}_s = (k_s \dot{\delta}_s dt) \mathbf{s}, \\ |\mathbf{F}_s| \le \mu |\mathbf{F}_n|, \end{cases} \quad (4) $$

where $\mathbf{n}$ and $\mathbf{s}$ are the unity vectors along the contact normal and shear directions respectively, $k_n$ and $k_s$ are the normal and tangential contact stiffnesses, $\delta_n$ is the normal penetration depth, $\dot{\delta}_s$ is the tangential displacement increment, $\mu$ is the local friction coefficient, $c_n$ is the normal viscosity coefficient and $dt$ is the time-step. The four physical parameters $k_n$, $k_s$, $c_n$ and $\mu$, are chosen to fit the behaviour of glass beads. The contact stiffness is reduced to decrease the total time of calculation, but the limit of rigid grains was systematically respected (see details in [28, 29]). The coefficient $c_n$ is set in the same manner as in [29] with a restitution coefficient $e = 0.5$, and $\mu$ was taken equal to 0.5.

## III. MEAN DYNAMICS

The present section investigates the time-averaged dynamics of the system: the streamlines within the granular sample and the vertical velocity profiles along the sample length (Sec. III B), the mean force $\overline{F}$ on the obstacle—the latter being the wall on the right side of the system, which faces the shearing velocity at the top (Sec. III A), and the local rheology within the granular sample (Sec. III C). In the rest of the paper, all time-averaged values are calculated over nine seconds, excluding the first second of simulations in order to avoid any non-permanent state.

### A. Bulk kinematics

The time-averaged velocity streamlines within the granular sample are displayed in Fig. 2, considering three different values of the macroscopic inertial number $I_M$ and one value of the obstacle height $h = 20d$. The norm $v_g$ of the local grain velocity, scaled by the top wall
192 shearing velocity $U$, is displayed on the colour scale. For all simulations the flow velocity 193 fits the top and bottom boundary conditions, ranging from zero near the bottom to $U$ in 194 the vicinity of the top wall. At low inertial numbers ($I_M \simeq 1.2 \ 10^{-3}$) the flow velocity 195 shows a symmetric pattern along the system length. Near both corners at the bottom of the 196 system, very low flow rates take place in front of the obstacle, thus producing streamlines 197 that are rather chaotic (due to the likely intermittent flows taking place here). The latter 198 are the markers of nearly triangular dead zones, forming near the obstacles (the wall of 199 interest on the right side and its duplicate for periodic conditions on the left side: see 200 Fig. 1) and whose heights gradually decrease while moving farther from the obstacle. In 201 a central region ($50d \lesssim x/d \lesssim 100d$) the streamlines are horizontal and parallel with each 202 other, revealing that the flow experiences a quasi-horizontal shear without any significant 203 vertical displacement. Back to past studies on free-surface gravity-driven granular flows 204 passing over a rigid wall [24–26], that homogeneous region found in our simulated system 205 here mimics the incoming flow, co-existing with a similar—nearly triangular—dead zone 206 formed upstream the obstacle (the wall on the right side).

207 The response of the velocity streamlines to the increase of $I_M$ is tricky. Figure 2(a) 208 exhibits a symmetric pattern at the lowest $I_M$ but not for higher $I_M$ in Fig. 2(b) and 209 Fig. 2(c). As the obstacle is also duplicated through the periodic boundary condition along 210 $x$, a vortex of height $h$ starts forming behind the wall on the left side of the system. This 211 phenomenon is caused by the grains flowing through the periodic boundary condition and 212 coming into the system—from the top of the vertical left sidewall—with high velocities. 213 The shape of this vortex varies with the macroscopic boundary conditions. In particular, 214 the vortex length increases with $I_M$, disturbing progressively the picture of the central region 215 where shear is quasi-horizontal. The latter region becomes narrower when increasing $I_M$, 216 and may totally disappear at the highest $I_M$. This peculiar feature of the periodic immersed- 217 wall system studied here is further investigated in the following by focusing on the vertical 218 velocity in the granular bulk.

219 Figure 3 depicts the vertical component $V_y$ of the grains’ velocity along the system length, 220 averaged over time and height, and scaled by the top wall shearing velocity $U$. This mea- 221 surement is presented for seven values of inertial numbers, from quasistatic to collisional 222 flow regimes.

223 For $I_M \le 1.2 \ 10^{-2}$ all profiles collapse, revealing that the global (scaled) grain velocity is
<figure><img src="image_3.png" /><figcaption>FIG. 2. Time-averaged streamlines for an obstacle height of $h = 20d$ (thin vertical bars on each border of the system), for three different values of $I_M$: $I_M = 1.2 \ 10^{-3}$ (a), $I_M = 6.1 \ 10^{-2}$ (b) and $I_M = 0.3043$ (c). The colour scale represents the norm of the local grain velocity scaled by the top wall shearing velocity, $v_g/U$.</figcaption></figure>

not influenced by the macroscopic conditions in the quasistatic regime. Near the (left side) incoming flow region ($x/L < 0.3$), $V_y$ is negative with a high slope gradient, induced by the grains that continuously fill the bottom of the system. Between $x/L = 0.3$ and $x/L = 0.85$, $V_y/U$ is nearly close to zero. Only a very small positive slope gradient is detected. This region, which we call a plateau, characterizes the homogeneous sheared zone described just above. Finally, from $x/L = 0.85$, the increase of the slope gradient marks the presence of the dead zone formed upstream of the obstacle (on the right side) overtopped by the grains. The presence of a plateau in the center of the curves shown in Fig. 3 demonstrates that the system length should be large enough to guarantee that the flowing direction of the granular
<figure><img src="image_4.png" /><figcaption>FIG. 3. Profile of the vertical component of the grains velocity $V_y$, averaged over both time and granular sample height $H$, scaled by the top wall shearing velocity $U$. The profile is plotted for $h = 20d$ and seven values of $I_M$ from $1.2 \ 10^{-3}$ to $0.6$.</figcaption></figure>

assembly is orthogonal to the obstacle. Thus, the force exerted on the obstacle of interest (on the right side) should not be influenced by its periodic duplicate (on the left side).

As shown in Fig. 2(b) and Fig. 2(c), the vortex formed in the left side of the system produces vertical movements of the granular flow that are inhomogeneous through the system height for simulations done at $I_M \ge 6.1 \ 10^{-2}$. Therefore, the averaging over height done in Fig. 3 to obtain the $V_y/U$ profile for the corresponding inertial numbers cannot be interpreted quantitatively. On a more qualitative note, the plateau is still observed up to $I_M = 9.1 \ 10^{-2}$ but the two curves obtained for the largest $I_M$ ($I_M \ge 0.3$) are highly disturbed by the vortex, thus preventing the occurrence of the plateau in the collisional regime. It will be shown in Sec. V that doubling $L/H$ does not change this observation. It is not possible to produce a plateau for $L/H = 10$ at high inertial numbers. The influence zone downstream of the duplicate obstacle is too large so that it disturbs the granular kinematics in the whole volume of the system up to the obstacle of interest. However, it will be shown in the following of the paper that an analysis of the time-averaged force and the force fluctuations both experienced by the obstacle on the right side can be done continuously, over the whole
range of $I_M$ covered by this study. In particular, a robust scaling for the mean total force exerted on the obstacle is established in the next section.

## B. Mean force on the obstacle

The key feature of the flow dynamics observed in the immersed-wall system investigated here is the following: some grains initially settled in the granular sample of height $H$ are trapped upstream the obstacle (the right sidewall) of height $h$ smaller than $H$, thus forming a dead zone, whilst other grains are able to flow over the obstacle (see details of the flow kinematics studied in Sec. III A). Based on the height $H-h$ of the sheared granular layer passing over the obstacle, a new macroscopic inertial number $I_M^*$ can be defined as follows:

$$I_M^* = \frac{d\sqrt{\rho/P}}{(H-h)/U} = d\sqrt{\Phi\rho_P} \frac{U}{(H-h)\sqrt{P}} \quad (5)$$

In the above definition, the dead zone is assumed to cover the entire height of the wall. The grains passing over the wall are sheared over the vertical distance between the top wall (at height $H$) and the top of the dead zone of null velocity (at height $h$). Introducing this macroscopic inertial number $I_M^*$ allows obtaining a remarkable scaling for the mean force $\overline{F}$ on the obstacle, scaled by the pressure force $Phd$ associated with the wall height, as a function of $I_M^*$ (see Fig. 4). Note that in the above definition of $I_M^*$, we also consider constant values of $H$ and $\Phi$ that correspond to the typical values before shearing the sample (see previous discussion for the definition of $I_M$).

All the data collapse into one single curve for different obstacle heights $h \in \{5d, 10d, 15d, 20d\}$. Considering the inertial number $I_M$ initially based on the total height $H$ of the sample does not produce such a good collapse. This clearly shows that the height difference $H-h$, corresponding to the thickness of the granular layer able to overtop the sidewall, plays a crucial role in the force applied on the obstacle, in addition to the combined effects of confinement pressure $P$ and shearing velocity $U$ already encapsulated into the macroscopic inertial number.

The scaled mean force, $\overline{F}/(Phd)$, shows a quasi-linear evolution with $I_M^*$ at low and intermediate $I_M^*$. Although a scattering of the data appears at high $I_M^*$, a saturation then occurs from around $I_M^* = 0.3$ (see the linear plot in the inset of Fig. 4). The following relation fits well the DEM data, as shown by the solid line drawn on Fig. 4:
<figure><img src="image_5.png" /><figcaption>FIG. 4. Time-averaged force on the obstacle $\overline{F}$ scaled by $Phd$ as a function of $I_M^*$, the inertial number corresponding to the height of grains overtopping the obstacle. Four obstacle heights are presented: $h \in \{5d, 10d, 15d, 20d\}$. The solid line shows the prediction of Eq. (6) with $r_1 = 1.4$, $r_2 = 3.3$ and $I_M^{*0} = 0.3$. Inset: same data with linear x-axis scale.</figcaption></figure>

$$ \frac{\overline{F}}{Phd} = r_1 + \frac{r_2 - r_1}{1 + \frac{I_M^{*0}}{I_M^*}}, \qquad (6) $$

where $r_1 = 1.4$ and $r_2 = 3.3$ are the two asymptotic values of $\overline{F}/(Phd)$ at low and high $I_M^*$, respectively. The constant $I_M^{*0} = 0.3$ is the inertial number for which $\overline{F}/(Phd) = (r_1+r_2)/2$. It is worth noting that the values of $r_1$, $r_2$ and $I_M^{*0}$ were obtained for the configuration presented in this paper, and may be influenced by the micromechanical grain parameters ($d$, $\mu$, etc.), as well as by the system configuration. In particular, the sensitivity of the mean force scaling to doubling $L/H$ will be discussed in Sec. V.

It is worthwhile to note that the data saturation (concomitant with some scattering) and the plateau which vanishes in the mean $V_y$ profiles (see Sec. III A) both occur at $I_M^* \simeq 0.3$. This observation appears to be a robust signature of the transition from the dense inertial regime to a more collisional regime, where enduring frictional grain-grain contacts are replaced by short-life contacts due to the strong increase of the inertia and the increase of the mean free path of particles. The role played by the local $\mu(I)$-rheology (earlier proposed
by [33, 34]) on the force transmision towards the wall-like obstacle is further analysed below.

## C. The local $\mu(I)$-rheology

The form of the empirical scaling showing a one-to-one relation between the scaled mean force $\overline{F}/(Phd)$ on the wall and the boundary macroscopic inertial number $I_M^*$, given by Eq. (6), suggests that the local $\mu(I)$-rheology [33] comes into play. To further investigate this, the time-averaged values of the local inertial number $I$ and the local effective friction $\mu_{loc}$ were measured inside the whole volume of the granular system.

For that purpose, the strain ($D$) and stress ($\sigma$) tensors at local (grain) scale were computed, using a spatial Kernel smoothing method and tesselation techniques. The technical aspects of those calculations can be found in [28]. The pressure $p$ within the granular medium was defined by:

$$p = \frac{1}{3} \mathrm{Tr}(\sigma). \qquad (7)$$

The local effective friction was calculated as:

$$\mu_{loc} = \frac{1}{3} \frac{\mathrm{Tr}(\sigma' D')}{p \|D'\|}, \qquad (8)$$

where $D' = D - \frac{1}{3}\mathrm{Tr}(D)I_3$ holds for the deviatoric strain tensor and $\sigma' = \sigma - pI_3$ is the deviatoric stress tensor. We define $I_3$ as the identity matrix of size 3 and $\|A\| = \sqrt{\mathrm{Tr}(AA^T)/2}$ for any $A$. Finally, the local inertial number was computed using the relation:

$$I = d \sqrt{\frac{\rho_P}{p}} \|D'\|. \qquad (9)$$

Although the strain field is rather complicated inside the immersed-wall system investigated here (see the streamlines drawn on Fig. 2), it was generally observed that the strain and the stress tensors were quite well aligned within most of the system volume. Figure 5 shows an example for $I_M = 9.1 \times 10^{-2}$. Nonetheless, there exist two regions inside which the principal directions of stress and strain tensors are perpendicular to each other. First, a region at the top of the duplicate wall (on the left side of the system) exhibits principal directions of stress and strain tensors that are orthogonal. It is worthy to note that for lower $I_M$ the stress and strain tensors were found to be collinear in that specific region (not shown here). This peculiar behavior at high $I_M$ might be attributed to the complicated interaction between the flux of grains arriving at high velocity (from the left side) and the
<figure><img src="image_6.png" /><figcaption>FIG. 5. Spatial maps of the principal strain (a) and stress (b) orientations averaged over time within the immersed-wall system: example for $I_M = 9.1 \ 10^{-2}$.</figcaption></figure>

vortex dynamics which takes place just below, behind the duplicate wall (on the left side). This could be another signature of the transition towards the dilute/collisional regime, as already discussed when interpreting the vertical velocity profiles depicted in Fig. 3. Second, there exists a region on Fig. 5, just above the lid and at about two thirds of the length of the system, inside which the principal directions of strain and stress tensors are not aligned but orthogonal, regardless of $I_M$ (the results shown in Fig. 5 are general for all $I_M$ for that second region). If one looks more closely at that region, this is exactly where the grains are permanently forced to climb over the quasistatic dead zone that forms upstream of the wall (which stands on the right side of the system and faces the direction of the imposed velocity $U$ at the top) and, at the same time, pushed down by the upper rigid wall boundary condition. Although this is not detected on the time-averaged kinematics (Fig. 2), this competition may lead to the formation of a strong discontinuity in stresses. As a result, in that shock region, the time-averaged principal directions of the strain and stress tensors may be not aligned but perpendicular. This peculiar behavior would merit further investigation in the future. Apart from those two specific regions, there is collinearity between stress and strain tensors within most of the system volume. The variation of $\mu_{loc}$ with $I$ is analysed in the following.

For each position $(x, y)$, the measured $\mu_{loc}$ was compared to a $\mu_{th}(I)$ derived from the
relation proposed by Jop et al. [33], which reads as follows:

$$ \mu_{th} = \mu_1 + (\mu_2 - \mu_1) \frac{1}{1 + I_0/I}, \qquad (10) $$

where $\mu_1$, $\mu_2$ and $I_0$ are parameters dependent on the mechanical grain properties. In our study, the best fitting parameters for $\mu_{th}$ were found to be (after an iteration process) $\mu_1 = 0.17$, $\mu_2 = 0.48$ and $I_0 = 0.18$, as shown in Fig. 6. Figure 6 shows the time-averaged value of $\mu_{loc}$ as a function of $I$. It is important to stress that the data points shown on Fig. 6 were extracted from a region centered on the system and limited to a length of 90d in the x-axis direction, thus excluding the two regions (of length 30d each) close to the vertical side-wall and its periodic image on the left side of the system. The very bottom (from 0 to 5d in the y-axis direction) and very top (from $H-5d$ to $H$) regions were also excluded from the analysis. What happens inside the regions close to the side-wall (non-local effects) is further discussed in the following. It is clear from Fig. 6 that the density of data points close to the continuous fit line is much higher than the density of the data points that discard from the main trend. The scatter of the data points observed is thus not so significant, in addition to probably be sensitive to some parameter choices made for the data processing through the tessellation method used, such as the choice of the kernel for coarse-graining (see technical details in [28]).

Note that the value $\mu_1 = 0.17$ is smaller than the typical values usually reported in the literature [35–37] and obtained from simulations with an interparticle friction of about 0.5 (as used in the DEM simulations presented here). An even smaller $\mu_1 = 0.12$ was found in a previous study on another nonstandard system, namely the lid-driven granular cavity [28], based on numerical tests with a grain-grain friction of 0.5 too. We do not have any clear physically-sound explanation for such results apart from the fact that those measurements were made on two systems with strong flow gradients for which non-local effects (existence of subyield flows), as reported and modelled by Kamrin and Henann [38], could be significant. The immersed-wall granular system with periodic conditions (studied here) and the lid-driven granular cavity [28] are indeed non-canonical systems that exhibit highly inhomogeneous kinematics in a finite-sized volume: strong variation of the streamlines, quasi-static dead zones, macroscopic vortex dynamics, etc. Note that the upper values of $\mu$ at low $I$ displayed in Fig. 6 match the classical values of about 0.25 for $\mu_1$ that were found under much more homogeneous flow conditions, such as planar shear and annular shear flow [35–37].
<figure><img src="image_7.png" /><figcaption>FIG. 6. The local $\mu(I)$-rheology for this study, extracted from all simulations. The data points come from a region of the granular immersed-wall system selected to approach the conditions of a simple sheared granular flow. That selected region is centered on the system and extends over 90d in the x-direction, by excluding the very bottom (from 0 to 5d) and the very top (from $H - 5d$ to $H$) layers from the analysis. The two regions (of length 30d each) that are close to the vertical side-wall and its periodic image on the left side of the system, and inside which non-local effects may appear (see Fig. 7), were not considered.</figcaption></figure>

Figure 7 displays the spatial maps of $(\mu_{loc} - \mu_{th})/\mu_{th}$ for three values of $I_M$, even if there exist some specific (yet small) regions inside which the principal directions of the strain and stress tensors were not aligned. First, it is generally observed that the granular sample obeys the local $\mu(I)$-rheology in most of the volume of the granular system, except in the regions close to the sidewalls in the two corners at the bottom of the system. Second, at low $I_M$ the patterns resulting from this $(\mu_{loc} - \mu_{th})/\mu_{th}$ mapping is symmetric, while it becomes asymmetric at the highest $I_M$. At the highest $I_M$, the region in the left corner at the bottom, in which the measured $\mu_{loc}$ discards from the $\mu_{th}$ derived from the $\mu(I)$-rheology, grows in size and tends to extend in the middle of the granular system. Those observations are fully consistent with the streamlines shown in Fig. 2. The regions with a dead zone (regardless of $I_M$) correspond to regions where the granular material does not obey the $\mu(I)$-rheology. In
<figure><img src="image_8.png" /><figcaption>FIG. 7. Spatial maps of $(\mu_{loc} - \mu_{th})/\mu_{th}$ for three values of $I_M$: $I_M = 6.1 \ 10^{-3}$ (a), $I_M = 9.1 \ 10^{-2}$ (b) and $I_M = 0.61$ (c). See text for the definition of $\mu_{loc}$ and $\mu_{th}$. The results are shown for $h/d = 20$.</figcaption></figure>

the latter regions, the values of $(\mu_{loc} - \mu_{th})/\mu_{th}$ are generally between 0 and -1 (blue-colored regions that appear in the two bottom corners of each map displayed in Fig. 7) and $\mu_{th} \sim \mu_1$ (small I), thus meaning that the effective local friction $\mu_{loc}$ exactly stands between 0 and $\mu_1$. This is fully consistent with expected creeping regimes (flow below the yield stress), as early reported and captured by the nonlocal granular fluidity model [36, 38].

The region (in the left corner at the bottom), where a vortex-like structure forms at high $I_M$, grows in size with $I_M$ and also corresponds to a region where the granular medium does not obey the local $\mu(I)$-rheology. Again, values of $(\mu_{loc} - \mu_{th})/\mu_{th} \sim (\mu_{loc} - \mu_1)/\mu_1$ are generally found to stand between 0 and -1, thus indicating that a flow occurs below the
yield stress due to nearby flow. We recall here that the data points from the regions close to
the lateral wall where significant non-local effects appear were systematically excluded from
the analysis shown in Fig. 6. Considering those regions in the analysis would increase the
number of data points below $\mu_1$ in the plot of Fig 6.

## IV. FORCE FLUCTUATIONS ON THE IMMERSED WALL

The present section tackles the force fluctuations on the obstacle. The autocorrelation of force time-series (Sec. IV A) and the distributions of grain-wall contact forces (Sec. IV B) are investigated. All data sets obtained for the following analysis were extracted from numerical simulations with an obstacle height $h = 20d$. Keeping the obstacle high enough to capture the collective effects of the grains (in order to avoid discretization issues), while the height available for the grains to overflow was one third of the total height $H$ of the sample, made the presentation of the results clearer. The results of simulations with $h/d = 5, 10$ or $15$ (not shown here) were very similar overall.

### A. Force time-series autocorrelation

In this section, the force fluctuations are investigated by analysing the autocorrelation of force time-series on portions $i$ of the wall. Information about the *memory of the system*, defined hereafter, is then extracted through a critical autocorrelation time called $\Delta t_c^i$.

The force signal record frequency is set to 10 kHz, this value being high enough to capture short-time events so that the results do not depend on this parameter (see more details in [28]). To study the potential heterogeneity of the force signals along the obstacle height, the sidewall is split into five portions $i \in \{1, 2, 3, 4, 5\}$ from bottom to top, of identical size $h/5$. Each slice has its associated force time-series from which the autocorrelation function $C^i(F(t), F(t + \Delta t))$ is then computed, where $\Delta t$ is the lag between two system states. The critical correlation time $\Delta t_c^i$ is extracted from each slice: it can be interpreted as the speed at which the system forgets a certain value of force imparted to the wall. In other words, it reflects the typical time during which a future state of the system keeps history of its past state regarding the force transmitted to the wall. The critical correlation time $\Delta t_c^i$ is defined
with an arbitrary low correlation threshold:

$$C^i(F(t), F(t + \Delta t_c^i)) = 0.15. \qquad (11)$$

As the local autocorrelation function $C^i$ may be a non-monotonic function, it is possible that the threshold value is reached many times. In this case the lowest value of $\Delta t_c$ is kept, avoiding the noise to have an influence on the result.

Figure 8(a) shows the autocorrelation of force time-series for each of the five identical portions of wall, for a low value of the macroscopic inertial number ($I_M = 1.2 \ 10^{-3}$). For this granular regime, the autocorrelations decrease slowly to reach a vanishing value between approximatively 0.2 s to 1.5 s, depending on the wall portion considered. Although an observable trend suggests that the bottom portions are correlated for longer times (lower slope gradient) compared to the top portions (higher slope gradient), no general rule could be extracted because this behaviour is not strictly monotonic.

The values of autocorrelation times $\Delta t_c^i$ for each portions of the wall, as well as their mean over $i$ which we called $\langle \Delta t_c \rangle$, are displayed on Fig. 8. The values of $\Delta t_c^i$ are scattered around $\langle \Delta t_c \rangle$, thus showing that autocorrelation times are inhomogeneous across the height of the sidewall. However, the average value $\langle \Delta t_c \rangle$ shows a remarkable decreasing trend with $I_M$. This behaviour may be interpreted as a decrease of the system memory while increasing the inertial number, which is in accordance with a collisional regime where short-life contacts are dominant compared to long-lasting enduring frictional contacts. Near $I_M = 1.2 \ 10^{-2}$ the mean autocorrelation time nearly vanishes, thus showing that the system is unable to keep memory of a past state. The granular medium, which is at the origin of the force transmission from the boundary conditions to the sidewall, handles here the signature of a remarkable cross-over from quasistatic to dense inertial regime. It is important to note that what we call the *memory of the system* is not an intrinsic memory of the granular material, but rather a system-dependent feature. This becomes obvious when those $\Delta t_c^i$ (and $\langle \Delta t_c \rangle$) are shown to be influenced by the length of the system, as it will be discussed in Sec. V.

## B. Force distributions

In the following, the distributions of force time-series evolving with the boundary macroscopic inertial number are analysed, considering two spatial scales: the entire wall (IV B 1)
<figure><img src="image_9.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b)</figcaption></figure>

FIG. 8. (a): evolution of $C^i(F(t), F(t + \Delta t))$ over time for five identical portions of the wall height $i \in \{1, 2, 3, 4, 5\}$ from bottom to top, and $I_M = 1.2 \times 10^{-3}$. The grey dashed line is the threshold used to obtain $\Delta t_c^i$. (b): critical time $\Delta t_c$ beyond which the force time series become weakly correlated ($\Delta t_c < 0.15$) for the five portions of the wall, as a function of $I_M$. $\langle \Delta t_c \rangle$ is the mean over the five $\Delta t_c^i$ values of a simulation.
and a portion of wall (IV B 2).

## 1. Force distributions on the entire obstacle

This section investigates the distribution of the total force experienced by the sidewall, on the whole surface facing the granular flow. This measurement gives information about the collective effect of the granular bulk on the obstacle. The total force time-series $F$ is obtained by summing all grain-wall contributions at each recording time. This force scaled by the mean force $\overline{F}$ is noted $\tilde{F} = F/\overline{F}$, and its probability distribution density is computed.

Figure 9(a) shows the distributions of $\tilde{F}$, for seven values of $I_M$ ranging from $1.2 \times 10^{-3}$ to 0.6. The distribution shapes follow a nearly Gaussian law at low $I_M$ ($= 1.2 \times 10^{-3}$). A slight asymmetry appears, while increasing $I_M$: the decreasing tail at high forces becomes more and more important. At the same time, from $I_M \ge 1.2 \times 10^{-2}$, the probability to have weak or vanishing force signals on the sidewall increases. A systematic fit of the distributions— actually measured in DEM—using a shifted (and truncated) log-normal function, called $P(\tilde{F})$, provides very good results over the whole range of $I_M$ tested, as depicted by the dashed lines on Fig. 9(a):

$$ P(\tilde{F}) = \frac{1}{S} \frac{1}{\tilde{F}_{\sigma}\sqrt{2\pi}(\tilde{F}+1)} \exp\left(-\frac{[\ln(\tilde{F}+1) - \tilde{F}_{\mu}]^2}{2\tilde{F}_{\sigma}^2}\right), \quad (12) $$

where $\tilde{F}_{\mu}$ is the scale parameter, $\tilde{F}_{\sigma}$ is the shape parameter, and the shift (truncation) is 1. $S$ is the normalization factor corresponding to the value at 1 of the survival function of the (untruncated) log-normal PDF. The evolution of these parameters with $I_M$ is displayed on Fig. 9(b).

Three distinct simulations, done for different macroscopic shear rates and confinement pressures $\{U, P\}$ but keeping the inertial number constant ($I_M = 6.1 \times 10^{-3}$), give identical distributions, as shown on Fig. 9(a). This reveals that the boundary macroscopic inertial number is the key parameter which controls the distribution of the total force on the wall. Such a remarkable result is further confirmed by the monotonic variation of the two fitting parameters $\tilde{f}_{\mu}$ and $\tilde{f}_{\sigma}$ [see Eq. (12)] with $I_M$, as displayed on Fig. 9(b), regardless of $\{U, P\}$. It can be concluded that Eq. (12), fed with the parameters controlled by $I_M$ (Fig. 9(b)), provides an empirical model capable of predicting the total force distribution for the granular flow-wall interaction system studied here.
Moreover, the DEM simulations show that all distributions of $\tilde{F}$ collapse when $I_M \ge 6.1 \ 10^{-2}$, as shown in Fig. 9(a). This observation becomes clearer in Fig. 9(b), as both $f_\mu$ and $f_\sigma$ reach a plateau from $I_M = 6.1 \ 10^{-2}$ up to the highest value of the inertial number achieved in this study. This saturation of the total force distribution is likely to reflect a fast flow state largely influenced by inertial effects. It occurs above a certain value of $I_M$ between $1.2 \ 10^{-2}$ and $6.1 \ 10^{-2}$, which closely matches the value of $I_M$ at which the loss of system memory was identified in Sec. IV A.

In order to analyse the distributions of the force at mesoscopic spatial scale (thus taking into account the potential heterogeneity over the system height), the next section tackles the distributions of force time-series recorded on portions of the obstacle height.

## 2. Force distributions on portions of the obstacle

Similarly to the approach described in Sec. IV A, the sidewall is split into five portions $i \in \{1, 2, 3, 4, 5\}$ of same height each. We note $F_i$ the force time-series associated with the portion $i$ and $\bar{F}_i$ its mean. The upper panel of Fig. 10 shows the distribution of $F_i/\bar{F}_i$, the force arising from each portion scaled by the corresponding mean, for $I_M = 6.1 \ 10^{-2}$. The distributions of the force time-series from the three portions located in the center of the sidewall are equivalent (see orange-colored intermediate lines for $i = 2, 3, 4$), unlike the two located at the top and bottom (red-colored line for $i = 5$ and blue-colored line for $i = 1$, respectively). This difference between the portions $i \in \{2, 3, 4\}$ on the one side, and the portions $i = 1$ and $i = 5$ on the other side, is observed for any $I_M$ (not shown here). This result reveals that the mesoscopic force distributions are homogeneous along a large area located at the center—along the y-axis (vertical) direction—of the obstacle, whereas some modifications near the weak forces appear at the top and the bottom. In the bottom corner, at the interface with the bottom rough wall, few grains can be trapped and exert enduring forces on the obstacle, thus decreasing the probability of weak forces on the time-series. Near the top of the obstacle, the increasing probability of small forces can be explained by the shortening of the dead zone thickness, so that the obstacle directly experiences the force of the flowing grains.

Because the collapse of the distributions for portions $i \in \{2, 3, 4\}$ is observed for all values of $I_M$, the mean of the distributions over these three central portions of the sidewall, which
we name $\langle F_i/\overline{F}_i \rangle_i$, was computed. The variation of $\langle F_i/\overline{F}_i \rangle_i$ with the macroscopic inertial number is shown on Fig. 10(b). Similarly to the distributions on the entire height of the obstacle (Sec. IV B 1), the resulting curves depict log-normal-like shapes. Near $I_M \simeq 10^{-3}$, the distribution displays a small amount of low forces and starts approaching a Gaussian-like shape. The number of weak forces increases with $I_M$, shifting the maximum probability density towards the low forces. This behaviour can be explained by a finite-size effect. At low $I_M$ the force signal may arise from the summation of several enduring contact forces, which is known to produce a Gaussian-like shape. By contrast, at higher $I_M$ the force signal arises from a lower number of contacts, thus increasing the probability of weak (or null) forces. Finally, all the force distributions on wall portions collapse for $I_M \ge 6.1 \cdot 10^{-2}$. The latter observation reveals that—excluding the effects related to the top and bottom ends—the force distributions at mesoscale become independent of $I_M$ as soon as a certain $I_M$ is reached. The collapse of distributions was interpreted as a signature of the transition from quasistatic to dense inertial granular regime, further confirming the result reported in Sec. IV B 1 for the entire height of the obstacle.

## V. EFFECT OF THE SYSTEM LENGTH

The results presented so far were extracted from simulations with a constant length of the system: $L/H = 5$ with $H = 30d$ and $d = 1$ mm. This section compares those results to additional results from simulations where the length of the system was doubled: $L/H = 10$, while keeping $H = 30d$ and $d = 1$ mm.

As already discussed in Sec. III A, the simulations at $L/H = 10$ did not allow producing a plateau around zero in the curves of $V_y/U$ versus $x/L$ (not shown here). The results were similar to the ones shown in Fig. 3. Doubling the length of the system did not allow reducing the effect of the vortex-like structure that is generated by the duplicate obstacle on the left side, and thus disturbing a lot the entire volume of the system for $I_M$ larger than 0.3. This proves to be a signature of the transition towards the collisional regime. Note that in the collisional regime (beyond $I_M \sim 0.3$), the analogy between the periodic system studied here and a free-surface gravity-driven granular flow passing over a rigid wall (which initially motivated the present study) becomes largely questionable.

Figure 11 shows the mean force $\overline{F}$ scaled by the pressure force $Phd$ associated with the
confinement pressure $P$ and the obstacle height $h$, as a function of $I_M^*$ for $L/H = 5$ and $L/H = 10$. The empirical scaling given by Eq. (6) is affected by the ratio $L/H$. Doubling the length of the system yields a slight increase of the scaled mean force, but the empirical scaling that relates $\overline{F}/(Phd)$ to $I_M^*$ is still robust, provided that $r_1$ and $r_2$ in Eq. (6) are considered as slightly dependent on $L/H$.

Figure 12 compares two extreme distributions of grain-wall forces for both $L/H = 5$ and $L/H = 10$ at macro scale ((a): entire height of the obstacle) and at meso scale ((b): portions of the obstacle). At high $I_M$ the results are nearly not affected by varying $L/H$ but at the lowest $I_M$ the distributions become sensitive to $L/H$, particularly when considering the entire height of the obstacle.

Overall, doubling the length of the system does small changes to the main results regarding the time-averaged dynamics and the distributions of grain-wall contact forces. However, moving from $L/H = 5$ to $L/H = 10$ has a significant effect on the variation of $\langle\Delta t_c\rangle$ with $I_M$. While $\langle\Delta t_c\rangle$ is largely positive at small $I_M$ for $L/H = 5$, it tends to approach zero for $L/H = 10$ at the same $I_M$, as shown in Fig. 13. This put emphasis on the fact that the positive value of $\langle\Delta t_c\rangle$ at small $I_M$ must be interpreted as a signature of a memory activated by the system (memory of the system) in terms of the force transmission towards the wall, but not of a memory that would be intrinsic to the granular material itself.

## VI. DISCUSSION AND CONCLUSION

The present contribution described periodic numerical simulations based on DEM of a granular flow passing over a wall-like obstacle, over a wide range of confinement pressure and shearing velocity imposed at the top of the granular flow. This concluding section summarizes the main results concerning the time-averaged dynamics (VIA) and the force fluctuations (VIB). The paper ends by exploring the potential implications of this basic knowledge for the specific problem of the interaction between full-scale granular flows in nature (such as avalanches) and civil engineering structures (Sec. VIC), which was one of the initial motivations of the present study.
A. Time-averaged dynamics

The mean force experienced by a wall (of height $h$) in contact with a dead zone subjected to a granular flow (of incoming height $H$) passing over the wall obeys the following scaling: $\bar{F} = f(I_M^*)Phd$, where $I_M^*$ is a macroscopic inertial number built on the boundary confinement pressure $P$ and the shearing velocity $U$ over the height difference $H - h$, and $f$ is a function which saturated at both low and high $I_M^*$. It was found that the bulk of the granular flow obeyed the local $\mu(I)$-rheology, while the dead zone formed upstream of the wall did not. At very high $I_M^*$, the $\mu(I)$-rheology became, however, much less relevant. A vortex-like structure, inside which the local $\mu(I)$-rheology was not verified, formed downstream of the duplicate obstacle (on the left side of the periodic system simulated) and grew in size, thus producing a disturbance that propagated over the entire volume of the system. This was interpreted as the transition towards the collisional regime.

Figure 14 shows how the mean volume fraction $\Phi_*$ of the system, measured in the DEM simulations, varies with the macroscopic inertial number $I_M^*$. Note that the volume fraction initially measured in two dimensions (spheres in the $(x, y)$ plane) was systematically transformed into an equivalent volume in three dimensions (3D), assuming (roughly) $\Phi_{3D} \approx \frac{2}{3}\Phi_{2D}$ if we compare a sphere of diameter $d$ included in a cube of identical size $d$ to a disk of diameter $d$ included in a square of size $d$. The curve drawn in Fig. 14 has the typical shape of the $\phi(I)$-constitutive law generally extracted from a number of other granular systems, such as plane shear flows [35], inclined flows down a slope [24, 39], or annular shear cell flows [40]. It is worthy to mention that those curves result from simulations at different obstacle heights: all data collapse into one single curve for $h/d$ between 5 and 20 (the case $h/d = 25$ discards a bit from the main trend). The inset of Fig. 14 shows $\Phi_*$ versus $I_M$. A noticeable scattering is observed at high $I_M$, thus confirming that the inertial number $I_M^*$ defined on the shearing height $H-h$ above the obstacle is the relevant parameter for the immersed-wall system studied here. Note that the relation $\Phi_*(I_M^*)$ was found to be slightly dependent on $L/H$ for high $I_M^*$ (not shown here).
## B. Force fluctuations

The analysis of the force fluctuations on the wall-like obstacle was performed through two main measurements. First, the autocorrelations of force time series were computed at a spatial mesoscopic scale, intermediate between the grain size and the entire wall height (considering five portions of the wall, identical in size). Those temporal autocorrelations were found to be inhomogeneous through the wall height. A criterion based on the mean of an autocorrelation threshold through the wall height ($<\Delta t_c>$) was however able to reflect the system memory. Second, the distributions of force signals were computed at both macroscale (wall height) and mesoscale (considering the same portions of the wall as used for the autocorrelations). The boundary macroscopic inertial number $I_M$ imparted to the granular system was found to be the control parameter of both types of measurements, thus highlighting remarkable trends that are summarized below.

For the lowest values of $I_M$ tested, the system exhibited a significant system memory provided that the length of the system was small enough ($L/H$ equal to 5), and narrowed distributions evolving from normal-like to log-normal-like shapes (depending on the spatial scale considered) were observed. This result typically reflects the quasi-static granular flow regime: the wall experiences a constant and large number of enduring contacts with the grains. This increases the time needed for a force value to be forgotten by the system on the one hand, and restricts the range of actual force values on the wall on the other hand. We stress here that this was a system memory which became nearly undetectable in our numerical simulations when the size of the system was increased up to $L/H = 10$.

While increasing $I_M$, the overall memory of the system decreased, and the distributions widened. The inertial forces started to come into play, with faster motions of particles causing shorter contact times of the grains with the wall. In comparison with the quasi-static regime, the broadening of the statistical range of force values is thus interpreted as the consequence of (i) shorter and scarcer contacts with the wall near the weak force values on the one hand, and (ii) more intense interactions near the highest forces on the other hand. The measurements of force fluctuations showed an asymptotic behaviour near $10^{-2} \le I_M \le 6 \ 10^{-2}$, materialized by both the collapse of force distributions and the vanishing of the system memory. In accordance with the well-established time-averaged dynamics of dense granular flows, as reviewed for instance in [39, 41], this asymptotic behaviour was interpreted
as a salient signature of the transition from quasi-static to dense inertial granular regime.

It should be stressed that the results presented in terms of both the distributions of grain-wall contact force and the autocorrelations of force time series are robust because they corroborate the ones we recently obtained on a different system named the granular lid-driven cavity [28]. Indeed, for that different system investigated in [28] with grains being trapped in the cavity and then forced to recirculate inside it (thus exhibiting a large vortex kinematics, different from the one observed in the nonstandard immersed-wall system studied in the present paper), very similar results were found in terms of force fluctuations on the wall facing the direction of the shearing velocity at the top. The distributions of grain-wall contact forces were also entirely controlled by the macroscopic inertial number $I_M$ (regardless of the values of $P$ and $U$ encapsulated in $I_M$), moving from nearly Gaussian-like to asymmetric log-normal distributions when $I_M$ was increased, and this change was also accompanied by a loss of system memory [28]. A detailed comparison between the two systems is beyond the scope of the present study and will be the topic of a future study.

## C. Insights for granular flow-wall interaction in geophysics

A free-surface flow of granular materials down a slope and coming across a perpendicular wall-like obstacle is a common situation in the field of natural hazards related to full-scale gravity-driven granular flows [42]. In particular, protection structures against snow avalanches or landslides are built to limit their propagation through the decrease of their kinetic energy [43]. In this situation, a roughly triangular stagnant zone of grains—a so-called dead zone—can form upstream the obstacle, and a more inertial zone made of flowing grains can overcome the dead zone (see [44] and references therein). One initial motivation of setting up the periodic immersed-wall system investigated in the present work was to gain basic knowledge on the mean force on the wall and the force fluctuations under a wide range of granular flow regimes. Although some results are affected by the length of the system relative to its height, the periodic immersed-wall system developed allowed us to draw a couple of conclusions helpful for the understanding of granular flow-wall interaction and the design of civil engineering structures impacted by natural granular mass flows.

First, the analysis of the mean dynamics suggests that there exists a relatively wide range of flow regimes under which the mean force of a granular flow on a rigid wall of
surface $S_o$ ($S_o = hd$ in the 2D system studied here) should be controlled to first order by the
confinement pressure: $\bar{F}/S_o \sim P$. Second, the ratio of the total time-averaged force $\bar{F}$ on
the wall to the typical pressure force $PS_o$ is not constant: it is a function of the boundary
macroscopic inertial number $I_M^* = \frac{U}{H-h}d\sqrt{\rho/\bar{P}}$. The latter is built on the confinement
pressure $P$ and the typical velocity gradient $U/(H-h)$ of the sheared layer above the wall.
The ratio $\bar{F}/(PS_o)$ typically increases from 1 at low $I_M$ to 3-4 when $I_M$ approaches unity.
The empirical relation between $\bar{F}/(PS_o)$ and $I_M^*$ gives a way to predict the variation of
the force on the wall, before reaching the faster inertial regime for which the mean force
will be controlled by the square of the flow velocity ($F \propto U^2$). Third, the distributions of
grain-wall contact forces (entire height of the wall or portion of wall) are controlled by the
macroscopic inertial number and are much larger at high $I_M$ than at low $I_M$. This result
can be summarized by introducing the coefficient of variation $\sigma_F/\bar{F}$ (where $\sigma_F$ denotes the
standard deviation) and considering its variation with $I_M$, as shown on Fig. 15.

The relative amplitude of the force fluctuations is about 25% in the quasi-static regime
and increases logarithmically with the inertial number. In the collisional regime, the relative
amplitude of force fluctuations reaches a constant value (independent of $I_M$) which becomes
very significant, about 80% for $h/d = 20$ in our study. In the dense inertial regime ($0.25 <
\sigma_F/\bar{F} < 0.7 - 0.8$) and the collisional regime ($\sigma_F/\bar{F} \approx 0.7 - 0.8$), the fluctuations play
a major role in the magnitude of the total force and are thereby a key ingredient for the
design of civil engineering structures (protection dams, buildings, etc.) likely to be impacted
by granular mass flows. It is worthwhile to remind that the results were obtained in two-
dimensional configurations (spheres stay in the $(x, y)$ plane on Fig. 1, by construction).
Three-dimensional simulations will need to be undertaken in the future in order to check
some factors that are likely to influence the force distributions. These may include the effects
caused by a drastic increase of the number of grains in contact with the wall, or the increase
of steric hindrance in the granular bulk.

673 ACKNOWLEDGMENTS

This work has been partially supported by the LabEx Tec21 (Investissements d'Avenir: grant agreement No. ANR11-LABX-0030). Thierry Faug and Mohamed Naaim are grateful to the financial support by the People Programme (Marie Curie Actions) of the EU 7th FP
under REA Grant Agreement No. 622899 (FP7-PEOPLE-2013-IOF, GRAINPACT).

[1] K. Wieghardt, Annual Review of Fluid Mechanics **7**, 89 (1975).

[2] R. Albert, M. Pfeifer, A. Barabasi, and P. Schiffer, Physical Review Letters **82**, 205 (1999).

[3] D. Chehata, R. Zenit, and C. Wassgren, Physics of Fluids **15**, 1622 (2003).

[4] T. A. Brzinski III and D. J. Durian, Soft Matter **6**, 3038 (2010).

[5] A. Seguin, Y. Bertho, F. Martinez, J. Crassous, and P. Gondret, Phys. Rev. E **87**, 012201 (2013).

[6] R. Bharadwaj, R. Zenit, and C. Wassgren, Physics of Fluids **18**, 043301 (2006).

[7] C. Wassgren, J. Cordova, R. Zenit, and A. Karion, Physics of Fluids **15**, 3318 (2003).

[8] I. Albert, P. Tegzes, R. Albert, J. Sample, A. Barabasi, T. Vicsek, B. Kahng, and P. Schiffer, Physical Review E **64**, 031307 (2001).

[9] J. Geng and R. Behringer, Physical Review E **71**, 011302 (2005).

[10] E. I. Corwin, E. T. Hoke, H. M. Jaeger, and S. R. Nagel, Phys. Rev. E **77**, 061308 (2008).

[11] E. Longhi, N. Easwar, and N. Menon, Phys. Rev. Lett. **89**, 045501 (2002).

[12] E. Gardel, E. Sitaridou, K. Facto, E. Keene, K. Hattam, N. Easwar, and N. Menon, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences **367**, 5109 (2009).

[13] D. M. Mueth, H. M. Jaeger, and S. R. Nagel, Phys. Rev. E **57**, 3164 (1998).

[14] F. Radjai, S. Roux, and J. J. Moreau, Chaos **9**, 544 (1999).

[15] J. H. Snoeijer, M. van Hecke, E. Somfai, and W. van Saarloos, Phys. Rev. E **67**, 030302 (2003).

[16] J. H. Snoeijer, T. J. H. Vlugt, M. van Hecke, and W. van Saarloos, Phys. Rev. Lett. **92**, 054302 (2004).

[17] A. Petri, A. Baldassarri, F. Dalton, G. Pontuale, L. Pietronero, and S. Zapperi, The European Physical Journal B **64**, 531 (2008).

[18] C. h. Liu, S. R. Nagel, D. A. Schecter, S. N. Coppersmith, S. Majumdar, O. Narayan, and T. A. Witten, Science **269**, 513 (1995), http://science.sciencemag.org/content/269/5223/513.full.pdf.

[19] F. Radjai, M. Jean, J.-J. Moreau, and S. Roux, Phys. Rev. Lett. **77**, 274 (1996).
[20] B. Miller, C. O'Hern, and R. P. Behringer, Phys. Rev. Lett. **77**, 3110 (1996).

[21] F. Radjai, Comptes Rendus Physique **16**, 3 (2015), granular physics / Physique des milieux granulaires.

[22] D. Howell, R. P. Behringer, and C. Veje, Phys. Rev. Lett. **82**, 5241 (1999).

[23] R. Soller and S. A. Koehler, EPL (Europhysics Letters) **80**, 14004 (2007).

[24] T. Faug, R. Beguin, and B. Chanut, Physical Review E **80**, 021305 (2009).

[25] B. Chanut, T. Faug, and M. Naaim, Physical Review E **82**, 041302 (2010).

[26] T. Faug, P. Caccamo, and B. Chanut, Physical Review E **84**, 051301 (2011).

[27] V. Šmilauer, A. Gladky, J. Kozicki, C. Modenese, and J. Stránský, in *Yade Documentation*, edited by V. Šmilauer (The Yade Project, 2010) 1st ed., http://yade-dem.org/doc/.

[28] F. Kneib, T. Faug, G. Nicolet, N. Eckert, M. Naaim, and F. Dufour, Phys. Rev. E **96**, 042906 (2017).

[29] F. Kneib, T. Faug, F. Dufour, and M. Naaim, Computational Particle Mechanics **3**, 293 (2016).

[30] P. Cundall and O. Strack, Geotechnique **29**, 47 (1979).

[31] T. Schwager and T. Pöschel, Granular Matter **9**, 465 (2007).

[32] S. Luding, Granular Matter **10**, 235 (2008).

[33] P. Jop, O. Pouliquen, and Y. Forterre, Nature **441**, 727 (2006).

[34] GDR-MiDi, European Physical Journal E **14**, 341 (2004).

[35] F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Physical Review E **72**, 021309 (2005).

[36] K. Kamrin and G. Koval, Phys. Rev. Lett. **108**, 178301 (2012).

[37] K. Kamrin and G. Koval, Computational Particle Mechanics **1**, 169 (2014).

[38] K. Kamrin and D. L. Henann, Soft Matter **11**, 179 (2015).

[39] Y. Forterre and O. Pouliquen, Annual Review of Fluid Mechanics **40**, 1 (2008).

[40] A. Fall, G. Ovarlez, D. Hautemayou, C. Mézière, J.-N. Roux, and F. Chevoir, Journal of Rheology **59**, 1065 (2015).

[41] P. Jop, Comptes Rendus Physique **16**, 62 (2015), granular physics / Physique des milieux granulaires.

[42] T. Faug, Phys. Rev. E **92**, 062310 (2015).
736 [43] T. Faug, P. Gauer, K. Lied, and M. Naaim, Journal of Geophysical Research 113, F03009 (2008).

738 [44] T. Faug, P. Caccamo, and B. Chanut, Geophysical Research Letters 39, L23401 (2012).
<figure><img src="image_11.png" /><figcaption>FIG. 9. (a): distribution of the total force on the right sidewall for seven values of $I_M$ ranging from $1.2 \times 10^{-3}$ to 0.6. The dashed lines show fits using a (truncated) log-normal PDF given by Eq. (12).
(b): variation of $f_\mu$ and $f_\sigma$, the two parameters of Eq. (12), as function of $I_M$.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 10. (a): distribution of force time-series experienced by five equal portions *i* of the sidewall for *N*<sub>U</sub> = 10 and *N*<sub>P</sub> = 1 (each portion time-series is scaled by its corresponding time-averaged force). (b): mean distribution over *i* for *i* ∈ {2, 3, 4} for seven values of *I*<sub>M</sub>.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIG. 11. Scaled mean force $\bar{F}/(Phd)$ versus $I_M^*$ for $L/H = 5$ and $L/H = 10$. For each $L/H$, four values of $h/d$ were considered: 5, 10, 15 and 20. The lines show the predictions of Eq. (6) with $r_1 = 1.4$ and $r_2 = 3.3$ (solid line) and with $r_1 = 1.5$ and $r_2 = 4.1$ (dashed line). For both predictions $I_M^{*0} = 0.3$ is used.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 12. Distributions of the total force on the entire height of the obstacle (a), and mean over $i$ (for $i \in \{2, 3, 4\}$) of the distributions on portions of obstacle (b), for $L/H = 5$ (thick lines) and $L/H = 10$ (thin lines). In each panel, we present the results for $I_M = 1.2 \cdot 10^{-3}$ (the lowest $I_M$ tested) and $I_M = 6.1 \cdot 10^{-2}$ ($I_M$ beyond which the distributions collapse).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIG. 13. Critical times $\Delta t_c^i$ (and their mean $\langle \Delta t_c \rangle$) beyond which the force time series become weakly correlated ($\Delta t_c < 0.15$), for $L/H = 5$ and $L/H = 10$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIG. 14. Volume fraction $\Phi_{*}$ actually measured in the DEM simulations as a function of $I_{M}^{*}$. Inset: $\Phi_{*}$ versus $I_{M}$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 15. Coefficient of variation $\sigma_{F}/\bar{F}$ as a function of $I_{M}$ for two values of $h/d$.</figcaption></figure>