Analytical nonlocal model for shear localization in wall-bounded dense granular flow

This work employs a Landau-Ginzburg-type nonlocal rheology model to account for shear localization in a wall-bounded dense granular flow. The configuration is a 3D shear cell in which the bottom bumpy wall moves at a constant speed, while a load pressure is applied at the top bumpy wall, with flat but frictional lateral walls. At a fixed pressure, shear zones transit from the top to the bottom when increasing lateral wall friction coefficient. With a quasi-2D model simplification, asymptotic solutions for fluidization order parameters near the top and bottom boundaries are sought separately. Both solutions are the Airy function in terms of a depth coordinate scaled by a characteristic length which measures the width of the corresponding shear zone. The theoretical predictions for the shear zone widths against lateral wall friction coefficient and load pressure agree well with data extracted from particle-based simulation for the flow.


Introduction
Granular flow is ubiquitous in many industrial applications (grain transport and storage) and geophysical events (landslides and debris flows), yet a unified description remains a challenge. Unlike usual fluids, granular material yields and flows when the ratio of bulk shear to normal stress exceeds a certain critical value [1]. Consequently, the material deforms nonuniformly by means of shear bands near boundaries [1][2][3][4][5]12,13]. Such a shear localization phenomenon is often described by a constant yield bulk friction for the flow to occur [1,3]. However, many studies have shown that solid boundary can have a strong influence on shear localization [2,4,5,12,13]. It is interpreted as nonlocal effects that wall disturbances can fluidize regions below yield through long-range transport along grain-scale contact network [5][6][7][8][9][10][11]. Several nonlocal models with diffusion order parameter have been proposed to capture wall shear localization in 2D confined flows [6][7][8][9][10][11].
Recently, Artoni and Richard show that shear localization can occur at several locations along the depth in a 3D rectangular parallelepiped periodic cell flow [12] or in a 3D torsional shear cell flow [13]. They found that the location of shear localization strongly depends on lateral wall friction coefficient. While a simplified stress analysis can provide a basic understanding of the effect of lateral wall friction on shear localization patterns [14], a full description for the flow rheology is desired.
In this paper, a Landau-Ginzburg-type nonlocal rheology model is used to describe shear localization

Shear localization transition
The configuration of the wall-bounded granular flow in [12] is shown in Fig.1(a), which is a rectangular channel with a periodic boundary condition in the main flow direction in x, with two bumpy walls at the top and at the bottom, and two lateral flat but frictional walls (in the ydirection). Gravity acts on the system along the zdirection. The flow is driven by applying a wall speed V at the bottom along the x-direction. The top wall is only free to move in the z-direction according to the pressure applied on it. By dimensional analysis, the flow is controlled by three parameters: (1) particle Froude number, Fr≡V/(gd) 1/2 (2) loading parameter M≡P0/ρgH, which compares the loading pressure on the top wall P0 and a pressure due to total particle weight (3) particlelateral wall friction coefficient μpw. Here, ρ denotes the bulk density, g the acceleration of gravity, d the particle diameter, and H the material height. In [12], profiles along z are computed by performing slice-averaging with a thickness 2d in the z-direction. The flow behavior is found insensitive to the Froude number but strongly depend on the load parameter and the lateral wall friction, which can be categorized into three regimes according to the depth velocity profiles shown in Fig.1 Fig.1(b), where σzx b is a bulk deviatoric shear stress on the z-face acting in the x-direction, and P≡(σxx b +σyy b +σzz b )/3 is a bulk isotropic pressure. It displays that increasing μpw weakens μ near the top but enhances μ near the bottom. Such a transition of the bulk friction gradient is speculated to cause the flow regime transition and has been qualitatively captured by a simplified stress analysis in [14] based on force balances between gravity and lateral wall friction. In [13], an inertial μ(I)-rheology model [1,3] is numerically solved to predict the detailed flow profiles in the shear cell while the predictions show rate-dependent depth velocity profiles, inconsistent with the rate-independent profiles observed in the simulations. Indeed, the rateindependence characterizes quasistatic granular rheology and suggests that the finite shear zones in the three regimes are generated by nonlocal processes driven by wall shearing along grain contact network. While several nonlocal approaches have been put forward to account for similar long-range transport phenomenon [6][7][8][9][10][11]15], here we adopt a simple yet effective Landau-Ginzburg approach which introduces a diffusion-reaction equation describing shear-induced fluidization processes, which has been shown to well capture the nonlocal rheology of 2D free-surface inclined flows and confined sheared flows [6,10,11].

