Stress and force measurement uncertainties in 3D granular materials

We have developed and employed a 3D particle stress tensor and contact force inference technique that employs synchrotron X-ray tomography and di↵raction with an optimization algorithm. We have used this technique to study stress and force heterogeneity, particle fracture mechanics, contact-level energy dissipation, and the origin of wave phenomena in 3D granular media for the past five years. Here, we review the technique, describe experimental and numerical sources of uncertainty, and use experimental data and discrete element method simulations to study the method’s accuracy. We find that inferred forces in the strong force network of a 3D granular material are accurately determined even in the presence of noisy stress measurements.


Introduction
Determining inter-particle forces in 3D, opaque granular materials has been a major pursuit in soft matter physics for decades. Techniques for determining interparticle forces include those leveraging photoelastic discs [1], compliant grains imaged with X-ray tomography [2], intra-particle speckle patterns combined with digital image correlation [3], and sand grains analyzed with X-ray computed tomography (XRCT) and 3D X-ray di↵raction (3DXRD) [4]. The combination of XRCT and 3DXRD has been used in natural and synthetic particle packs to study stress distributions [5], particle fracture mechanics [6], energy dissipation [7], and wave propagation [8]. Although prior work has described sources of experimental uncertainties in stresses and inferred forces [9,10], limited work has systematically examined uncertainty propagation.
In this paper, we discuss stress and inferred force uncertainties in 3D granular materials studied with XRCT and 3DXRD. Uncertainties arise from (1) experimental noise and (2) image processing and algorithmic uncertainties. The first type of uncertainty arises from the resolution limits of X-ray di↵raction detectors and associated raw data processing algorithms. These uncertainties cause "noisy" particle strain and stress tensors, as well as missing particle stress tensors, and are the main focus of this contribution. The second type of uncertainty arises from algorithms employed in image processing, such as particle segmentation, and data processing. These uncertainties bias inter-particle contact point locations and orientations and may lead to the presence of false contacts or absence of true contacts in the force inference procedure. These latter uncertainties will be addressed in future work. ⇤ e-mail: rhurley6@jhu.edu

Particle Strain and Stress
We have employed 3DXRD measurements at various synchrotron radiation facilities, including the European Synchrotron Radiation Facility (ESRF), the Cornell High Energy Synchrotron Source (CHESS), and the Advanced Photon Source (APS), each with unique X-ray area detectors including the GE 41RT (APS, ESRF, CHESS) and Dexela 2923 detectors (CHESS). Distinct data processing algorithms including ImageD11 [11], heXRD [12], and MIDAS [13] are used at these beamlines (ESRF, CHESS, APS, respectively). In this section, we discuss possible sources and magnitudes of noise in individual particles' strain and stress tensor components.
To estimate the noise in each particle's average stress tensor, which is employed in calculating inter-particle forces, we first identify noise in particle strain tensors. Uncertainties in particle strain tensors arise from the finite resolution of Bragg peaks on X-ray area detectors and the propagation of peak position and magnitude uncertainties through data processing algorithms (ImageD11, heXRD, MIDAS). Strain tensor errors have been well characterized as close to 1x10 4 for the detector geometries used in our past work [14,15]. We previously verified, by examining the standard deviation of strain tensor components in a packing of 360 single-crystal alpha-quartz particles that was not under external load (i.e., any apparent strains were those due to measurement errors) [10], that standard errors in diagonal and o↵-diagonal strain tensor components are 1x10 4 and 5x10 5 , respectively. We further studied standard errors in stress tensor components of the same single-crystal alpha-quartz particles by employing the known uncertainties in single-crystal sti↵ness tensor components, C i jkl , of synthetic alpha-quartz [16]. By computing the standard deviation of 400 possible stress tensors per particle from these sti↵ness tensor uncertainties, we found standard errors around 10 MPa and 5 MPa for on-and o↵-diagonal stress tensor components [10], as Table 1. Strain and stress tensor standard deviations, , for quartz (SiO2) and sapphire (Al2O3), as described in text. Strain in units of 10 4 and stress in units of (MPa). While the stress tensor component standard deviations reported in Table 1 are often half as high as the mean stress level in the samples studied in [4,10], we note that the stresses experienced by particles within the strong force network tend to be much larger than the mean. We show in section 3 that when the magnitude of noise in the stress measurements is half of or one times the mean stress level in the sample, the strong force network is minimally affected by measurement noise. These strong forces are the most salient for examining heterogeneity, fracture, and other phenomena [4,10].
Here, we repeat the procedure used to examine strain and stress tensor noise in single-crystal sapphire data described in [7]. We compute the standard deviation of particle strain tensors for more than 1,700 strain-free (before mechanical loading) single-crystal sapphire particles using identical nominal lattice parameters for each particle. The resulting strain tensor standard deviations are close to those obtained for alpha-quartz (Table 1). Standard deviations for sti↵ness tensor components of sapphire are not available in [17]. We therefore examine the standard deviation in the stress tensor components in the first load step after beginning to strain the samples. Each particle's strain tensor in the strain-free state is subtracted from its strain tensor computed after sample compression to eliminate natural di↵erences in lattice parameters between crystals (e.g., due to defects). Therefore, stress tensor standard deviations in the first load step should provide insight into the scatter of stresses caused by noise in the data acquisition (in addition to some scatter due to natural sample heterogeneity). The resulting stress tensor component standard deviations are around 50 MPa and 23 MPa for on-and o↵-diagonal components, respectively, as shown in Table  1. These stresses are again similar to what would be obtained from directly multiplying nominal sti↵ness tensor components with standard deviations of strain tensor components. We note that these uncertainties are very close to quoted uncertainties from the developers of 3DXRD software [14,15].

