Compaction of granular material inside confined geometries

In both nature and engineering, loosely packed granular materials are often compacted inside confined geometries. Here, we explore such behaviour in a quasi-two dimensional geometry, where parallel rigid walls provide the confinement. We use the discrete element method to investigate the stress distribution developed within the granular packing as a result of compaction due to the displacement of a rigid piston. We observe that the stress within the packing increases exponentially with the length of accumulated grains, and show an extension to current analytic models which fits the measured stress. The micromechanical behaviour is studied for a range of system parameters, and the limitations of existing analytic models are described. In particular, we show the smallest sized systems which can be treated using existing models. Additionally, the effects of increasing piston rate, and variations of the initial packing fraction, are described.


Introduction
When granular materials are placed in confined geometries, we often observe a significant portion of the stress being redirected towards the confining boundaries. This phenomenon has been systematically studied for many systems [2,10,16], most notably in silos, beginning with [8]. Force redirection has been attributed to the granular nature of the material, and has in many cases been shown to be well represented by a constant coefficient, known as the Janssen coefficient K, defined in one spatial dimension as where σ r is the redirected stress due to some applied normal stress σ n . Here we investigate the development of stresses within a granular packing, confined between two horizontal plates, subjected to a rigid piston impacting it from one side. As the piston moves, granular material is compacted near the piston, and with increasing displacement of the piston, the size of the packing increases. Such an accumulation process is known to occur in the petroleum industry, where sand is liberated from the host rock during extraction, altering the underground morphology of cracks [18,19]. This may also be relevant for understanding proppant flowback in propped fractures [12]. Additionally, this geometry is representative of a number of recent experimental studies in Hele-Shaw cells [4,9,14,15] where the validity of Janssen stress redirection has not been ascertained.
There are a number of interesting patterns which form when a granular material is displaced by a flexible interface in such a geometry [9,15]. The nature of the patterns have been shown to depend on many factors, primarily the initial packing fraction and the rate of displacement [14]. For this reason, we here investigate the microstructural and mechanical evolution of such a system under these conditions. To reduce the complexity of the system, we consider only a rigid piston.
We are interested in systems which are highly confined. In common experiments with granular material inside Hele-Shaw cells, there are in general fewer than 20 grain diameters between the two Hele-Shaw plates, typically down to around 5 grain diameters [14]. As the confinement increases, i.e. as the gap spacing decreases, we expect a transition from three dimensional behaviour towards a behaviour governed by the boundaries, as demonstrated for vertical silos in [1]. It is then of interest to study the changes that result from increasing confinement. We expect that altering the confinement will affect the force redistribution. A transition may occur for extremely confined systems where such an assumption concerning force redirection may not be valid.
Existing analytic models for the micromechanics of such a system generally reduce the problem to one spatial dimension (x), assuming that variations in both remaining directions (y and z) are small, although recently curved interfaces have also been described [4]. For simplicity, we consider a flat interface, and validate the analytic description with a discrete element model.
Firstly, in Section 2 describe the numerical model that has been used to simulate this system. In Section 3, we establish continuum properties which correspond to the analytic formulation, and show comparisons between the two. In particular, the limitations of current analytic models are identified. Finally, a parameter set is proposed that best fits the analytic theories for a wide range of system variables.

