Jonathan Barés

Yuchen Zhao

Mathieu Renouf

Karola Dierichs

Robert Behringer

Structure of hexapod 3D packings: understanding the global stability from the local organization

  L'archive ouverte pluridisciplinaire

Introduction

The use of aggregates made from non-convex particles or particles with anisotropic contact laws is an emerging area for both the research in physics [1][2][3][4][5][6][7] and the design of functionalized materials by tuning local properties of the grains [8][9][10][11]. This new field in granular matter science has already shown to be extremely promising for future lightweight and reversible architecture [12]. It permits to design easy-to-shape and remove strong structures. Indeed the shape of the grains [2,13,14] and the contact anisotropy [8] have been recently used to help crossing the jamming transition and to increase the stability/rigidity of granular systems. Different particle shapes, such as long and thin rods [15], staples [4], Z-shaped [6] and starshaped [2,5] particles have been shown to bring additional rigidity to 3D packings, compared with spherical shapes [16].

Two of the most striking characteristics of non-convex designed granulates, in terms of the granular system and possible design applications are their ability to form vertical structures with a 90 • angle of repose, and to sustain small tilting or loading perturbations [2]. Even if this tremendous stability illustrated in the fig. 1-A has already been evidenced in several experiments [2,5,17,18] and if packing of such particles has been widely studied [19][20][21], few is know about the fundamental reasons of why non convex particle piles are much more stable than bead columns [16,22,23].

In order to understand the origin of this rigidity, in this paper, we investigate the local organization of columns � e-mail: jb@jonathan-bares.eu made of star-shape or hexapod particles (see fig. 1-A) using X-ray CT-scan. We have performed experiments by randomly pouring identical star-shape particles (2.5 to 10cm wide) into hollow cylinders, and removing the cylinders to get stable star columns, for different preparation conditions. Then, by mean of a X-ray CT-scan and a home made post-processing code, we have reconstructed the full structure of the packing: (i) position, (ii) orientation and (iii) contact points of each particles with a very good accuracy (0.1mm). We have analysed the coordination number and the distributions of the contact point positions.

In a first part we present the experimental set-up to prepare and scan hexapod piles. In a second part the method to measure the position, orientation and contact points for each particle is explained. Then, analysis tools are detailed before presenting the results.

Experimental method

The cylindrical columns have been prepared in two different ways and we used three kinds of hexapods particles as shown in fig. 1-A. These last consisted of six orthogonal beams with square cross section, that tapered from a thickness of 2mm at the center of the particle to 1mm at the tips. The end-to-end size, L, varied between 2.5cm, 5cm and 10cm. They are made of laser-sintered white nylon PA2200, which has a friction coefficient 1.0 ± 0.3. In this paper we will only talk about 2.5cm, 5cm star packings because for 10cm stars the number of particles in the structure is not large enough to give statistically meaningful results. The protocol to make particle piles consisted in filling thin-walled PVC tube with particles by dropping them one-by-one at a steady rate from an overhead hopper as shown in fig. 1-B. The inner diameter of the PVC tube was 15cm. In some cases, during the filling process, the system could be vibrated by an eccentric that was attached to the outside of the tube as presented in the fig. 1-B. This has been shown to increase the pile stability [2]. We tuned the speed of the motor to create a ∼ 1700m/s 2 cyclic acceleration of the 16g eccentric. As shown in fig. 1-C, the tube was then carefully and slowly removed by lifting it vertically. Friction between the particles and the cylinder was low enough not to significantly perturb particles inside the cylinder. This protocol is similar to what has been used by [2].

The hexapod column was then set inside a X-ray mi-croCT scanner Nikon XT-H225 and scan with a X-ray source of 185kV and 235µA. Raw data were then postprocessed using Nikon's Feldkamp [24] cone based CT algorithm to get 2000 density slices of size 2000 × 2000, with a spacial resolution of 0.1mm and a 16bits digital resolution. A slice in the middle of the column is shown in fig. 2-A.

Data post-processing

Collecting the density slices we then formed a full 3D density matrix M 0 (see cartoon for two particles in fig. 2-B) which we post-processed with a homemade Python algorithm to extract the position x, y, z and the orientation θ 1 , θ 2 , θ 3 of each particle. Because of the large dimension of the density matrix we first resize it to M a 400 × 400 × 400 binarized matrix choosing a 60% density threshold (100% is the highest density, 0% the lowest). A virtual smaller binary matrix m(θ 1 , θ 2 , θ 3 ) is then formed. As presented in fig. 2-C this matrix is null everywhere except for the pattern of an hexapod centered in the middle of the matrix and oriented along three axis: θ 1 , θ 2 , θ 3 . Taking advantages of the particle symmetry, for all the possible families of direction with θ i going from 0 • to 90 • with 2 • steps (91125 families), we convoluted the binarized density matrix M with m(θ 1 , θ 2 , θ 3 ) to get a convolution matrix C(θ 1 , θ 2 , θ 3 ) of the same size as M. This computation took 80h parallelized on a 60 nodes cluster. Then, a maxima matrix D of the same size as M is formed. For each coordinate i, j,k it took the maximum convolution value C i, j,k (θ 1 , θ 2 , θ 3 ) over the 91125 possible orientations and keep the record of the corresponding (θ 1 , θ 2 , θ 3 ) family. As presented in the cartoon fig. 2-D, this matrix shows high convolution values 3D islands at the position of the hexapod center. Using this property, a 3D local maxima research algorithm applied on the matrix D gave the position of each hexapod x i , y i , z i and the record of their orientation gave θ 1i , θ 2i , θ 3i .

