Particle scale force sensor based on intensity gradient method in granular photoelastic experiments

The intensity gradient method (G2 method), namely computing the light intensity gradient-squared, is a widely used non-invasive experimental method to extract stress information from quasi-two-dimensional photoelastic granular materials. Previous works show that calibrated G2 is an accurate measure of global stress. However, whether it can be used at the particle scale aside from the special case of diametric loading remains unclear. We test here the applicability and limitations of G2 as particle scale stress indicator and specify its dependence on relevant experimental parameters of the particles, light conditions, and imaging system. We first propose an explicit formula to calculate the relationship between the G2 value and stress based on the linear elasticity and photoelasticity theories, and then validate our formula by numerical and experimental tests. We find that G2 is proportional to ∑ i ∣ F ⃗ i ∣ , the sum of magnitudes of the contact forces, for disc particles when forces are not large. We also observe that, for large enough resolution, G2 does not change with the number of contacts as well as the direction of the contact forces under same ∑ i ∣ F ⃗ i ∣ value. However, we find that this relation between G2 and ∑ i ∣ F ⃗ i ∣ is not universal for any particle shape. As an example, we show that a square particle can have dramatically different values of G2 under the same contact forces with different contact types (point-edge contact and edge–edge contact).


Introduction
The transition between jammed and unjammed states of granular materials is very important in many industrial processes [1] and natural hazards [2,3]. Observing and studying the microscopic, i.e. particle-scale, stress is crucial to understand the macroscopic transition in granular materials. The challenge is how to experimentally extract particle scale stress information non-invasively. Currently, there are several experimental techniques up to this task, including X refractive ray/computed tomography scan [4][5][6], reflectively index matching [7], and photoelastic imaging [8][9][10].
The photoelastic technique has been widely used in granular materials research, and is one of the most effective experimental methods, for two-dimensional (2D) granular systems [11][12][13]. Two prevalent methods well developed in Behringer's lab [11,14] can be used to extract the particle scale stress from the experimental data: the force resolving method [9,11] and the G 2 method [14]. The force resolving method is based on a nonlinear fitting algorithm that can reveal the contact vector forces in quasi-2D granular materials composed of discs with different sizes [11]. Details of this nonlinear fitting algorithm have been well described in a recent review paper [10]. However, the implementation of this technique depends on the accuracy of the initial guess for the contact forces before fitting [9], and is highly computationally consuming [10]. By contrast, while unable to solve the contact vector forces of each particle, the G 2 method can well reveal the particle scale pressure/stress in a reasonable stress range, which is first reported by Howell et al [14]. The G 2 method is still widely used because of its computational efficiency and simplicity [8,10,12,[15][16][17][18][19][20][21][22][23]. G 2 is defined by the averaged square of the intensity gradient of the region of interest in the experimental image recording the photoelastic fringes, and monotonically grows with the boundary stress of the granular packing [8,14,15,22]. Even though several recent works assumed that a similar relation holds at the grain scale [18,21], the particle scale validity of this empirical method has been neither systematically tested nor proven except for diametric loadings.
In this work, we try to answer questions including what stress information is correlated with G 2 at the grain scale, what determines this correlation, and whether the correlation is the same for different particle shapes. We first propose a general formula that specifies the dependence of G 2 on all the relevant experimental parameters: (1) the background light intensity I 0 , (2) the stress-optic coefficient f σ , (3) the camera resolution N, defined as the number of pixels per meter, and (4) the disc particle radius R. Then we provide both numerical and experimental tests to validate the formula. Finally we consider the validation of the formula on non-circular particles. Our work provides regulations and information to setup and analyze granular experiments using photoelastic materials.

