Influence of the feeding mechanism on deposits of square particles

In a previous paper [Hidalgo et al., Phys. Rev. Lett. 103, 118001 (2009)] it was shown that square particles deposited in a silo tend to align with a diagonal parallel to the gravity, giving rise to a deposit with very particular properties. Here we explore, both experimentally and numerically, the effect on these properties of the filling mechanism. In particular, we modify the volume fraction of the initial configuration from which the grains are deposited. Starting from a very dilute case, increasing the volume fraction results in an enhancement of the disorder in the final deposit characterized by a decrease of the final packing fraction and a reduction of the number of particles oriented with their diagonal in the direction of gravity. However, for very high initial volume fractions, the final packing fraction increases again. This result implies that two deposits with the same final packing fraction can be obtained from very different initial conditions. The structural properties of such deposits are analyzed, revealing that, although the final volume fraction is the same, their micromechanical properties notably differ.


I. INTRODUCTION
In a static granular sample, forces are transferred from particle to particle through their contacts, giving rise to a network of force chains [1][2][3][4] which is behind the microand macromechanical properties of the packing.In deposits formed under the action of gravity, this particular stress transmission leads to characteristic phenomena such as the pressure dip below the apex of a pile [5,6], or the Janssen effect [7], i.e., pressure saturation in silos.Moreover, it is well known [8][9][10] that these phenomena strongly depend on the packing formation history.In particular, the pouring process has been shown to have a remarkable influence on the global mechanical properties of the packing [8,11].
In the past few years, interest is growing in the study of the effect that particle shape has on the properties of both piles [12] and silos [13][14][15] as well as in other configurations [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31].Among these works, some have focused on the role of faceted particles which have revealed conceptually new scenarios and special features.Very recently, the deposit of flat faceted particles in a silo was studied, in the limiting case of very low feeding rates [13,14].There, it was found that when elongated particles are deposited under the action of the gravity, they tend to orient horizontally.This orientation favors vertical transmission of the stress, which increases almost linearly with the height of the layer.Hence, pressure saturation at the bottom of the silo is not achieved, at least for the height of the columns investigated.Nevertheless, when the aspect ratio of the particles is decreased, their orientation deviates from the horizontal, giving rise to stress transmission towards the walls, which implies a faster pressure saturation.For the limit case of squares, a well marked preferred orientation with the diagonal parallel to the gravity is observed.This surprising configuration promotes strong transmission of stress towards the sides and a quick saturation of the pressure at the bottom of the silo.
In this paper we explore the effect that the pouring mechanism has on the structural properties of deposits of square particles in a rectangular silo.The different pouring consists, basically, of feeding the grains with a fixed vertical velocity but from initial configurations with different volume fractions.The volume fraction has been chosen as a control parameter for simplicity's sake.In addition, this protocol has the advantage of allowing a straight comparison with existing works in other similar scenarios such as tapping of grains [32][33][34][35] or deposits of particles at controlled velocities [36,37].
The paper is structured as follows.First we will explain the experimental procedure in Sec.II and the numerical simulations in Sec.III.Then, in Sec.IV we describe the morphology of the packings obtained from different pouring mechanisms by analyzing the final volume fraction and the orientation of squares within the deposits.Afterwards we focus on the effect that the pouring mechanism has on the micromechanics of the deposits.We conclude in Sec.V, drawing the main conclusions and implications of our results.