However, this measurement is not accurate enough since it has been done on a density matrix whose size has been divided by 5 (0.5mm accuracy on the particle center position) and orientation is only 2 • accurate. To improve the measurement, for each detected particle center a submatrix of the full density matrix M 0 is extracted and binarized. A larger virtual binary matrix is made and a Powell optimization algorithm was implemented to improve the match between both matrices and give more accurate position (0.1mm accuracy) and orientation (0.5 • accuracy) values. This part of the algorithm is reminiscent of the method used by [25]. Contact between the hexapods arms were finally detected from the particle position and orientation, and the density matrix M. As presented in fig. 2-E, for each arm of each particle, the distance between the segment formed by the arm and the the segment formed by each arm of the neighbour particles was computed. If this distance were lower than 4.5mm, a 10×10×10 submatrix (corresponding to 1mm 3 ) centred around the hypothetical point of contact was extracted from the density matrix M. This submatrix was binarized as previously and the number of connected domains computed. If 2 distinct domains were observed there was no contact. If there was only one large domain the two arms were in contact. This method gives contacts with a 0.1mm threshold accuracy. Further investigations are currently carried out to shows that this threshold is small enough not to detect too much fake contacts. A reconstruction of the 3D packing as well as the position of the contact is presented in fig. 3.

Results

A coordination number (average number of contact per grain) of 6.58 have been measured for non vibrated 2.5cm particles. This value is larger than the theoretical coordination 6 at the limit of stability, so this explain the stability of such a system. When particles are vibrated, as emphasized by [2], we expect to have a more stable structure. Indeed, we measure a coordination number of 7.13 for such a configuration. For particle with longer arms (5cm) we expect something even more stable [2] which corresponds with the coordination number we measured: 9.18. Hence, this explain the increase is stability for vibrated system and large particle observed by [2] but also quantify this stability by giving the hyperstaticity of the different configurations.

In fig. 4-A, we present the probability density function (PDF) of the contact position along a particle arm. For each system configuration, 2.5cm, 5cm particles, vibrated and non-vibrated, for each arm of each particle, the distance of the contacts to the particle center has been measured. Then the PDF of this distance rescaled with the arm length L has been computed. We see that in the case of large particles, contacts are homogeneously spread along the whole arm, with an over-representation close to the particle center. On the contrary, for shorter particles the contacts are more concentrated at the arm tip. This explain the more important stability of large particle piles since hexapods are more intricated so more difficult to destabilize.

Fig. 4-B present the PDF of the local packing fractions. For each system configuration, from the position coordinate of the particles, we have computed the 3D Voronoi tesselation. The volume of each tile is measured and rescaled with the an effective hexapod volume L 3 . The PDF of this quantity is then plotted. For each case it vanishes before 10%, then presents a sharp peak and vanishes more progressively for larger volumes. We remarks that the vibration does not seems to have any significant effect on the statistical repartition of the local particle densities.However, the PDF peaks correspond to a larger scaled volume for 2.5cm particles than for 5cm ones. Indeed, for larger hexapods, the harm aspect ratio is thinner so particle can entangle deeper. It is remarkable to notice that this property corresponds with a more stable system configuration.

Fig. 4-A shows the evolution of the average local density measured as the volume of the Voronoi tiles scaled with an effective hexapod volume L 3 . This quantity has been averaged vertically and circularly to give the radial evolution of the local particle density for each packing configuration. We see that the pile is slightly denser in the inside and for larger particles but constant for 5cm ones. Averaging along the vertical dimension, we get the vertical evolution of the density. Surprisingly, within the error-bar, there is no variation of the density as a function of the pile height.

Conclusions

We presented an experimental method to investigate the local structure of aggregates made from non-convex particles. This method is based on X-ray scanning and on an efficient and 0.1mm accurate homemade particle detection code using parallelized matrix convolutions. This permits to measure the position and orientation of different sized hexapods forming large piles in vibrated and nonvibrated configurations. This local scale information is in good agreement with the macroscopic stability measurement made by [2]. Moreover, contact position statistics and local packing fraction variability provided interesting information about the effect of the pile preparation. However, these experimental results do not permit to detect the exact nature of the contact (point, line or plane) nor to measure the force. Hence, to test the accuracy of our experimental measurement as well as to improve the contact analysis, numerical simulations will be carried out.