Theory
For an experimental image showing fringes of photoelastic particles under stress, G 2 for a particular disc is defined as å å º -+ - I  I  I  I   I  I  I  I   1  4  2  2 2   2  2  which is the discrete version of calculating the averaged inner product of the intensity gradient at different pixels inside the disc. I i j , in equation (1) is the intensity of light at (i, j) pixel [9,14]. More information on the set-up of the lightening and imaging system utilizing photoelastic technique can be found in [10]. Viewed with an ideal polariscope, the light intensity is determined by Ct , where λ is the wavelength of the illumination light, t is the height of the particle, and C is a material constant, (2) the background light intensity I 0 , and (3) the eigenvalues of the stress tensor ŝ (denoted as σ 1 and σ 2 ) [9]: where s 1 and σ 2 depend on contact forces acting on the particle. We consider a disc under Z number of contact forces, each contact force is specified by three parameters b a ) 0, 2 is the impact angle, and  | | F is the magnitude of the contact force  F (figure 1). In order to get the distribution of light intensity inside the disc, we need to calculate the distribution of the stress tensor. We show a brief derivation below and a more completed version can be found in [9]. In the following derivation, linear elasticity is assumed, and the results are applicable only for small strain, which is the case for typical experiments in granular physics where the deformations of particles are usually small (e.g. [9,20]).
To calculate the stress tensor caused by one of the Z forces (denoted as  F ), we consider a spatial position inside the disc (c, r) (figure 1). The infinitely large plate solution gives the only non-zero component of stress where r real is the distance in meter from the contact point to the point (c, r). The pixel unit distance is obtained by = r Nr dis real , where N is the number of pixels per meter. Also note here we imply a polar coordinate system, and the base vectors at point (c, r) is shown by the small red arrows in figure 1. In order to rotate this axis by angle a a q ¢ º + to the Cartesian axis system, we apply the rotation matrix a¢ where R is the particle radius in meter. Due to the geometry, the complete stress tensor at (c, r) is therefore For multi-contacts, the resultant stress tensor at point (c, r) is Note that in the above derivation we assumed the strain is infinitesimal. Moreover, the model does not assume any relation between contact force magnitude and strain at contact (e.g. Hertzian or linear). From the above equation we see each element of the stress tensor is a linear combination of the contact force magnitudes divided by a length in meter. We define ¢  F i for ith contact force as ¢ º at a particular pixel position. Note any parameter inside g function except ¢  F s i is either dimensionless or has unit of pixel length. We derive a general formula to calculate G 2 below: The last step uses , which is an assumption for small enough forces. We note in principle c in equation (8) does not have to be constant, although we assumed it to be a constant for simplicity and we will test this assumption later. Equation (8) provides a general formula that relates G 2 with contact forces of a particular disc particle using all the relative experimental parameters I 0 , N, R and f σ . In the following sections we test if this relation really holds. Since I 0 dependence in equation (8) is straightforward from the derivation, we will focus on N, R and f σ in this paper.