Numerical model
This paper is an investigation into the micromechanics of a system which is highly constrained by external boundaries. For this reason, it is ideal to use a particle based method to model the behaviour, as the total number of grains in the system is small. Towards this end, we use a conventional soft sphere discrete particle approach, implemented in the open source code MercuryDPM Figure 1. Particle positions during a single test, shown at six piston displacements s. Top to bottom: s = 0, 10, 20, 30, 40 and 50. Labels refer to the coordinate system x, y and z, the gap height D, the plug length L and the initial packing fraction φ 0 . Particles are coloured by size, darker colours representing smaller particles.
(www.mercurydpm.org) [17,21]. The geometry considered here, shown in Figure  1, consists of a rigid piston, oriented in a space r = {x, y, z}, with normal along the x axis, which pushes particles between two rigid walls, separated by a spacing D and having normals in the ±z directions, with two periodic boundaries in the remaining perpendicular direction, y. The coordinate system moves with the piston, such that it is located at x = 0 at all times. As the piston moves horizontally at a velocity u towards the grains, its displacement at any time t is then s = ut.
We work in a system of non-dimensionalised units with the following properties; length and mass have been non-dimensionalised by the length d m and mass m m of the largest particle in the system, respectively, where the prime indicates that the quantity has dimension. The particle diameters, d, we use are therefore d ≤ 1, with material density defined by 4π(1/2) 3 ρ p /3 = 1, or ρ p = 6/π. Time is non-dimensionalised by the time taken for the largest particle to fall from rest its own radius under the action of gravity, so that a unit time is t = d m /g, which requires that g = 1. Other values are non-dimensionalised by a combination of these three scales, for example stress is non-dimensionalised by m m g /d  Figure 2. Contact laws used in the discrete element model. Interactions are characterised by (a) normal, (b) tangential and (c) rolling laws, parameterised by stiffnesses k, k t and k r , viscous dissipation γ, γ t and γ r and friction coefficients µ t and µ r . Full details are given in [11]. As shown in Figure 2, the particles' material properties are described by normal, tangential and rolling damped spring sliders [11,17] with the properties contained in Table 1. These values have previously been calibrated to mimic micron sized silica beads [5]. The walls are implemented such that they are rough; when a particle contacts a wall, the piston, or both, it is prohibited from rotating, i.e. µ r = 1. Otherwise, the interaction properties are the same as between two particles, except that the walls and piston are of infinite mass. We therefore have a well defined macroscopic sliding friction of µ = 0.4 that does not depend on the rate of loading.
In the following Section we will firstly detail the important macroscopic quantities measured from a single simulation. We will then investigate the effect of three controlling parameters on the evolution of the system: the Hele-Shaw spacing D, the initial packing fraction φ 0 and the velocity of the piston, u. The piston rate, however, is not a priori a governing quantity, so we choose to control the piston rate via the inertial number, I, which is defined as I =γd m / P/ρ p , which is the ratio of inertial to imposed stresses, where P is a typical pressure andγ is a typical shear strain rate [6]. Takingγ = u/D, and P = ρ p gD, gives