Figure 1 .

 1 Figure 1. A: Columns of 2.5cm, 5cm and 10cm wide hexapods made of laser sintered nylon PA2200. Hexapods can form very stable structures capable of sustaining their own weight, tilt and external loading [2]. B: Cartoon where the non-convex particles are dropped one-by-one from an overhead hooper in a cylinder of diameter 15cm. The system can be vibrated by an eccentric, driven by a DC motor which is attached to the cylinder. C: The cylinder is then carefully and slowly removed by sliding lifting it vertically. Hexapods form a stable pile of height ∼ 20cm.

Figure 2 .

 2 Figure 2. A: Slice of a density matrix in the middle of the hexapod pile. Blue is for low density (surrounding air) while red is for high density (particle nylon). B: Cartoon of a density matrix with only 2 hexapods. Matrices are 16bits 2000 × 2000 × 2000 voxels. C: Cartoon of the convolution density matrix. This matrix is used for hexapod detection by convolving it with the density matrix (B) for particle with a certain orientation: θ 1 , θ 2 , θ 3 . D: Synthesized convolved matrix. Blue high value islands show position and give orientation of the particles in the density matrix (B). See text for details. E: Cartoon of two neighbour particles. The red segment show the shortest distance between the two blue hexapod arms.

Figure 3 .

 3 Figure 3. A: Full reconstruction of an hexapod packing from the X-ray scan of a non vibrated column of 2.5cm stars. The drawing has been generated using Vapory, the Pov-Ray wrapper in Python, knowing the hexapod geometry, and the measured position and orientation. Red spheres show the position of the contacts. B: Zoom in the specified area presented of A. A complete movie of the reconstruction for each particle size is given here: https://www.youtube.com/watch?v=2HUUAPet2rs.

Figure 4 .

 4 Figure 4. A: Probability density function of the position of a contact along a particle arm for the different system configurations: 2.5cm, 5cm particles, vibrated and non-vibrated. x-axis is rescaled to compare 2.5cm and 5cm stars. Vibration has no effect on the average contact position but the longer the arms, the closer the contact points are from the particle center. B: Probability density function of the particle local density. For each configuration the probability density function of the rescaled Voronoi volume around each star is plotted. The larger the particles or the more vibrated the system, the wider the statistical repartition. C: Averaged scaled Voronoi volume along the radial (horizontal plot) and height (vertical plot) directions for the 3 different packings. The particle density just slightly increases radially. Vertical error-bars give the 95% confidence intervals of the measurements.

Acknowledgements

This work was supported by the Labex NumEv (anr-10-labx-20), the William M. Keck Foundation, by NSF grant DMR1206351 and NASA grant NNX15AD38G and ITASCA Educational Partnership Program.


Yiqiu Zhao

Jonathan Barés

Hu Zheng

Robert Behringer

Tuning strain of granular matter by basal assisted Couette shear

We present a novel Couette shear apparatus capable of generating programmable azimuthal strain inside 2D granular matter under Couette shear. The apparatus consists of 21 independently movable concentric rings and two boundary wheels with frictional racks. This makes it possible to quasistatically shear the granular matter not only from the boundaries but also from the bottom. We show that, by specifying the collective motion of wheels and rings, the apparatus successfully generates the desired strain profile inside the sample granular system, which is composed of about 2000 photoelastic disks. The motion and stress of each particle is captured by an imaging system utilizing reflective photoelasticimetry. This apparatus provides a novel method to investigate shear jamming properties of granular matter with different interior strain profiles and unlimited strain amplitudes.

Introduction

Granular matter can be jammed by shear. A jammed granular system can sustain finite stress without irreversible deformation [1]. If the packing fraction φ of a granular system is larger than a particular value φ S , the system can transition from an unjammed state to a shear jammed state after a finite amount of shear strain [1][2][3][4][5]. In simulations, the shear strain is usually defined as the affine deformation driving individual particles and can be easily tuned [2,5]. However, in the experimental case the shear strain is usually given by the deformation of the confining boundaries [1,3,6,7]. This boundary strain, although a practical tool, is not always useful to describe the deformation inside the granular matter. The shear jamming transition is achieved by percolation of a force network inside the granular system under shear [1]. The evolution of this force network is strongly related to the interior strain of the granular system [1,6]. Therefore the analysis of a stress-strain relation based on boundary strain may not be able to probe the intrinsic constitutive relation of granular matter undergoing a shear jamming transition.