II. EXPERIMENTAL PROCEDURE
The experimental setup consists of a two-dimensional silo (Fig. 1) which is filled with square particles starting from initial configurations with different volume fractions.The square particles are DIN 557 nuts whose lateral sides are 6.9 mm and the width is 3.16 mm.The whole setup (Fig. 1) has two compartments: the proper silo (D) at the bottom, and a reservoir (B) at the top from which the particles are poured into the proper silo.Each compartment was built with two glass plates of 6 mm thickness separated by two stainless steel strips 3.3 mm thick.Note that the strips are 0.14 mm thicker than the squares.As a result, the particles stay confined in a monolayer between the plates.The horizontal separation between the two strips is W = 258 mm, which is, indeed, the width of the two-dimensional silo.
The proper silo (D) was built with two glass plates of 900 mm height and 500 mm width.The bottom is flat and formed by two metal pieces whose edges touch each other.This configuration allows the discharge of the particles by separating the pieces.The number of square particles necessary to fill the silo is around 5 × 10 3 .The top reservoir (B) was built with two flat glasses of 750 mm height and 500 mm width.At the bottom of the reservoir (C), a removable metal piece of 258 mm length and 1 mm width is placed to prevent the particles from passing into the proper silo.
In order to change the pouring rate, a number of obstacles were placed in a regular arrangement at the position (C).Hence, four obstacle configurations have been considered.Moreover, an extra experiment was performed by pouring the grains homogeneously along the hopper at the top.In this way, we obtained a very low initial volume fraction and, consequently, a very low feeding rate, which is comparable with the one reported in Ref. [15].
The experimental procedure is as follows.First the metal piece is placed in the lower part of the reservoir, which is then filled through the hopper at the top (A).When all the particles are deposited, the metal piece is removed and the grains fall freely under the action of gravity to the proper silo.Once the deposit is generated, an image is recorded with a standard 10.2 megapixels camera.The region registered covers the whole width of the silo and a height that goes from 290 to 410 mm from the bottom of the silo (Fig. 1).The images recorded are analyzed in order to obtain the center of the particles and their orientation.For each initial configuration, 100 deposits were performed to obtain good statistics.The initial volume fractions are measured using a high speed camera at 3000 frames/s.The recorded area covers a height from 330 to 370 mm and the whole width of the silo.For each experimental condition, five depositions are recorded for about 1 s.From the images, we have estimated the initial volume fraction (φ 0 ) of the deposition process and the mean particle velocity in the vertical direction (v y0 ).In these experimental conditions, the vertical velocity of the particles at the region of measurement is practically constant, v y0 = 350 ± 10 cm/s, for all the obstacle configurations.In that sense, these experiments resemble the ones reported in Ref. [37], although in our case gravity is at play.

III. NUMERICAL SIMULATIONS
We have performed discrete element modeling of a twodimensional (2D) granular system composed of spheropolygons with four vertices [38][39][40].Specifically we used undeformable squares of side d and area A = d × d.The system is confined within a rectangular box of width W = 36 × d and with lateral and bottom boundaries built of fixed particles.To mimic the real experimental scenario, the particles are continuously added at a given position y = h 0 with orientation θ 0 and volume fraction φ 0 .The initial velocity on the y axis is set assuming that the particles have been thrown without an initial velocity from a height of 144 × d with respect to the bottom of the silo; and the velocity is considered to be zero on the x direction.In each case, the values of h 0 , φ 0 were carefully chosen to be as close as possible to the experimental conditions.Once the particles have been added in the system, they settle under the effect of gravity and relax until their mean kinetic energy is several orders of magnitude smaller than its potential energy.
To better describe the deposit formation process and its micromechanical properties, we executed simulations starting from two different initial conditions, concerning the angular degree of freedom of the squares.In protocol I the initial angular orientation of the squares is drawn from a uniform distribution, p(θ 0 ) = 1/2π .In protocol II the angular orientation was fixed to θ 0 = π/4, i.e., the diagonal of the squares is oriented with gravity.In both cases the initial rotational energy was set to zero.In every case, we have simulated 4 × 10 4 squares and, to get good statistics, the results presented for each pouring mechanism constitute averages over at least eight different initial configurations.
In the simulation, each particle i (i = 1 . . .N) has three degrees of freedom, two for the translational motion and one for the rotational one.The particles' motion is governed by Newton's equations of motion.For calculating the particles' interaction F ij we use a very efficient algorithm proposed recently by Alonso-Marroquín et al. [38,39], allowing the simulation of a large number of particles.This numerical method is based on the concept of spheropolygons, i.e., a polygon i is defined by the set of vertices V i and edges E i .Moreover, the force calculation is then executed only over neighboring particles and a link cell algorithm is used to allow rapid examination of each neighbor list [41].
Following this method, the interaction is equivalent to the interpenetration between the two neighboring disks (see the sketch in Fig. 1).Thus, the force F ij exerted on particle i by the particle j is defined by where F (V ,E) represents the force between the vertex V and the edge E of each contacting pair ij .Hence, the local interaction between two contacting particles is only governed by the overlap distance δ between the vertex and edges.In Eq. ( 1), each term F (V ,E) can be decomposed as where F N is the component in the normal direction n to the contact plane.Complementary, F T is the component acting on the tangential direction t.To define the normal interaction F N , we use a linear elastic force, proportional to the overlap distance δ.Moreover, to introduce dissipation, a velocity dependent viscous damping is assumed.Hence, the total normal force reads as where k N is the spring constant in the normal direction, m r = m i m j /(m i + m j ) = m/2 stands for the pair's reduced mass, γ N is the damping coefficient in the normal direction, and v N rel is the normal relative velocity between i and j .The tangential force F T also contains an elastic term and a tangential frictional term accounting for static friction between the grains.Taking into account Coulomb's friction constrain, which reads as rel is the tangential component of the relative contact velocity of the overlapping pair.ξ represents the elastic elongation of an imaginary spring with spring constant k T at the contact [42], which increases as dξ (t)/dt = v T rel as long as there is an overlap between the interacting particles [42,43].μ is the friction coefficient of the particles.
In order to model hard particles, the maximum overlap must always be much smaller than the particle size; this is ensured by introducing values for normal and tangential elastic constants, k t k n = 0.1 and k n = 10 4 N/m.The ratio between the normal and tangential damping coefficients is taken as ν n ν t = 3, while gravity is set to g = 10 m/s 2 .For these parameters, the time step was set in t = 10 −6 s.We have carried out numerical simulations in which we changed the damping coefficients and the static friction coefficient.Thus, we have ensured that the kinetic energy loss and the dynamics of sediment formation are analogous to those seen experimentally.We have converged to ν t = 10 2 s −1 and μ = 0.4 as best fit parameters.In all the simulations reported here, we have kept constant the previous set of parameters and only the volume fraction of the initial configuration has been modified.