Results
As depicted in Figure 1, three distinct regions exist inside the system. These are termed the plug zone, where particles are densely packed near to the piston (x ≤ L), the undisturbed zone, far from the plug, and the transition zone, where the plug is accumulating. To define these regions systematically, we must first measure the solid fraction, ν = V s /V t , which is a local measure of the ratio of the volume of solids to the total volume. The solid fraction is coarse grained in one spatial dimension, x, as where N is the number of grains in the system, W is the width of the system, V i is the volume of the i-th grain, x i its centre of mass and W is the coarse graining function, chosen in this case to be a 1D normalised Gaussian function [7,21]. Such a coarse graining method allows the coarse graining width to be defined, so that either the macro-or micro-structure is visible. We choose to use a coarse graining width equal to the maximum particle diameter, such that small scale variations are minimised, and smooth continuum fields are obtained [20]. All other continuum quantities are defined using the same coarse graining technique, presented in [7,21]. Using the definition (1), we denote the maximum solid fraction, ν m , as an average over the solid fraction close to the wall at some time when the transition zone is far from the piston as where in practice a numerical integration is done over the discrete coarse grained cells, and the limits of 5D and 10D are chosen arbitrarily. We define the normalised packing fraction, φ(x, t), as φ = ν/ν m , and note that in the absence of volumetric expansion or dilation of the granular material, φ is directly proportional to the height of the packing between the Hele-Shaw walls. The undisturbed zone is that which is maintained at the initial packing fraction φ 0 , which is defined at time t = 0 as To delineate the plug, transition and undisturbed zones, at each time t we use linear regression to find the best linear fit to the points in the range φ = φ 0 + 0.1 to φ = 0.9, such that φD ≈ b − x tan θ, for some value of b and slope angle θ. Five examples of this are shown in Figure 3A. Using this best fit, we can define the point at which the plug zone meets the transition zone as the intersection of the best fit regression line and φ = 1, such that L = (b − D)/ tan θ, as shown in Figure 3A. The value of L grows monotonically with piston displacement, as shown in Figure 3B. Conservation of mass implies that on average, if there is no volumetric change in the packing, and no slip between the grains and the walls, this relationship can be expressed as which is shown as the dashed line in Figure 3B. For all cases reported here, the value of θ does not appear to vary with increasing plug length L. The point at which the undisturbed zone meets the transition zone can then be defined in a similar manner as above, by the intersection of the best fit regression line with φ = φ 0 . The coordination number, Z, (average number of contacts per particle), is fairly constant in the plug zone, (Figure 3C), increasing with compaction, as rearrangement occurs. Additionally, at large values of L the stresses are high enough to cause significant overlap of the particles (up to 1%). In the transition zone, significant particle rearrangement lowers the coordination number. Coarse graining techniques in general cause measured fields to converge towards zero near boundaries [21]. While it is feasible to reconstruct these fields in general near individual boundaries, near the piston we have three distinct boundaries which all interact. To access stresses in this region, it is then preferable to directly measure the forces applied to the rigid boundaries of the system. Towards this end, we denote σ p n as the normal stress measured at the piston. This is shown as a function of piston displacement in Figure 3D. The stresses measured from coarse graining within the packing are shown in Figure 3E. In both cases, the stresses grow exponentially both with increasing L and decreasing x, as shown by the linearity in a semilog space in Figure 3D and E.
An analytic expression to describe this stress evolution was first derived in [9], assuming that: (a) the stress redirection in the z direction is described by a constant K z = σ zz /σ xx , where σ zz = σ zz − Dρ p ν m g/2 is the component of the vertical stress not due to gravity and (b) friction at the side walls is µ = σ xy /σ zz . It can be shown using force balance in the x direction that under these conditions, if there is no net acceleration, By further assuming that the stress at x = L is a constant, i.e. σ T ≡ σ xx (x = L), Previously, the threshold stress σ T has been modelled as either estimated from experimental data [4], or as a constant by assuming sliding of a wedge of material [9,15]. Here, we choose to model the threshold stress σ T as a φ 0 dependent quantity by considering limit equilibrium of a wedge of material being displaced into the undisturbed zone, as shown in Figure 4. We assume a noncohesive Coulomb failure of the material internally, along a failure plane parallel to and meeting the surface of the transition zone. We additionally assume that the internal friction angle is also defined by µ. Limit equilibrium in the x direction can then be expressed using the notation defined in Figure 4 as F x = T b + T i cos θ, where F x = Dσ T , T b = µW 0 = µD 2 ρ p ν m g/(2 tan θ) and T i = µW 1 cos θ = µD 2 ρ p ν m gφ 2 0 cos θ/(2 tan θ) per unit length in the y direction. After some rearrangement this implies that This assumption of the failure surface introduces no new parameters into the model, and as will be shown in the following, closely predicts the measured value of σ T for a large range of system parameters. In the limit where φ 0 → 0, this definition reduces to that used in [9] and [15]. Including this new definition of the threshold stress, Equation (5), in Equation (4) gives A best fit estimate is used to find K z from (6), and is shown as the dashed line in Figure 3D at x = 0, using the measured values of ν m , θ and φ 0 , which adequately captures the behaviour of the system past L/D = 2. Before this limit, stress redistribution has not saturated, and σ xx is less than the predicted value. We find the same transition value of L/D ≈ 2 for all cases reported here.
The measured value of apparent friction µ = σ xz /σ zz inside the packing, shown in Figure 3F shows that the system is far from failure inside the plug zone, and increasingly unstable in the transition zone. The out of plane Janssen coefficient, K y = σ yy /σ xx , shown in Figure 3G, has large variations in x, but is not significantly affected by the formation of the plug. The in plane Janssen coefficient ( Figure 3H), K z = σ zz /σ xx , however, is strongly influenced by the formation of the plug, and is relatively constant with increasing L inside the plug zone.   An underlying assumption of the Janssen stress redistribution is that when averaging over the width of the system (here in the y and z directions), the principal stress directions are parallel to the system geometry [8]. For this reason we measure α, the absolute value of the x component of the eigenvector of the major principal stress, which is shown in Figure 3I. When α ≈ 1 the major principal stress points along the x-axis, and when α ≈ 0, the principal stress lies in the yz plane. For x ≤ L we find that the major principal stress is collinear with the system geometry, and the Janssen stress model fits well. In a traditional silo problem, gravity acts parallel to the applied stress, and averaging across the width of a silo ensures the validity of this assumption. For this case, however, because the direction of gravity has broken the inherent symmetry of the silo problem, its validity is not ensured [13]. We do, however, find that this assumption holds well for this and all simulations reported here.