To understand the interior strain of the granular matter, we present a novel apparatus that can generate boundary Couette shear with programmable interior strain field. The interior strain of the granular system corresponds to the displacements of each particle. This displacement field, or strain field will be called in this paper the resultant profile. And the collective motion of boundaries and bottom rings used to apply shear will be called the driving profile.

In this work, we compare the resultant and driving profiles for three kinds of shear: (i) boundary Couette shear (ii) linear shear and (iii) power-law shear with exponents n = 0.5 and n = 2. We show that this novel Couette shear apparatus leads to resultant profiles keeping the features of the driving profiles. In particular, we show that a linear shear results in a homogeneous density field even for very large boundary shear strains (∼ 500%). By extracting stress information from the photoelastic patterns, we track the evolution of force network and probe the shear jamming transition. This apparatus provides a novel method to investigate the role played by the interior quasistatic strain of granular matter during shear jamming transition.

Experiment

The granular system studied here consists of about 2000 frictional (friction coefficient: µ ∼ 0.6) bi-disperse photoelastic disks. The bigger ones have a radius Rb = 7.95mm, while the smaller ones have a radius Rs = 6.35mm. The size ratio is kept constant for all experiments: N big /N small = 1/3. The particles have reflective bottoms (see fig. 1(b)), making it possible to view the photoelastic patterns from the top using reflective photoelasticimetry [8]. We use a cooling system to ensure the room temperature remains the same to get rid of any possible thermal deformation of the particles. The system consists of 12 cooling fans to exhaust hot air heated by step motors and one thermal meter to monitor the temperature (fig. 1(d)). The room temperature is kept at 70 • F, at which the particles are free from any residual stress.

The particles are placed randomly onto the Couette shear apparatus shown in fig. 1. This Couette shear appara- tus is made of two boundary wheels and 21 independently movable rings. Each ring has a width that is the same as the diameter of the smaller particles. As the ring rotates, it carries the particles on it. However, the friction between particle and ring (∼ 4 × 10 -3 N) is smaller than any interparticle forces are visible as photoelastic patterns. The inner wheel and the inner-most ring of the apparatus are joined and the outer wheel moves with the out-most ring. Each ring but the inner-most one is controlled by a stepper motor beneath. Both boundary wheels have racks attached to apply boundary shear (see fig. 1(c)).

To apply shear to the system, both outer boundary and bottom rings move to drive the system. The shear is quasistatic, which means we shear the system step by step and take pictures between steps when the rings are stopped. For each step, the shear rate applied by the rings is in the order of 10 -3 s -1 . The motions of rings are controlled by 20 stepper motors. We call the relation between the number of steps a motor moves, N(r), and the radial position of the ring controlled by the motor, r, the driving profile. The displacement of the ring is N(r) • d, where d = 0.01431(9)mm is the same for all rings. The driving strain applied by this ring is defined as (r) = N(r) • d/r (see fig. 2 insert). The strain applied by the outer-most ring will be called the boundary strain.

The driving profile determines the collective motion of all the rings and boundaries. If N(r) = c • r, we call this driving profile a linear profile, and the shear applied to the system in this way a linear shear. Similarly, if N(r) = c • r n , we call the driving profile a power-law profile with exponent n, and this type of shear a power law shear with exponent n.

Each experimental run is composed of a series of quasi-static shear steps. Fig. 2 presents two driving profiles for one quasi-static steps. The driving profile is the same for all quasi-static steps of a given experimental run.

For every quasi-static steps, the motors move with different speed so that they start and stop at same time. The boundary shear displayed in fig. 2 means that only the outer-most ring and the outer boundary fixed to it are moving. This boundary shear can be regarded as a power-law shear with a very large exponent. As an example, the red curve demonstrates a fit with exponent n = 100. The imaging system of the experiment is shown in fig. 3(a). A 20Mp CCD camera (Canon R 70D) lays 3m above the Couette shear apparatus to image the whole system. Between two consecutive quasi-static steps, the camera is triggered to take two pictures. The first one is taken under UV light (see fig. 1(a)). The particles are labeled by UV sensitive marks on their top, which permits us to find their position during picture post-processing. Such a picture is presented in fig. 3(b). Another picture is taken using reflective cross-polarized light coming from the green LED light surrounding the camera. This green light picture shows the photoelastic response of each particle (see fig. 3(c)).

Result