IV. MORPHOLOGY AND MICROMECHANICS
In Fig. 2 we present photographs of the experimental deposits of square particles obtained when pouring the grains from different initial volume fractions φ 0 .Clearly, there is an influence of the initial volume fraction on the final deposit morphology.For very low values of φ 0 the particles arrange oriented with their diagonal in the direction of gravity in very well ordered configurations.As φ 0 is increased, this order is lost and the preferred orientation is less well defined.In what follows we present a detailed quantitative analysis of the morphology of these deposits.The first question that should be tackled is if the deposits have homogeneous properties along the whole width of the silo.In a previous work we observed that, apart from a region adjacent to the lateral walls where the particles tend to be aligned with them, the spatial dependence is negligible [15].Here, we check if this behavior holds for the different pouring mechanisms.In order to do so, we measure the probability that the squares placed at different distances from the lateral walls of the silo are oriented vertically.In practice, from the images, we obtain the particles whose center is within vertical slices of 1 cm width.From all these particles, we measure the probability that their orientation (the orientation of one face) is at θ = π/2 ± 0.035 with respect to the horizontal.Note that this orientation is equivalent to θ = 0 ± 0.035.In Fig. 3 we present these probabilities for slices at different distances of the wall (x) for three different initial volume fractions.Clearly, the higher the initial volume fraction, the stronger is the vertical alignment of the squares.In any case, independently of the pouring conditions, the effect of the lateral wall becomes negligible at 4-5 cm apart from it.Hence, from now on, in order to avoid wall effects, the region where we will perform our analysis will comprise the central part of the silo excluding the 4 cm adjacent to each of the walls.
In Fig. 4, we show the results of the volume fractions obtained in the final deposits φ f for different values of the initial volume fraction φ 0 .From the experimental data, it becomes obvious that, up to φ 0 ≈ 0.25, the higher the initial volume fraction, the lower is the final one.This result is in excellent agreement with experiments performed in a belt where particles were deposited at constant speeds [36,37].Interestingly, when the squares are thrown at an initial volume fraction of φ 0 ≈ 0.5, the final volume fraction increases again.
An important consequence of this nonmonotonic behavior of φ f with respect to φ 0 is that two deposits with the same final packing fraction can be obtained by using different initial conditions.This nonmonotonic behavior has also been reported in Ref. [37] and is strongly reminiscent of the packing fraction behavior in tapping experiments of granular layers.Indeed, it has been shown, both experimentally [33,34] and numerically [32,35], a nonmonotonic behavior of the volume fraction with the intensity of the tap.Note that in Ref. [35] the numerical simulations were also performed with .For the simulation we illustrate, as indicated by the legend, the results obtained from a system of particles deposited with an initial random orientation distribution (protocol I) and from an ordered system of particles with θ 0 = π/4 (protocol II).faceted particles.Finally, let us stress that the experimental results reported in Fig. 4 are nicely reproduced by the two pouring protocols implemented numerically.There is a good quantitative agreement for small values of φ 0 , while for high values of φ 0 the differences appear.
A simple explanation for the nonmonotonic behavior of the final packing fraction versus the initial one can be borrowed from experiments of tapping, considering that high intensities of tapping correspond to low values of φ 0 and low intensities of tapping are related to high values of φ 0 .For very low initial volume fractions, particles have enough time to explore several possible configurations.Hence, they minimize their potential energy and settle in a well ordered arrangement.If the initial volume fraction is increased, however, the time that the particles have to explore the most stable configuration is reduced, giving rise to less ordered structures [44].Nevertheless, when the initial volume fractions are exceedingly high, particles seem to have a certain "memory" of their initial configuration and the final volume fractions start to increase with the initial ones.Interestingly, the final deposits attained with the two numerical protocols (where the particles are thrown with different initial orientations) are very similar for small values of φ 0 .This indicates a weak (or null) dependence of the final packing on the initial configuration.However, for very high values of φ 0 , both protocols display differences and considerably higher values of φ f are obtained with protocol II than with protocol I. Clearly, this may be due to the initial order with which the particles are thrown in protocol II (recall that the initial orientation of all the particles was fixed at θ 0 = π/4) and could be related with the existence of high correlations between the particles during the deposition process.
The next variable that is experimentally accessible is the orientation of the squares in the final deposits.This variable has been shown to be strongly correlated with the force transmission within the silo [13,14].In Fig. 5 we present the angular distributions of the squares obtained, both experimentally and numerically with protocol I, using three different initial volume fractions.The excellent agreement of the results validates the numerical scheme which will be used to analyze the contact forces among the particles, a parameter which is difficult to obtain experimentally.The main result revealed by Fig. 5 is that, independently of the initial volume fraction, the squares tend to align at θ = π/4, i.e., with their diagonal parallel to gravity.In addition, the lower the initial volume fraction is, the higher is the probability that the squares orient at θ = π/4.This result can be understood if we assume that θ = π/4 is the most stable configuration for the square particles.Hence, when φ 0 is low, the particles have enough time to explore different configurations and settle at θ = π/4 preferably.On the contrary, when the particles fall from denser initial states, they cool down without time to explore different configurations and the probability of finding a particle oriented at θ = π/4 is reduced.
In order to further explore this result, in Fig. 6 we report the values of the maximum value of the probability density [f (θ = π/4)] as a function of the initial volume fraction.Focusing on the experimental results, it becomes clear that the amount of squares oriented with their diagonal parallel to the gravity is very high when the initial volume fraction is very low.Then, up to φ 0 ≈ 0.25, the number of squares oriented at θ = π/4 is reduced.Interestingly, for higher initial volume fractions, f (θ = π/4) seems to increase again.In principle, one can expect that this behavior is a consequence of the very particular configuration from which the experimental squares were thrown.Indeed, as the particles fall from a previously formed deposit at the top reservoir, it could be expected that they fall with their diagonal parallel to gravity.Hence, for very dense initial configurations, the squares may not have space to modify their orientation and their tendency to orient at θ = π/4 could increase again.However, this assumption was invalidated by the numerical results as the increase in the values of f (θ = π/4) is observed with both protocols, when the particles are thrown with random orientation and with their diagonal parallel to gravity.
This surprising result encouraged us to study in detail the orientations of the squares in the deposits obtained for very high initial volume fractions.In Fig. 7 the distributions of orientations are reported for three different cases: experiments [Fig.7(a)], numerical protocol I [Fig.7(b)], and numerical protocol II [Fig.7(c)].The behavior obtained for the case of very low initial volume fractions is very similar for all the cases.The only apparent difference is that protocol II-where the particles are thrown with θ 0 = π/4-seems to enhance the sharpness of the peak at θ = π/4.For high initial volume fractions the differences among protocols are more evident.Surprisingly, for high values of φ 0 , numerical protocol II gives rise to a local minimum at θ = π/4 in the orientation distribution.This unexpected result is also behind the fact that, in Fig. 6, the values of f (θ = π/4) for φ 0 = 0.30 are higher for protocol I than for protocol II.
It is a fact that the two numerical protocols give rise to similar angular distributions for low φ 0 .High values of φ 0 , however, reveal important differences between the two protocols.Those outcomes suggest that, in the packing formation, the initial conditions are specially relevant for very high feeding rates.Furthermore, we observe that, by using the same protocol, different values of initial volume fractions can lead to the same value of final volume fraction.For example, the final volume fractions obtained for φ 0 = 0.05 and φ 0 = 0.25 (Fig. 4) are very similar.However, these deposits are clearly different as they display distinct orientation distributions (see Fig. 6).
Thereby, to highlight the spatial disorder induced by the feeding mechanism, we have examined the radial distribution function, which reads as and where φ f is the average packing fraction in the deposit.N (r + δr) accounts for the number of particles with their center of mass at a distance r within a differential of area δS = 2πrδr.In Fig. 8 the radial distribution function g(r) is illustrated for experimental and numerical depositions of squares.In all cases the first result that becomes apparent is the existence of peaks at r/d = 1,2,3 . . ., which correspond to structures where squares are perfectly aligned with their faces parallel to each other.As the initial volume fraction gets larger, the final deposits present a more disordered structure and consequently the intensity of the peaks decreases and their location deviates from the exact values.Nevertheless, for very large initial volume fractions, the configuration of perfect aligned squares is again favored.Complementary, we have carefully explored the impact that the feeding rate has on the packing properties at contact scale.In Fig. 9 we present, for both numerical protocols, the number of face-to-face contacts as analyzed in previous works [26,45].For comparison the data values are rescaled with the mean number of contacts N c in the region of study.For low values of φ 0 the higher the initial volume fraction, the lower is the number of face-to-face contacts.Moreover, in this region the results obtained for both protocols are very similar.After certain values of φ 0 , to increase the initial volume fraction results in a further increase of the relative number of aligned contacts.As we pointed out earlier for such large initial volume fractions, the mean distance between particles in the yield front is very low.As a result, the cooling process is mainly governed by the local correlations between neighboring particles, which gives rise to the formation of well defined faceto-face interactions.This behavior is enhanced in protocol II, where the initial particle orientation (θ 0 = π/4) gives rise to an increase in the number of face-to-face alignments within the final deposit.Interestingly, those findings correlate with the behavior of the global volume fraction of the packing, which denotes the remarkable influence of the local correlations of the global properties of the deposits.
Microstructure and morphology changes induced both by particle geometry and the feeding mechanism might determine the stress transmission within the column.Nowadays, it is well established that the stress acting on a single particle i can be defined in terms of the contacting forces F c i and its corresponding branch vectors l c i [46,47].Hence, the local stress tensor can be calculated as where l c i is the branch vector related to the contact c and V p accounts for the particles volume.In Eq. ( 3), the sum runs over all contacts of particle i.
We numerically describe the stress transmission within deposits of squares, which have been formed following the two protocols explained above.Employing Eq. ( 3), the stress tensor on each particle within the representative volume element (RVE) has been computed.Therefore, the stress eigenvalues and their corresponding eigenvectors have been examined.To clarify the preferential directions on which the stress is transmitted within the silo, we have studied the polar distribution of the stress eigendirection n11 .We have chosen the eigendirection related to the largest eigenvalue σ 11 of the stress tensor and its polar distributions obtained for different pouring mechanisms are illustrated in Fig. 10 (numerical protocol I) and Fig. 11 (numerical protocol II).Since the studied packing has been built under the action of gravity, the polar distributions of the stress eigendirections show an intrinsic anisotropy.That is, they describe the preferential directions for the stress transmission with respect to the external field, in our case the up-down direction.
The statistics corresponding to very low initial volume fractions [Figs.10(a) and 11(a)] reveal that forces are mainly transmitted along the π/4 and 3π/4 directions.It is noticeable that the outcomes are very similar for both deposition protocols.As explained above, when the initial volume fraction is low, the squares cool down, exploring many energy states until reaching the most stable configuration.This process seems to be independent of the initial angular orientation of the squares.Increasing the initial volume fraction, however, leads to final packings which display a preference for the stress transmission parallel to gravity.This behavior seems to be caused by the formation of less dense packings and the enhancement of the  disorder.Furthermore, this correlates with the decrease of the number of face-to-face interactions (see Fig. 9).Moreover, it is again interesting to compare two deposits that have a similar final packing fraction but have been built from different initial volume fractions.It is noticeable that the outcomes presented in Figs.10(c) and 11(c) indicate a stronger vertical stress transmission than those illustrated in Figs.10(b) and 11(b).These differences are notably enhanced for protocol II, when particles initially have their diagonal parallel to gravity.Finally, the stress polar distributions obtained for very large initial volume fractions [Figs.10(d) and 11(d)] suggest again a certain prevalence of the φ = π/4 direction.As we pointed out above, this is related with the formation of better define face-to-face contacts, which is mainly induced by the initial conditions and the strong local correlations.