Gap spacing
As motivated in Equation (6), the gap spacing D largely controls the magnitude of the stresses within the system. For this reason, we here vary this spacing systematically from D = 1 to D = 15, while maintaining φ 0 = 0.5 ± 0.05 and I = 0.01 (except for the case of D = 1, where φ 0 ≈ 0.66), as depicted for four values of D in Figure 5. To make a reasonable comparison between these simulations, in each case the area of the piston is kept constant, such that its area is W D = 100. Select measures of the behaviour of the system are shown in Figure 6. Slope angles, θ, are averaged over the times corresponding to 2D ≤ L ≤ 10D, whilst K and σ T are averaged temporally over values in the range 2D ≤ L ≤ 10D, where at each time we measure spatially in the range L/4 ≤ x ≤ 3L/4. In Figure 6A and B, we observe increasing solid fraction in the plug, ν m , and slope angle, θ, with increasing gap spacing, as the effect of the boundaries on the system decreases. For D < 8, we note that the effect of the For each simulation, the measured normal stress at the piston, σ p n , is fitted with Equation (6), and a best fit estimate of K z is shown as crosses, with the standard deviation of the error of the regression used as error bars, in Figure  6C. The mean and standard deviation of the measured values of K z and K y from the continuum data between L = 2D and L = 10D are shown as dots and triangles, respectively. For D ≥ 2 we find that both the measured values of K z and K y are independent of D, and have mean values of K y = 0.67 ± 0.04 and K z = 0.58±0.05. Best fit estimates of K z are also independent of D, with mean values of K z = 0.56 ± 0.07. Figure 6D shows the mean of σ T also from L = 2 to L = 10. Triangles represent the prediction from Equation 5 using the measured values of ν m , θ and φ 0 . In all cases, we find the measured and fitted values of K z to be in agreement, whereas the values of σ T agree only with D ≥ 3. We note, however, that σ T depends strongly on θ, and we have as yet no means for predicting this quantity. The dependence of σ T on θ is in contrast to studies on fold and thrust belts [3], where there is no confinement vertically above the material.

Initial packing fraction
Existing models [4,9] for the evolution of the system have neglected any effect of the initial packing fraction φ 0 . To test this assumption, we here vary φ 0 from 0.1 to 0.6, while maintaining D = 5 and I = 0.01. As shown in Figure  7A and B, the packing fraction inside the plug and the slope of the transition zone are independent of the initial packing fraction. In Figure 7C, we observe that the measured and best fit values of K z are in agreement for a wide range of φ 0 , and that these values are lower than the measured values of K y . The prediction of threshold stress from Equation (5), which includes a dependence on φ 0 , slightly under-predicts the threshold stress at φ 0 = 0.1. Nevertheless, both the measured and predicted values of the threshold stress in Figure 7D are in agreement with observations from [4], where a non-dimensional threshold (C) Slope of the best fit value of L(s) denoted by black dots against prediction using incompressibilty shown as the shaded region, which denotes ± one standard deviation around the mean value for each case. (D) Crosses represent best fit value of K z from measurement of the stress on the piston head, σ p n using Eq (6). K z and K y , represented by dots and triangles respectively, are calculated directly from the coarse grained granular packing. (E) Threshold stress σ T . Dots represent the mean value of the coarse grained continuum field σ xx at x = L, and crosses represent predicted values from Eq (5) using measured values of ν m , θ and φ 0 . stress of σ T = 10.7 was found to reproduce the observed pattern formation behaviour in the quasi-static limit, at D = 5, for a range of values of φ 0 ≤ 0.5.