We call the displacements of particles resulting from the driving profile a resultant profile. The resultant profiles of the driving profiles shown in fig. 2 are displayed in fig. 4. In these experiments, the granular system under shear has a packing fraction of φ = 0.7304. The resultant profiles show the displacements of particles after 398 quasi-static shear steps, corresponding to 500% boundary strain. As shown in fig. 4, the azimuthal displacements of particles generated by linear shear are very well fitted by a straight line. However, this linear fit does not exactly goes through origin. This means there is a small global rotation of the whole system. The red circles in fig. 4 show the displacement of particles under boundary shear. For azimuthal displacements, the power law-fit gives an exponent n = 16.9 (7), which is significantly less than the driv-ing profile (n = 100). But it keeps the feature of driven profile: the deformations are localized close to the outer boundary. The right panel of fig. 4 shows the radial displacements of particles. The green curve shows the average displacements from the linear shear resultant profile. It fluctuates around zero which indicates there is no collective flow along the radial direction. However, as shown by the black curve, which is for boundary shear, the average radial displacement is non-zero and positive (towards ring center). This collective radial motion results in a shear band near the outer wheel. Fig. 5 shows the radial density (packing fraction) profile of the sample granular system before and after linear shear (a) and boundary shear (c). Blue squares show the density distribution before shear while red ones show the density distribution after 500% shear strain. The data points are averaged from the local packing fractions defined as particle area divided by the Voronoi area associated with the same particle. Error bars are at a 95% confidence interval. We remark that after 500% shear strain, the boundary shear generates an obvious shear band (lower density area) with width around 10 • Rs, while the system under linear shear remain homogeneous everywhere except for small fluctuations due to the intrinsic randomness of the granular systems.

We present the green light images for the sample system after 500% boundary strain under linear shear (see fig. 6

Displacements Steps

Figure 7: Driving and resultant profiles for a granular system with packing fraction φ = 0.7257 after 125% boudnary shear strain. The left panel shows the number of motor steps for each ring needed to generate the power-law driving profiles. The right panel shows the resultant azimuthal displacements for individual particles (Unit: Rs). The power-law fits show that driving profiles with exponents 0.501(2) and 2.000(6) generate resultant profiles with exponents 0.526(1) and 1.879(9) respectively.

two key observations from fig. 6: (i) the apparatus is able to generate resultant profiles retaining the features of the driving profiles even if the inter-particle forces are much larger than basal friction forces; (ii) by changing the interior strain we generate force networks with significantly different structures under the same boundary strain. Both observations confirm the ability of our apparatus to be used to study the relation between interior strain and shear jamming, and shear flows.

Besides linear and boundary shear, we can generate arbitrary programmable driving profiles. We present two more examples. Fig. 7 presents the driving and resultant profiles for a system with packing fraction φ = 0.7257 under power-law shear with exponents n = 0.5 and n = 2. The left panel shows the driving profile for an individual quasi-static step. And the right panel shows the resultant particle azimuthal displacement. Power-law fitting of the resultant profile shows that an n = 2 driving profile results in an n = 1.879(9) resultant profile, and an n = 0.5 driving profile results in an n = 0.526(1) resultant profile. Although the resultant exponents deviate slightly from the driving exponents, they are close enough to retain the main features of the driving profiles.

Conclusion and outlook

We have described a novel experimental system that can tune the interior deformation of granular systems under boundary Couette shear with basal assistance. The shear is applied by a novel Couette shear apparatus consisting of 21 independently movable rings and two boundary wheels. We present resultant profiles for three kinds of driving profiles: (i) boundary Couette shear, (ii) linear shear, and (iii) power-law shear with exponents 0.5 and 2. The resultant profiles keep the features of the driving profiles. In particular, we showed that under linear basal assisted shear, the granular system does not form shear bands even under very large boundary strains (∼ 500%). This setup makes it feasible to investigate shear jamming properties of granular system by changing interior strains.

Figure 1 :

 1 Figure 1: Experimental setup. (a) The Couette shear apparatus with outer boundary diameter 70.9cm and inner boundary diameter 17.5cm. Concentric rings can be rotated independently. (b) The photoelastic particles. The upper particle is up side down, showing its painted reflective bottom. (c) Blue rack attached to to inner and outer boundaries to increase their roughness and apply boundary shear. (d) Cooling system. Left: fans to exhaust hot air. Right: Stepper motors controlling the motion of rings.

  •r Power Law Fit N=c•r^100 Distance To Inner Wheel: r (Unit:Rs) Number of Motor Step: N

Figure 2 :

 2 Figure 2: The driving profiles of linear and boundary shears. Main panel: The motor steps needed to generate linear and boundary shear for one quasi-static shear step. Inset: The azimuthal displacement of a ring at radial distance r is N(r) • d, where d = 0.01431(9)mm. The driving strain caused by this ring is (r) = N(r) • d/r.

Figure 3 :

 3 Figure 3: (a) Imaging system of the experiment. Particles placed on the shear apparatus can be luminated by two kinds of lights: (i) UV light from the UV LEDs surrounding the shear appratus and (ii) Green polarized light from the green LEDs surrounding the camera. (b) UV light image. Labels drawn on the particles using UV sensitive ink are luminated. Those labels are used to find centers and orientations of the particles. (c) Green light image.Reflective photoelasticity permits to observe stress inside particles. Green light is used because it is optimal for the circular polarizers we used.

Figure 4 :

 4 Figure4: Resultant profiles of linear and boundary shears after 500% boundary strain. This example system has a packing fraction of φ = 0.7304 (see also fig.2). Left panel shows the azimuthal displacements of particles v.s. their initial position. The linear fit slightly deviates from the origin. The power-law fit gives an exponent n = 16.93(34). The right panel shows the radial displacements of particles v.s. their initial positions. The positive direction is centripetal direction. Boundary shear results in a collective centripetal motion while linear shear does not.

Figure 5 :

 5 Figure 5: Shear-band study. (a) and (c) present the radial distribution of packing fractions before and after 500% boundary strain for systems under linear shear (a) and boundary shear (c). The blue squares show the density distribution along the radial direction before shearing, while the red ones display the density after shearing. The blue lines indicate the system packing fraction (φ = 0.7304). In (a) no low density area is created. But in (c) a localized low density band (shear band) is created after shear. (b) shows the local packing fraction computed from Voronoi cells after 500% linear shear. It is compared to the case where the system is sheared from boundaries in (d).

  (a)) and boundary shear (see fig.6(b)). There are