Different contact number
Previous experiments confirmed that under diametric loading (Z = 2 and both forces are normal), is the pressure acting on a disc particle [18,20]. It is more convenient to study this relation for different number of contacts using numerical simulation. For simplicity we first simulate a particle under Z normal contact forces that have the same magnitude and are distributed evenly over the circular boundary. The simulated disc has f σ =100, R=10 −2 m, and I 0 =1. The values of f σ and R are close to the values actually used in typical experiments (e.g. [9,10,12]). Figure 2(a) shows how G 2 changes with pressure P (and å It is clear that the linear dependence holds for different Z, but the slope of G 2 drops a little when Z increases. We attribute this discrepancy to the finite resolution N. To show this, we consider the relative difference at P=100 N m −1 between the values of G 2 with Z=2 and Z=10, quantified by 1 . The relation between Δ and N with other parameters fixed is shown in figure 2(b), which shows a power law decay without an upper cutoff as the resolution increases. So we expect that when the resolution goes to infinity, all Z cases will have the same linear behavior with the same slope. This fact is even more striking when comparing their very different photoelastic patterns (figure 2(c)).

Different tangential forces
We then test if the linear dependence holds when the contact forces are not normal, which is common for frictional granular materials. For a disc with friction coefficient μ, the largest possible contact angle β=arctan (μ) is referred to as the friction angle. We consider a pair of forces with the same force magnitude  | | F acting at the friction angle for different μ. An example is shown in the inset of figure 3(b). Figure 3(a) plots the G 2 changes with  | | F at frictional angle β(μ) for different μ. For small forces, different μ seems to have similar linear slopes. However, the difference between different μ does not vanish as the resolution goes to infinity. The inset of figure 3(a) shows the relative difference Δ of G 2 at F=2 N between μ=0 and μ=2 changes with the camera . The relative difference Δ goes to a finite value, Δ≈10%, as  ¥ N , but this difference is not unacceptable experimentally. Still, the discrepancy of G 2 at the same pressure becomes very large ( figure 3(b)).The above observations hold when the contact numbers change. For example, we consider another case when Z=4 which includes 2 pairs of forces with same magnitude  | | F that all act on the friction angle (shown in figure 3(d) inset). Figure 3(c) plots the G 2 versus å  | | F i i for Z=2 cases (") and the Z=4 cases (!). Figure 3(c) shows that the linear relation still holds for different Z when å  | | F i i is not large. Figure 3(d) plots the same G 2 data versus pressure, which shows the linear relation does not change with contact number even for large tangential forces. But with large μ the slope between G 2 and pressure can be very different. Particles used in most previous experiments have friction coefficients less than 1 (e.g. [12,15,18]).

A real granular packing
In the previous section we show that in special cases the relationship between G 2 and å  | | F i i is linear with a universal slope. Now we test in a general granular packing, where particles can have random contact numbers, positions and forces. Figure 4(a) shows a reconstructed image using the fitted forces extracted in a real jammed . Δ goes to finite value as  ¥ N .  The system contains particles with two different radii R big =0.56 cm and R small =0.44 cm. Each disc is in force and torque balance. The photoelastic fringes are calculated using equation (2). The friction coefficient between particles is μ≈1. (b) The relation between rescaled G R 2 2 and pressure for each individual particle in the jammed granular system. R is the radius of the particle being calculated. (c) The relation between rescaled G R 2 2 and å  | | F i i for each individual particle in the jammed granular system. The data collapses much better than (b). granular system [9]. The system contains particles with two different radii: R big =0.56 cm and R small =0.44 cm. Each disc is in force and torque balance. The photoelastic fringes are calculated using equation (2). The system contains about 1000 particles. The friction coefficient between particles is μ≈1. We check the dependence of G 2 on å  | | F i i and pressure respectively for those particles. Figure 4(b) plots the rescaled G R 2 2 versus pressure for each individual particle in the jammed granular system, where R is particle radius. The data are scattered as expected from the simulations we showed above. Figure 4(c) plots the relationship between rescaled G R 2 2 versus å  | | F i i for each individual particle in the jammed granular system. The data collapses nicely on a straight line. This test confirms that G 2 always has a linear dependence on å  | | F i i under loading situations with any contacts.

Test of the general formulation
Now we have shown that with high resolution and small forces, equation (8) holds for different Z and different contact angles with same R, N and f σ values. We then move forward to test the validation of equation (8) when R, N and f σ changes. We study how the linear slope = ¶ ¶ å  ( | |) k G F i i 2 changes with f σ , N, and R. Figure 5 for example shows our analysis of G 2 when changing f σ while keeping all other relevant parameters I 0 , N, R the same. Figure 5(a) plots G 2 versus å  | | F i i using different f σ . The black lines show linear fits for the linear regime with different f σ . The red dots mark the points where the linear regime breaks down. The corresponding å  | | F i i value is denoted as the saturation force F S . Figures 5(b) and (c) plot the change of the linear slope k and the saturation force F S with f σ . Both figures show good power-law dependence as expected from equation (8). We perform same analysis on N and R as well, which all follow the expected scalings from equation (8).
The overall validation of equation (8) can be proved by looking at the slope, k, for different f σ , N, R values changing together over large ranges, as shown in figure 6(a). Here the background light intensity is set to be I 0 =1. Figure 6 , which confirms that equation (8) holds for the numerical simulations. We note that the range of values we tested in numerical simulation for parameters R, N, and f σ is not only large (over 3 decades for each) but also covers the typical values the most up-to-date experiment may have. For example, recent works typically involve ∼10 3 number of particles with R∼10 −2 m recorded by camera having~10 px 6 2 resolution [10,[18][19][20]. This typical setting gives N∼10 4 . Figure 6(b) plots the change of saturation force F S against the same set of parameters. Unexpectedly, F S scales as 1/k. This information is very important when designing and analyzing the photoelastic granular physics experiments.

Comparing to the numerical results
We have shown by numerical simulations that in ideal cases equation (8) holds for discs under point contacts in general. We now test if equation (8) holds for real experimental photoelastic discs. Figure 7(a) plots a comparison between the G 2 versus å  | | F i i for both numerical and experimental discs. Both numerical and experimental tests have all parameters I 0 , N, R, f σ and the loading conditions (Z = 2 with both forces normal and balanced) the same. We see the numerical solution gives a much larger G 2 at the same value of å  | | F i i . This may be due to the non-ideal response of the material close to the contacts. Figures 7(c) and (d) plots the color-map of the calculated G 2 distribution from numerical solution and experimental test respectively. Red means large G 2 value and blue means small G 2 value. It is clear that the value of G 2 around contact points in the real experiment is significantly smaller than the numerical solution, which confirms our argument. However, the averaged intensity á ñ I , where áñ means averaging over the pixels inside the disc, from both simulation and experiment are similar as shown in figure 7(b). It should be noted that the monotonic regime of á ñ I versus å  | | F i i is significantly smaller than G 2 . Thus á ñ I is not a good measure of stress unless the contact forces are very small. Figures 7(e)-(g) plot the intensity distribution of the same particle under point loading, plane loading and from numerical simulation under the same pair of forces with

Validation of the scaling relation
To test the scaling relation equation (8), we perform diametric loading to discs with 4 different radii: 0.92, 1.25, 1.57, and 2.26 cm under the same light conditions. Different N is achieved by blurring the experimental image by resizing. Figure 8(a) plots the G 2 dependence on å  | | F i i with different parameters. All the curves collapse in their linear regime after rescaling using equation (8), as shown in figure 8(b). Therefore although the value of c in equation (8) is different in experiment due to the non-ideal response of the material close to the contacts, the scaling relations still holds, which is very useful in designing and analyzing experiments.

General contact situations
The above analysis confirms that equation (8) holds for experimental diametric loadings. We expect this relation holds for more general cases: particles with different contact numbers, positions, and forces. To confirm this expectation, we show a test with a real jammed experimental system. The system is the same packing shown in figure 4. Figure 9(a) shows the experimental photoelastic fringes captured by camera for this packing. The white circles are plotted to highlight the particle boundaries. Figures 9(b) and (c) plot the relation between the rescaled G R 2 2 versus pressure and å  | | F i i respectively. å  | | F i i and pressure are obtained by the vector force fitting algorithm [9,10]. The data points in figures 9(b) and (c) are both scatted around a straight line. Due to the scattering of the data, no measurable advantage is observed for å  | | F i i over pressure. This is probably because the force fitting algorithm to measure å  | | F i i itself has around 5%-10% uncertainty [9], which overwhelms the small distinction between pressure and å  | | F i i shown in figure 4. Thus experimentally, for this system, G 2 may still be used as a pressure indicator considering the experimental noise. However, we expect that when the friction coefficient μ is large enough, G 2 in experiment will no longer be a good indicator of pressure, but still a good indicator of å  | | F i i . We note in this packing, μ≈1. Previous experiments have usually had μ<1 [12,15,18].
Most importantly, figure 9 shows G 2 for particles with different contact situations collapse reasonably well after proper rescaling. Thus equation (8) is verified in general cases experimentally.

Other particle shapes: contact area matters
Particles in real world granular systems and industrial products are not spherical. It is important to test if the G 2 method can be used to analyze non-circular particle systems.
Consider a square particle under a uniform edge-to-edge uni-axial compression along y direction. The stress tensor everywhere inside the square should be the same: , which gives  = I 0. Therefore, although the boundary force experienced by this square particle is not zero, the value of G 2 is zero. Consider if the same forces are applied through point contacts; then the stress field would be significantly different (similar to disc cases). To verify these theoretical arguments, we experimentally compress a square photoelastic particle using two point contact forces with the same magnitude (figure 10(c)) and using two edge contact forces resulting in a uniform stress field (figure 10(d)). Figure 10(a) plots G 2 versus å  | | F i i for point contacts (red circles) and edge contacts (blue circles). In the two cases, the value of G 2 is significantly different even at the same å  | | F i i . For completeness we calculated the averaged intensity á ñ I versus å  | | F i i in figure 10(b), which also shows a significant discrepancy between point and edge contact loading cases. Thus, G 2 in angular particles, where edge to edge contacts are possible, can not be used to determine the forces on the particle. However, we note that calibrations from previous work show that the proportionality of G 2 and boundary stress holds at the global scale even for such particles, which may due to the fact that in stable packings the ratio between edge-edge and pointedge contacts is more or less consistent [8].

Conclusion
We provide numerical and experimental tests for G 2 as a particle-scale stress measurement. For disc particles, numerical tests show that G 2 is proportional to the sum of contact force magnitude å  | | F i i , instead of the pressure experienced by the particle. The proportionality between G 2 and pressure only roughly holds when tangential forces are not large. We also show that the dependence between G 2 and å  | | F i i is independent of the number of contacts as well as the directions of the contact forces when resolution is high and forces are small. We describe Figure 9. (a) The photoelastic fringes of an experimental jammed packing captured by camera (same packing as shown in figure 4). The white circles are plotted to highlighted the particle boundaries. The system contains particles with two different radii R big =0.56 cm and R small =0.44 cm. The friction coefficient between particles is μ≈1. (b) The relation between rescaled G R 2 2 calculated from the experimental image and pressure calculated by fitted contact forces for each individual particle in the jammed granular system. (c) The relation between rescaled G R 2 2 calculated from experimental image and å  | | F i i obtained from the fitted forces for each individual particle in the jammed granular system. The data collapses slightly better than (b), but no significant improvement can be evidenced. . This formula is very useful to guide the design of photoelastic granular experiments. It is particular useful in experiments with poly-disperse systems, where multiple R value exist in one image, as well as in comparing experiments with different imaging systems, where I 0 and N are usually different. The value of c obtained from experiments is much smaller than that obtained from simulations. Since the meaning of c value is not clear, we suggest readers to perform diametric calibrations to measure c in real experiments instead of just using the value we reported in this work. Nevertheless, the scaling relation between G 2 and other parameters still holds in experimental tests. Finally, using a square particle as an example, we show that G 2 does not work as particle scale stress indicator for particles with edge contacts. Our work for the first time gives a complete analysis for an important experimental technique in granular physics, which helps to guide future experiments.