FEM modeling of the sandpile dip effect

Previously, it was shown that the pressure “dip” in the stress profile observed at the bottom of sandpiles prepared by successive avalanches can be related with the skewness effect observed inthe stress response of sheared granular layers to a localized vertical overload. Here, we present a FEM study of the sandpile dip effect using an alternative geometry (“wagon”) to build the pile, considering a simplified anisotropic elastic model, with an orthotropy angle. We explore the effect of different boundary conditions and grid coarsening on numerical stress response profiles. Our results indicate that the wagon geometry combined with partially fixed boundary condition may describe better experimental results than the conical geometry, a remark which can contribute for better comprehension of this effect.


Introduction
One of paradigmatic features in the physics of static granular materials is the sandpile stress dip effect [1][2][3][4].Following an intense debate in the literature along several years, first, to recognize this effect [5][6][7] and, later, to furnish a fair theoretical description [8][9][10].Nowadays, it is well established that this dip effect is due to the preparation history of the sandpile [11][12][13][14][15][16] which lead to a spontaneous pattern formation on the force chains along the layer.Aside the debate along the mechanical origins for the dip, a rich discussion concerning how forces are transmitted along granular materials have emerged in this context, forming a theoretic panorama for stress propagation in granular layers: hyperbolic partial differential equations (PDE) are associated with the stress propagation in ordered materials, while elliptic PDE are associated with the stress propagation in disordered assemblies [17], but transition between these behaviors can occur under given conditions of friction or force magnitude [18,19].
In this context, a tool initially designed to distinguish between the very different mathematics of hyperbolic and elliptic PDEs was proposed: the stress response function (SRF) [8,20].It consists basically to measure the pressure profile at the bottom of a layer in response to a localized overload at its top surface.Elliptic PDEs predict a single broad peak for this stress response while hyperbolic models would yield two peaks in 2D or a ring in 3D.Several experiments [11,14,[21][22][23] and simulation studies [15,16,[24][25][26] have addressed this issue and the above scenario for the stress propagation in granular materials is well established today and fairly verified [4,27].
Nevertheless, there are still several unresolved questions.For instance, Atman et al. [15] have proposed a relation of the stress dip observed in sandpiles prepared with successive avalanches with the skewness angle observed in SRF profiles of sheared granular slabs [14], a remarkable example of strain induced anisotropy.However, the comparison of experimental results and numerical data did not presented a good match, except when one of the control parameter was set to infinity [15].In that paper, the authors have suggested that the geometry used in the construction of the 3-D conical piles maybe could have some influence in the disagreement observed in the results, since they considered a 3-D conical pile produced by rotation of a triangular element around the vertical axis.This geometry led to a sharp tip, a point in the mathematical sense, where the overload was applied on the top surface of the sandpile.Here, we present an extension of this work considering a trapezoidal geometry to build the pile -Figure 1, and study the dependence of the skewness angle β with the orthotropy angle τ of the layer.We also study the effect of using different boundary conditions for the substrate in the SRF: completely fixed or partially fixed, and also different grid coarsening.
The paper is structured as following: the next section is devoted to present the basic methodology used to build the trapezoidal layer ("wagon") in CAST3M T M [28] language, an open source for Finite Element Method (FEM).We also present the analytical treatment used to obtain the SRF function with this geometry and the details of the numerical model.Next, the results obtained for the range of parameters and boundary conditions tested are shown and discussed.Finally, we draw some conclusions and perspectives to future investigations and our acknowledgements.

Methods
Here we present briefly the elastic calculations considered to build the model and obtain the stress response function profiles.Due to didactic purposes, we first will present the 2D case before to generalize the result for the 3D numerical model.In both 2D and 3D, isotropic elasticity describes a given material by means of two constitutive parameters, the Young modulus E and the Poisson ratio ν.A third constitutive parameter, the shear modulus G, can be written in function of E and ν in the isotropic case.The anisotropic case is a bit more complicated, since the number of parameter grows fast with dimension.Here, we will consider an anisotropic elastic model with an orthotropy angle, which leads to five independent parameters in 2D [17,22], as shown below.A meaningful fitting study of stress response profiles with this model was published recently [26].
The motivation to use an anisotropic model with an orthotropy angle is the possibility to reproduce a well known experimental phenomena of granular assemblies submitted to plane shear, the tilting of chain forces at an angle of 45 o with respect to the shear direction.This rearrangement of force chains lead to a reinforcement of contacts in that direction and, consequently, a softening of those contacts in the perpendicular directions.Thus, the Simple Elastic Model (SEM) in 2D [16] considers a stiff direction (direction 1) which makes an orthotropy angle τ relative to the vertical axis z (x is the horizontal direction), and a softer direction (direction 2) with a smaller Young modulus E 2 < E 1 .Along with two Poisson coefficients, ν 12 and ν 21 , and the shear modulus G, we can write the strain-stress relation in the matrix form, expressed in these principal stress directions: The symmetry properties of the compliance matrix leads to the following E 2 /E 1 = ν 21 /ν 12 , and the requirements for a proper definition of the elastic free energy are: E 1 , E 2 , G > 0 and 1 − ν 12 ν 21 > 0. The isotropic limit is obtained with Otto et al [17], have studied analytically the stress tensor components for the SEM model in the case of a semiinfinite medium (z > 0) loaded with a unitary point force at the origin x = z = 0.For a given anisotropy angle τ, they have shown that the SRF cna be described with only two parameters: For a layer of finite thickness, an additional parameter is needed to be describe the profiles [15,22].The extension of the model to consider an arbitrary orthotropy angle τ and finite thickness for the layers was presented recently [25,26].To recover the stress components at the normal and tangential coordinates, a rotation matrix should be applied, where W † is the compliance matrix.The explicit form of W can be found in [26].The generalization of this approach for 3D is immediate, considering two axes (directions 2 and 3) orthogonal to direction 1 as isotropic and soft directions, in the principal stress reference system.In the real space, besides the horizontal direction x along which the shearing is applied, and the vertical direction z, a third direction, y, perpendicular to these two, is now considered.The 3D equivalent of the compliance matrix has a similar structure to the 2D case, but with three Young moduli E i , three shear moduli G i and six Poisson ratios ν i j (i and j must be different).Again, the symmetry of the matrix gives three relations of the form ν i j /E i = ν ji /E j .In the final analysis, there are 9 independent coefficients besides the anisotropy angle τ.Stability now requires that E i , G i > 0, 1 − ν i j ν ji > 0 and 1 − i j ν i j ν ji − ν 12 ν 23 ν 31 − ν 21 ν 32 ν 13 > 0 .
To our knowledge, there is no closed analytical expression for the stress response profiles in 3D up to now, and we can use a Finite Element open source code developed in soil mechanics called CAST3M T M [28] to compute the stress components.We test three grid coarsening values to observe the sensibility of results with the grid resolution: "Model 1" ( 41 × 13 × 13 ), "Model 2" ( 41 × 15 × 15 ) and "Model 3" ( 43 × 13 × 13 )cells in the x, y and z directions respectively.We also considered two boundary conditions for the substrate: a completely blocked condition -"Condition 1", where no displacement is allowed in any direction (U z = U x = U y = 0); and a partially blocked condition -"Condition 2", where only the corners of the substrate and the vertical direction (U z = 0) are blocked, while in (x, y) directions small displacements are allowed.
Following the reasoning of Atman et al., we limit the parameter space to better understand the role of each player: the average stiffness E = (E 1 +E 2 +E 3 )/3, the shear modulus G and τ, the orthotropy angle, measured with respect to the vertical direction.Thus, we fixed the same value of G for all directions - ), with ν i j = 0.3 and E = 150MPa.In the spirit of the 2D analytical results, we use two rescaled parameters to present our data, i.e. the stiffness ratio t = E 2 /E 1 , and the shearing ratio u G = E 2 /G.This last parameter is analogous to the parameter r in 2D (see equation 2).For the anisotropic samples, we have considered t = 1.5, and keep all other parameters fixed.

Results
An essential result from shearing experiments in 2D and in 3D was that the orthotropic direction is tilted with an angle τ = 45 o with the shear direction, leading to a skewness angle of ∼ 22o in 2D and around 8 o in 3D.By means of a fitting procedure Atman et al. have obtained t = 1.5 and u G = 2.14 for the best fit of SRF in shearing experiments with sandbox.When applying these values in a conical pile numerical finite element model, build a pile from the rotation or a triangular around the z axis, the authors did not found a good fit (skewness angle η = 3 o ).The best fit for the SRF from sandpile experiments seems to correspond to the limit u G → ∞ [15], a result that still is obscure.Here, we test a different geometry, that we called "wagon", in the hope to contribute to the understanding of the reasons for this mismatch.
In Figure 1 we show the "wagon" geometry that we have considered to build the pile, and in Figure 2 the corresponding SRF profiles.These results were obtained considering the SEM model with t = 1.5 and u G = 2.08, but τ = 0 o to test only the influence of the grid coarsening and boundary conditions.It is clear that the coarse-graining has a small influence on the SRF profiles, but the boundary conditions plays a less important role.The partially blocked condition furnished results which seem to better reproduce the profiles obtained in experiments or simulations.The comparison of the results for the different grid coarsening is shown in Figure 3.
The dependence of the SRF profiles with the orthotropic angle is shown in Figure 4 for each grid coarsening considered and different boundary conditions.Figure 5 show the summary of the results obtained, displaying the dependence of the skewness angle with the control parameters considered here.The main result is that the skewness angle can be considered enhanced with a proper combination of boundary conditions and orthotropy angle.The grid coarsening plays a role on the stress profiles, it is related mainly to the steepness of the response peak.This result indicates that the mathematical feature of the conical pile, with a perfect sharp tip maybe was the cause for the mismatch of the results.However, we believe that we have to better improve the aspect ratio of the wagon geometry in order to reproduce better the experimental SRF profiles of avalanched sandpiles.Moreover, we hope to test other values for t and u G in order to understand better what is the reason for the mismatch in the skewness angle, which still is verified here.

Conclusion
We present an extension of previous study of the FEM modeling of stress dip effect in sandpiles prepared by successive avalanches.We test different boundary conditions and grid coarsening combinations to study their influence on the stress response profiles.Our results indicate that a proper combination of boundary condition, grid coarsening and elastic parameters could furnish results which could fit better the experimental and simulation results.We hope to explore the parameters space and also test other grid coarsening combinations to understand better the reason for the mismatch in the skewness angle which persists with the wagon geometry.

Figure 1 .
Figure 1.The "wagon" geometry.We show two different side views of the prismatic pile considered to measure the SRF profiles.Above: x − z plane.The SRF profiles are measured along the x axis at the bottom of the pile.Below: y − z plane.Note the translational symmetry along the y direction and the increasing of resolution (finer grid coarsening) at the center of the pile, where the SRF is measured.

Figure 2 .
Figure 2. Comparison of stress response profiles for different boundary conditions.Panels a), b) and c) show results for Model 1, Model 2 and Model 3, respectively, with τ = 0 o .Panel d) show the result for Model 3 and τ = 15 o .

Figure 3 .
Figure 3.Comparison of the SRF considering different boundary conditions.Left panel show results for the three models considered and fixed boundaries.Right panel show the results considering partially fixed boundaries.

Figure 5 .
Figure 5. Skewness angle variation in function of the orthotropy angle for the different models considered, as indicated in the insets.