Figure 6 :

 6 Figure 6: The green light cross polarized image for system with φ = 0.7304 under linear (a) and boundary (b) driven profile after 500% boundary strain. The particle with brighter intensity bears larger force.

e-mail: yiqiu.zhao@duke.edu e-mail: jb@jonathan-bares.eu e-mail: hz64@duke.edu e-mail: bob@phy.duke.edu

Acknoledgement

This work is supported by NSF-DMR1206351, DMS1248071, NASA NNX15AD38G.


Serge Mora
email: serge.mora@umontpellier.fr

Thi-Lo Vu
email: thi-lo.vu@umontpellier.fr

Jonathan Barés

Saeid Nezamabadi
email: saeid.nezamabadi@umontpellier.fr

Highly deformed grain: from the Hertz contact limitation to a new strain field description in 2D

Most of the 2D studies on granular materials consider an assembly of cylindrical particles in the small deformation regime. Investigating the range of validity of the small deformation hypothesis is therefore of crucial importance in order to determine the quantitative relevance of the predictions emerging from these 2D studies. We introduce here a novel experimental set-up based on high resolution imaging and image correlation, and capable of measuring the strain field inside the particles with a cutting edge accuracy. We measure the strain field inside a unique circular cylindrical elastic particle undergoing an uniaxial compression. The significant deviations arising from the comparison with the predictions of the small deformation limit predicted by the Hertz's theory, and numerical simulations, give insight to improve the modeling of 2D granular materials.

Introduction

When local forces are measured or computed in granular matter physics, whether it is for numerical [1][2][3] or experimental studies [4][5][6], collections of circular grains with low polydispersity and in small deformation regime are very commonly considered. This permits to conform with the Hertz contact law [START_REF] Johnson | Contact mechanics[END_REF] which, along with the superposition principle, provides the most common mechanical description of a grain contact. Unfortunately in real life most of the granular problems involve unevenly shaped grains, on a very broad range of sizes, in highly stressed situations [START_REF] Romeo | [END_REF][9][10][11]. The lack of a net demarcation line separating the regime of small deformations with a regime where it is necessary to take into account large deformations within the grains is a source of confusion. For example in the case of jammed granular systems some grains (the one involved in force chains [1]) can be far above the small deformation hypothesis which ruins the macroscopic description of the whole system.

Hence, in order to observe, model and simulate properly realistic granular media it is necessary to get a complete description of a grain behaviour beyond the linear elastic limit. This would permit to check already existing numerical models [12][START_REF] Topin | selection of Papers from the Symposium Powder Science and Technology -Powders and Sintered Material STP-PMF 2007 Symposium on Powder Science and Technology -Powders and Sintered Material STP-PMF 2007[END_REF][START_REF] Nezamabadi | [END_REF] and to provide inputs for improvements.

By means of novel 2D experimental approaches, we investigate when the Hertz strain field [START_REF] Johnson | Contact mechanics[END_REF] is not valid any more and how it has to be modified beyond linear elasticity. To do so we introduce a new experimental measure-ment method suitable for granular material. We set-up a step by step single disk radial compression experiment attached to a very high resolution imaging apparatus (5µm accuracy over a 6cm diameter disc particle). The disk is one base surface of a thin cylinder made in an hyperelastic material. The displacement field in the plane of the disk is then precisely measured both spatially and in intensity via Digital Image Correlation (DIC). Strain is directly deduced and compared with theoretical results. The limitations arising from the linear elastic theory are then highlighted, for this common mechanical situation. These results are also compared to what is obtained by finite element (FE) simulations.