Inertial number
Finally, we wish to comment on inertial effects in such a system. Towards this end we systematically vary the inertial number I from 10 −2 to 10 while maintaining D = 5 and φ 0 = 0.5 ± 0.05. As shown in Figure 8, a transition occurs at I ≈ 0.1, where the quasistatic behaviour begins to be dominated by inertial effects, and the system is fluidized. In Figure 8A, we notice a jump in the maximum solid fraction, as the particles begin to flow and rearrange due to the increased piston velocity. This is accompanied by an increase in the slope angle θ (Figure 8B), and a decrease in the accumulation rate ( Figure 8C), as the grains begin to slip relative to the walls. With increasing piston velocity, the anisotropy of the system is lost, as shown in Figure 8D, and both K y and K z tend towards a mean value of 0.56 ± 0.02 for I ≥ 1. As shown in Figure 8E, the threshold stress, σ T , also diverges above I = 0.1 away from the theoretical prediction. For values of I ≥ 0.1, we have therefore used the measured value of σ T in the best fit estimation of K z shown in 8D, rather than the value predicted from (5), as used in all other cases.

Conclusions
We have here described a large number of simulations of granular material which have been compacted in a confined geometry. For all cases, we observed that the stress distribution within the packing is well approximated by previous models, once a more rigorous definition of the threshold stress is used. This is true for a wide range of gap spacings, initial filling fractions and piston rates.
In this study we have used a rough boundary condition, where macroscopic friction at the piston and walls is equal to the inter-particle friction. However, in many systems we expect the roughness at the boundaries to be lower than that between particles. It is unclear how this difference will affect either the accumulation of material near the piston head, or the stress distribution within the packing.
Below a gap spacing of 3 particle diameters, the stress distribution is not well represented by this model. We conclude that D = 3 represents the smallest system size which may reasonably be described by the one dimensional Janssen stress model. In addition, at inertial numbers of I ≥ 0.1, we find that there is significant slip at the boundary, and the threshold stress diverges from the model prediction. In all cases, we cannot as yet predict the slope of the free surface in the transition zone, but we observe that this slope approaches the angle of repose for large systems at low piston rates. Janssen stress coefficients for this system are well represented by K z = 0.6 ± 0.1 and K y = 0.7 ± 0.1 for a wide range of system parameters. A model for the threshold stress has been presented using limit equilibrium, and this holds well for systems with D ≥ 3 and I < 0.1.
The slope angle, θ, has been measured for different system parameters to lie in the range of 2 • -18 • . A priori, we could only assume that this angle must be less than or equal to the angle of repose, which for these grains is approximately 20 • . The wide variability of θ is as yet unexplained, and is in stark contrast to the case where a top boundary does not exist, for example in fold and thrust belts [3]. We do note, however, that at large values of D the slope angle approaches the angle of repose.
With regards to the two Janssen parameters, we can clearly distinguish the values of K y and K z in Figure 6C, Figure 7C and Figure 8D, for I < 0.1. The reason for the difference between these two quantities may either be due to anisotropy in the granular packing, or due to the differing boundaries in the y and z directions. As the Janssen parameters are relatively insensitive to the gap spacing D, we conclude that this anisotropy is due to the accumulation process, which creates a preferential direction within the packing. This distinction is important when considering models which account for more complex geometries, such as in [4].

Disclosure/Conflict-of-Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Author Contributions
BM conducted the simulations, post-processing and authored the paper. BS conceived the idea for the simulations. BS, GD, JAE and KJM assisted with the interpretation and analysis of the simulations, and editing of the manuscript.