Force Uncertainty
We now examine the errors in inter-particle forces obtained through the optimization procedure proposed in [4]. Rather than study a specific experiment presented in previous papers, we study error propagation by using synthetic discrete element method (DEM) data for which forces are known exactly. Particle stresses are calculated directly from inter-particle forces in DEM. Therefore, when the average coordination number is approximately 8 or less, forces are exactly determinable from DEM [4,9]. We evaluate the propagation of errors as a function of the ratio of stress tensor noise (added to the exact stress tensors) to a macroscopic stress level. We first provide a brief review of the force inference technique. We note that the use of DEM restricts our attention to spherical particles. However, particle stress measurements from 3DXRD and force inference are each independent of particle shape and rather a function of particle volume. We therefore expect that while the nature of force transmission may change in packings of angular particles, the results of our analyses will not.

Force Inference Technique
The force inference technique involves three governing equations: linear momentum balance, angular momentum balance, and a force-stress relation. For a particle ↵ with N ↵ c contacts with neighboring particles or system boundaries, momentum balance is given by and the force-stress relation is given by where x (i) ↵ is the contact point, f (i) ↵ the force vector, ⌦ the tensor product, V ↵ the particle volume, and ↵ its stress tensor. Equation 2 is derived using the divergence theorem applied to the volume-average stress expression for a single particle. Combining Eq. 1 into matrix form K eq f = 0, and Eq. 2 into matrix form K st f = b st , where K eq and K st contain only components of x (i) ↵ , f contains all force vectors, and b st contains V ↵ ↵ , we minimize where is an optimal tradeo↵ parameter, chosen to weight stress and equilibrium to an optimal extent to minimize the cost function. Constraints to enforce f (i) ↵,t < µf (i) ↵,n (Coulomb friction) and f (i) ↵,n < 0, where f (i) ↵,t is tangential force, f (i) ↵,n is normal force, and µ is inter-particle friction coe cient, are often enforced. We perform minimization in Matlab using CVX [18]. The stress tensor errors considered in this work a↵ect the first of the two objective functions in Eq. 3.

DEM Model and Force Inference Without Noise
We simulate the compression of a 594 monodisperse spheres in a cylindrical geometry. The cylindrical particle packing was simulated in LIGGGHTS [19] using 1.4 mm diameter particles with Young's modulus 1 GPa, Poisson's ratio 0.2, inter-particle friction coe cient 0.4, restitution coe cient 0.1, and a "hertz tangential history" contact model [20]. Figure 1(a) shows the stress-strain curve and coordination number for the DEM simulation. Insets illustrate the distribution of zz at steps labeled 2 and 3. Figure 1(b) illustrates the distribution of zz at steps 2 and 3, illustrating the spread of particle stresses. Figures 3(a) show the interparticle forces at steps 2 and 3. Force vectors are rendered as lines, centered at contact points, oriented parallel to the force vector, and scaled in width and length linearly with contact force magnitude. These forces are employed to determine particle stress tensors using Eq. 2. Stress tensor, particle volume, and contact location information is then used in the inter-particle force inference using Eq. 3. The number of unknowns is the number of contacts times three (5,937 for step 2 and 6,435 for step 3) and the number of equations is 12 times the number of particles (7,128).
Force inference without noise results in a nearly exact match with the actual forces. A small number of inferred forces di↵er from the actual forces, with the average difference between inferred and actual forces being less than 1% of the mean force in the system. These di↵erences likely reflect errors in the optimization procedure due to machine precision.

Force Inference with Noise
Noise was introduced into each particle's stress tensor at steps 2 and 3 by adding a noise tensor, e i j , with each diagonal component drawn from a separate normal distribution with mean zero and standard deviation¯ zz /2 and each o↵-diagonal components drawn from a normal distribution with mean zero and standard deviation¯ zz /4 (preserving symmetry). The same procedure was repeated for twice these standard deviations. The use of normal distributions is supported by our previous work [10] (e.g., by visual inspection of Fig. 3 in that reference). Force inference was then performed using Eq. 3. Figure 2(a) shows the normalized error in inferred force magnitudes with noisy stresses. This error was calculated by first computing the absolute value of force magnitude error at each inter-particle contact and then dividing by the mean force in the system. Three aspects of the result are noteworthy. First, an increase in noise level relative to the mean stress in the system increases errors in inferred force magnitudes. Second, force inference in systems with increasing coordination number (step 3 versus 2) do not demonstrate an increase in force magnitude error. This indicates that force inference error does not increase in systems for which the number of unknowns approaches the number of equations. Finally, the most important observation from Fig. 2(a) is that the normalized error in inferred force magnitude is generally smaller than the mean force, and decreases as the ratio between the stress noise level and the sample stress decreases. This finding supports the notion that the statistics of the strong force network are insensitive to noisy stress measurements, particularly as the sample stress is increased relative to the noise level. Figure 2(b) shows the orientation error for all inferred forces shown in Fig. 2(a). The orientation error was calculated by finding the angle between inferred and actual forces. We find that orientation error increases only minimally as the noise level increases. Figure 3(b) and (c) show the inter-particle forces inferred in the presence of the two noise levels. Although noisy stresses cause errors in inferred forces, careful inspection of the figures reveals that the salient features of the force network, particularly the strong force network, are still captured well.

Conclusion
We have shown that particle stress tensor noise is fixed and is a function of the strain resolution of 3DXRD measurements (hardware and algorithms). We have further shown that the strong force network is accurately inferred using the force inference procedure, and that errors in force magnitude and orientation generally decrease as the sample stress increases relative to the stress tensor noise level. In our previous work, we have examined systems in which the peak sample stress is slightly greater than