In the first part of this article, we present the experimental set-up, how we compress a particle and how we image it. Then, in a second part we give some image processing details to measure the displacement field inside the particle. Strain field is computed and compared to what is obtained with FE. Finally, in a last part, we present some results, conclusions and perspectives.

Experimental set-up and material

A manufactured in-house black silicone cylindrical particle [START_REF] Smooth-On | Mold Star R 15 SLOW colored with few drop of black Silc Pig R[END_REF] (Young modulus 2.7MPa) coated on one side with shiny very thin glitter [START_REF]Smooth-On R , Cast Magic R (silver bullet)[END_REF] (average size 30µm) is led on the top glass of a flatbed scanner as presented in fig. 1-A. The initial diameter of the particle is 5.9cm and its thickness is 9.5mm. As presented in fig. 1-B the glitter coating forms an homogeneous thin random pattern on the bottom face of the particle with a high black/white contrast. Step by step, the flat cylinder is squeezed in its radial direction in a quasistatic manner by a manufactured in-house uniaxial compression apparatus (see fig. 1-A). This device Figure 1. A: Sketch of the experimental set-up. A cylindrical particle (5.9 cm diameter and 9.5 mm thickness) lays on a flatbed scanner. A manufactured in-house compression device squeezes the particle step by step while the particle is imaged from below and the compression force is measured. B: Black and white scanner bottom view of the particle. A thin pattern has been deposited on the particle. C: Zoomed view of the particle pattern. D: Zoomed view to the size of a correlation length. The pixel size is 5.29 µm.

is accurately actuated by means of a stepper motor and a threaded rod with a very low distance/step ratio. The force is measured with two 1kg force cells. Vegetable oil with a very low viscous coefficient (60 mPa.s) is coated on the scanner glass to almost vanish static basal friction and improve optic transmission.

After each 0.5mm compression step, the system is imaged with the flatbed scanner. As emphasized by the results presented in this paper, it constitutes a very stable and accurate measurement apparatus both in space (tunable from 70dpi to 9600dpi) and numerical depth (tunable from 8bits to 16bits for each color) on a large area (maximum 21 × 29.7cm 2 ) without any optical distortion on a flat area. For the results presented in this paper, we used a scanner CanoScan 9000F Mark II with resolution 4800 × 4800dpi 2 on a 6 × 6cm 2 area on 8bits for black and white picture. As presented in fig. 1-C & D, we image the particle bottom on a wide area with an extreme accuracy.

Image post-processing

Image correlation for granular material has been widely used both in 2D and 3D at the continuum scale [START_REF] Hall | [END_REF][18][19][20][21] to measure the grain flow but to our knowledge it has never been used at the scale of each grain to measure their deformation. However, as emphasized earlier in this paper it is of high interest to study stress networks in the most common granular media: unevenly shaped, highly stressed or deformed grains with any rheological behaviour. Until now, the only way to measure force between grains was to use photoelasticity [4,22], thermoelasticity [5] or boundary deformation [6]. However for each of these techniques, particle materials and shapes are very restricted and it can be performed only in the small deformation regime. Here we propose a method capable of measuring the full stress field for any kind of grain and any loading in 2D.

By means of digital image correlation (DIC) technique [23,24], adapted to high image size and large deformations we measure the evolution of the displacement field inside a compressed silicon particle. A regular mesh of 260 × 260µm 2 (50 × 50px 2 ) is first built on the undeformed initial picture. For each cell of the mesh, also called correlation cell, the displacement along x and y directions is measured incrementally between the picture of each compression step and the undeformed picture. To do so, a first measurement is done, with a 1px accuracy, using a Fourier transformed based algorithm [24]. Then, the measurement is improved to subpixel accuracy by means of computationally heavier algorithm based on an optimization algorithm [25] maximizing the correlation between deformed pictures and the undeformed one. To take the large deformations into account both algorithmic steps are performed sequentially for correlation cells going from 3 times bigger to 50 × 50px 2 . Parallelized on 12, 3GHz nodes, the computation take around ∼ 6 hours for 30 compression steps.

Results and discussions

As shown in fig. 2-A & B for a global compression of 8% strain, the x and y components of the displacement, u x and u y (y is for the direction of the global compression, x is for the in-plane orthogonal direction) are measured with Using symmetries only 1/4 of the the particle is simulated. The material is supposed to be linear, elastic and isotropic for a global 8% compression. Fields are plotted on the deformed system (Eulerian). Results are qualitatively similar to what is found experimentally. Difference is explained because the system is at the limit of the small deformation hypothesis.