Model formulation 3.1 Landau-Ginzburg order-parameter model
A phenomenological Landau-Ginzburg order parameter model is introduced and simplified for the current flow condition as follows. The basic idea of the model is to treat the transition of granular non-sheared and sheared states as a solid-fluid phase transition. In this sense, thermal temperature in the original Landau-Ginzburg framework is replaced with the effective friction coefficient μ which controls the occurrence of granular "phase" change due to material yield. Hence, a Landau-Ginzburg differential equation describing how a shearinduced fluidization state evolves in space and time under a given μ-field can be written as Here, λ denotes an order parameter (OP) measuring the degree of shear-induced fluidization at a location r and time t, t0 represents a relaxation timescale, l a microscopic lengthscale (~d). The source term corresponds to the derivative of a coarse-grain Landau free energy density, which can be expanded to be where μc denotes a critical friction coefficient at which the material yields, and α and β are positive exponents with β>α>0. Note that different forms of Floc have been proposed in literature [6,10,11]. For a steady, quasistatic flow considered here, the unsteady term and the higher-order source terms in Eq.(2) are dropped to get This simplified equation describes how a steady flow zone above yielding μ>μc fluidize neighboring zones below yielding (μ<μc). The nonlocal Laplacian term in Eq.(3) hence extends the classical Mohr-Coulomb plasticity theory which indicates zero thickness for shear zone. Note that in order to predict flow profiles the OP has been linked to a variety of rheological parameters such as the fraction of fluid stress contribution [6], a granular fluidity, the inverse of viscosity scaled by pressure [8,10], and the inertial number I, a scaled shear rate by pressure indicating the quasistatic and the inertial regimes [9,11]. Also note that the form of Eq. (3) is identical to the steady-state-only granular fluidity theory in [8] without the inertial-rheology contribution. Yet the effectiveness of these OP definitions in the current geometry needs to be further investigated via systematically comparing the predicted profiles to the flow and the stress data from the simulations. As a preliminary study here, we aim to employ the OP model (3) to derive the scaling lengths of the shear zones and examine the predictions using the width data of the shear zones extracted by exponential fit to the depth velocity profiles in simulation.

Width-averaged stress model
To solve the nonlocal model (3), the effective friction profile is required. To simplify the problem, the widthaveraged technique is employed to reduce the momentum balance equations. This treatment is reasonable as long as the flow height, H, is far greater than the channel width, W. In our simulations, the ratio K≡H/W=4.5 which is acceptable. The width-averaged momentum equations in x and y are given by where σyx w is a deviatoric shear stress on the lateral walls acting in the x-direction, σzx deviatoric shear stress on the z-face acting in the x-direction, and σzz normal stress in the z-direction, and <…> denotes the width-averaging in the y-direction. To avoid confusion with the lateral walls, z=H and z=0 are referred to as the top and bottom boundaries hereinafter. The averaged normal stress <σzz> is assumed isotropic and is prescribed to be <σzz>(H)=P0 at the top boundary. The effective friction and the lateral wall friction coefficient are defined as respectively. For the sake of simplicity, here we assume the lateral wall friction is equal to the particle-wall friction, i.e., μw=μpw. Note however that this assumption may be inappropriate for large μpw where significant wall friction weakening, i.e., μw<μpw, is observed everywhere in the system as shown in [12]. With this assumption, Eq.(4a,b) and (5a,b) can be solved to be

Near-boundary asymptotic analysis
The asymptotic solutions near the top and bottom boundaries are sought. The boundary conditions for the OP are prescribed to be no-flux at the two solid boundaries, meaning that there is no flow of shearinduced fluidization across the boundaries. In addition, the OP is assumed to vanish as approaching the central zone. Similar OP boundary conditions have been adopted in 2D nonlocal modelling for other confined flows [16]. Assuming the widths of the shear zones (scaled with l) are asymptotically smaller than the flow height, Eq. (3) and (6)  Here, ΔμT and ΔμB measure the significances of the spatial nonuniformity of μ in the shear zones at the top and bottom, respectively: Here, qc=-1.019 is a numeric value at the first local maximum of Ai(q) along the negative q axis [Ai'(qc)=0], and λT and λB denote the OP at the boundaries. Significantly, Eq.(11) and (12) indicate two lengths associated with the widths of the boundary shear zones: , .
In order to check these expressions, we performed molecular dynamics simulations of the confined shear flow with the LIGGGHTS software [17], for 10000 spheres (corresponding roughly to H=45d), one value of the Froude number (Fr=1), and different values of the dimensionless load (M=0.1, 1, 10) and the particle-wall friction coefficient (μpw =0-0.4). A linear spring-dashpot contact model with friction was used, with normal stiffness kn = 8×10 5 mg/d , tangential stiffness kt = 2/7 kn, restitution coefficients en =0.7, et =0, and interparticle friction μp = 0.5. From these simulations, we estimated the top and bottom shear zone widths as the exponential decay lengths obtained by exponential fit to the velocity data. Figure.

Conclusion
We have analytically solved a Landau-Ginzburg nonlocal rheology model to account for the shear localization in a 3D wall-bounded sheared granular flow shown in Fig.1. With the steady quasistatic flow condition and the width-averaged effective friction profiles, the boundary asymptotic solutions for the fluidization order parameters are solved to be the Airy function of the first kind. The solutions reveal two decay lengths associated with the shear zone widths at the top and bottom. We have shown that the width solutions agree well with the data extracted from the particlebased simulations and can account for the observed flow regime transition.
A key conclusion draw from the present study is that the formation of the shear zones strongly depends on the extent of nonuniformity in the effective friction coefficient μ near the solid boundaries, as shown in Eq. (15) and (16). According to the nonlocal model (3), if μ is significantly nonuniform near a boundary (large ΔμT or ΔμB), it allows the coexistence of a near-boundary flow region above yielding (μ>μc) and a far-field creep region below yielding (μ<μc). On the contrary, when μ becomes uniform near the boundary, μ tends the yielding point μc in the two regions so that the material near the boundary becomes solidlike and hence the shear zone vanishes. For the wall-bounded flow considered here, increasing the lateral wall friction and the loading weight result in a locally uniform μ profile near the top (smaller ΔμT) but enhances the gradient of μ near the bottom (higher ΔμB), causing the shear zone to transit from the top to the bottom [see Fig.1(b)]. The current nonlocal model is consistent with the physical arguments based on the stress analysis in [14] and further captures the rheological features beyond the scope of the inertial μ(I) rheology [13].