V. CONCLUSIONS
In this paper we have experimentally and numerically investigated the role of the pouring mechanism on the properties of deposits of square particles in a silo.The first noticeable result is that squares tend to align with their diagonal parallel to gravity, independently of the pouring mechanism utilized.This behavior confirms that such an orientation is the one preferred by the squares, even for very high feeding rates.The robustness of the orientation θ = π/4 over the a priori expected one (faces aligned with the walls) correlates with the fact that the diagonal parallel to the gravity is the orientation that maximizes the contact area between squares [13].This surprising alignment promotes arching, which enhances the stability against perturbations of this configuration in relation with the horizontal one.Certainly, there is a region near the bottom of the silo where the horizontal wall induces the alignment of squares with it.Nevertheless, the wall effect is lost as we move upwards within the silo.In fact, at the height where we perform our analysis, we make sure that the pattern is fully developed.
For low values of initial volume fraction, it is observed that the higher is the initial volume fraction, the lower is the final volume fraction of the deposits.This could be explained by the reduction of the mean distance between particles in the yield front.This would imply that the time that particles take to cool down decreases as the initial volume fraction increases.Hence, the disorder in the final packing increases with the initial volume fraction [44].Surprisingly, very dense initial configurations reverse this tendency, leading to higher values of the final packing fraction.This phenomena can be understood if we consider the existence of strong spatial and temporal correlations related with the initial particle arrangement.This statement is supported by the fact that the results of the two numerical protocols display stronger differences for high initial volume fractions than for low ones.An important consequence of this nonmonotonic behavior of the final packing fraction with respect to the initial one is that deposits with the same final packing fraction can be achieved starting from different initial conditions.Throughout the paper it has been shown that, despite the fact that two deposits obtained with different initial volume fractions may have the same final packing fraction, their mechanical properties are very different.In particular, in deposits obtained for high initial volume fractions the stress transmission in the vertical direction is favored when compared with the deposits obtained for low initial volume fractions.