an excellent smoothness and accuracy compared to what is classically done in the small deformation domain for disc compression [23]. Nevertheless around the narrow particle boundary, data are noisy since correlation cells overlap both particle and unpatterned surrounding space: correlation information in the particle is mixed with surrounding black area where no correlation is possible. Hence for these correlation cells displacement information is less accurate than displacement for cells in the middle of the particle. Deriving these raw displacement data u x and u y , as shown in fig. 2-C & D we then compute the strain inside the particle, ε xy = ∂u x ∂y + ∂u y ∂x and ε yy = ∂u y ∂y respectively. It is important to note that these data are obtained without any smoothing nor filtering at a scale higher than the cell scale. As expected from Hertz's theory [START_REF] Johnson | Contact mechanics[END_REF], we find (from ε yy ) that the strain parallel to the compression direction is symmetric and maximum in a thin ellipse joining the two contact points. In fig. 2-D we also observed slight horizontal changes in the blue shade. They correspond with small instability of the scanner motor speed. This is the only artefact we have noticed coming from the scanner and it can be put in perspective to the fact that no noise filtering has been applied to the pictures nor to the post-processed data. Similarly, for ε xy , results qualitatively in agreement with Hertz's theory are observed since the shear strain is anti-symmetric with two opposite lobes going from the contact areas. Then, assuming a linear rheology, the stress field is computed from the measured strain field with a Young modulus of 2.7MPa and a Poisson ratio of 0.49 (not shown).

In parallel to this, numerical simulations have been carried out, by means of the FE code Cast3M [START_REF]Cast3[END_REF]. Taking advantage of the symmetries, 1/4 of the particle is simulated in the same mechanical condition as the experiments except the material is considered as linear, elastic and isotropic for the full deformation range. As presented in fig. 3-A & B the results are in qualitative agreement with what is measured experimentally but strain intensity do not seem to be the same.

In fig. 4-A, we show ε yy (y) the compression strain value along the loading direction. This has been plotted for different compression steps. We see that this value vanishes at the contact point, rapidly increases toward the center of the particle and then decreases before reaching a plateau in the middle of the sample. As presented in fig. 4-B, this is consistent with what is predicted by Hertz contact law for small deformations. Also, even for large deformations the minimum value of the strain near the contact points is well predicted. However when the global loading is larger than few percent, fig. 4-AB shows that the strain is quantitatively different in the center of the disc: the strain is largely overestimated. This can be explained by the fact that small deformation hypothesis is not valid any more or because out of plan displacement cannot be neglected any more in the experiment.

Conclusion and perspectives

In conclusion, in this paper we have presented an original experimental apparatus capable of compressing cylindrical granular particles while imaging them with an excellent accuracy. Along with a manufactured in-house fast DIC algorithm adapted to take large deformation into account, the evolution of the displacement field inside the particle is measured with an edge cutting accuracy. Strain field has been deduced and compared to what is expected from Hertz's contact law. A fairly good qualitative agreement have been observed as long as the system stays in the small deformation regime.

This paper validates an experimental method for the measurement of the full strain and stress field at the grain scale for any loading and opens a new front of experimental exploration for a wide range of granular material: soft, brittle, mixed, etc. Hence in future works, the experimental measurement will be extended to different model materials: plastically deformable material, brittle material, zero Poisson ratio material,etc. The results will be compared to the ones obtained with different numerical models (finite element method for example) for different boundary conditions. The goal will be to show the limitation of these numerical models either by the fact they do not take the full complex rheology into account or by the fact they are limited in terms of spacial discretisation of the particles (for the bulk or for the contact model). This method can also be extended to 3D measurements in granular matter using tomography imaging.

Figure 2 .

 2 Figure 2. A,B: Displacement field measured from the local image correlation along the x and y direction respectively. C, D: ε xy and ε yy components (respectively) of the strain field. Fields are plotted on the undeformed geometry for a 8% total deformation of the particle. Measurement are noisy at the boundary because the correlation cells overlap patterned and unpatterned domains. The different fields are plotted on the undeformed coordinates (Lagrangian).

Figure 3 .

 3 Figure 3. Computed strain.A, B: ε xy and ε yy components (respectively) of the strain field computed by FE using Cast3M[START_REF]Cast3[END_REF]. Using symmetries only 1/4 of the the particle is simulated. The material is supposed to be linear, elastic and isotropic for a global 8% compression. Fields are plotted on the deformed system (Eulerian). Results are qualitatively similar to what is found experimentally. Difference is explained because the system is at the limit of the small deformation hypothesis.

Figure 4 .

 4 Figure 4. Strain evolution during the compression process along the y-direction. A,B: Evolution of the strain component ε yy (y) as function y, along the diameter parallel to y, for different compression levels going from 0% to 10% in the experimental case (A) and computed from Hertz's contact law (B). y = 0 corresponds to the center of the specimen. For the experiment, data close to the contact point have been cropped because they are noisy. The curve shape is sensibly different for large deformations.

Acknowledgements and References

This work is supported by the Labex NumEv (anr-10-labx-20).