FIG. 1 .
FIG. 1. Scheme of the experimental setup.A, hopper; B, top reservoir; C, slot where a metal piece is introduced to keep the particles accumulated at the reservoir; D, proper silo.The shadowed region in the middle of the proper silo is the region where pictures are taken and measurements performed.The black region corresponds to the closeup of the deposit of squares presented in the figure at the right.At the top right corner a sketch of the model is presented.

FIG. 2 .
FIG.2.Pictures of the deposits of square particles obtained when pouring the grains at different initial volume fractions (φ 0 ), as indicated by the labels.

FIG. 3 .
FIG. 3. (Color online)Probability that the particles align vertically (with the direction of the lateral wall of the silo) at different distances to the wall (x) for different values of the initial volume fraction, as indicated in the legend.This probability is calculated as the number of particles aligned with θ = π/2 ± 0.035 (or what is the same θ = 0 ± 0.035) with respect to the total number of particles found at the different positions of the silo.

FIG. 4 .
FIG. 4. (Color online)Volume fraction of deposits of squares (φ f ) obtained from different initial volume fractions (φ 0 ).For the simulation we illustrate, as indicated by the legend, the results obtained from a system of particles deposited with an initial random orientation distribution (protocol I) and from an ordered system of particles with θ 0 = π/4 (protocol II).

25 FIG. 5 .FIG. 6 .
FIG. 5. (Color online) Orientation distributions of squares in the final deposits obtained in experiments and simulations following protocol I, for different initial volume fractions, as indicated in the legend.

FIG. 7 .
FIG. 7. (Color online) Orientation distributions of squares in the final deposits obtained in experiments (a) and simulations following protocol I (b) and protocol II (c) for different initial volume fractions, as indicated in the legends.

GFIG. 9 .
FIG. 8. (Color online) Spatial correlation function obtained for packing of squares in experiments (a) and simulations following protocol I (b) and protocol II (c) for different initial volume fractions, as indicated in the legends.

FIG. 10
FIG. 10. (Color online) Polar distribution of the principal direction (larger eigenvalue of the stress tensor) for deposits of squares obtained with protocol I: Squares are thrown with completely random orientations.Results are presented for different initial volume fractions, as indicated in the graphs.
FIG. 11. (Color online) Polar distribution of the principal direction (larger eigenvalue of the stress tensor) for deposits of squares obtained with protocol II: Squares are thrown with their diagonal aligned with the gravity θ 0 = π/4.Results are presented for different volume fractions, as indicated in the graphs.