Enlightening force chains: a review of photoelasticimetry in granular matter

TL;DR

This review explores the history, applications, and protocols of photoelasticimetry in studying granular materials, highlighting its scientific, educational, and public outreach uses over the past 50 years.

Contribution

It provides a comprehensive overview of photoelastic techniques in granular matter research and offers practical protocols for diverse users, fostering broader understanding and application.

Findings

Photoelasticimetry reveals internal stresses in granular materials.

The method has been influential in physics research and public education.

Protocols are provided for artists, educators, and scientists.

Abstract

A photoelastic material will reveal its internal stresses when observed through polarizing filters. This eye-catching property has enlightened our understanding of granular materials for over half a century, whether in the service of art, education, or scientific research. In this review article in honor of Robert Behringer, we highlight both his pioneering use of the method in physics research, and its reach into the public sphere through museum exhibits and outreach programs. We aim to provide clear protocols for artists, exhibit-designers, educators, and scientists to use in their own endeavors. It is our hope that this will build awareness about the ubiquitous presence of granular matter in our lives, enlighten its puzzling behavior, and promote conversations about its importance in environmental and industrial contexts. To aid in this endeavor, this paper also serves as a front door to a detailed wiki containing open, community-curated guidance on putting these methods into practice.

Figures

Equations

$$n_1-n_2=C({\sigma }_1-{\sigma }_2),$$

$$\alpha =\frac{2\pi Cd}{\lambda }({\sigma }_1-{\sigma }_2),$$

$$I=I_0{\sin }^2\frac{\alpha }{2}={\sin }^2\left[\frac{\pi Cd}{\lambda }({\sigma }_1-{\sigma }_2)\right].$$

$$\begin{array}[]{l}G^{2}\equiv\braket{}{\nabla I_{i,j}}{{}^{2}}=\frac{1}{4}\langle[(I_{i-1,j}-I_{i+1,j})^{2}/4\\ +(I_{i-1,j-1}-I_{i+1,j+1})^{2}/8+(I_{i,j-1}-I_{i,j+1})^{2}/4\\ +(I_{i+1,j-1}-I_{i-1,j+1})^{2}/8]\rangle\end{array}$$

Full text

††journal: Granular Matter

1 Introduction

Most of the transparent objects we encounter are photoelastic: their degree of birefringence depends on the local stress at each point in the material Frocht (1969); Cloud (1995). This property can be used to visualize, and even quantitatively measure, what is usually invisible to our naked eye: the stress field. When such materials are subjected to an external load, and placed between crossed polarizing filters, each different region of the material rotates the light polarization according to the amount of local stress Daniels et al. (2017). This creates a visual pattern of alternating colored fringes (see Fig. 1) within the material which, on top of their aesthetic and pedagogical aspects, permits us to quantify the stress field within the material.

Photoelastimetry has its roots in engineering practice, where it was widely used to design parts before the rise of computational finite element methods Cloud (1995). It also provided the first glimpse of the internal forces within granular materials, at first qualitatively Wakabayashi (1950); Dantu (1957); A and Josselin (1972); Liu et al. (1995) and later quantitatively Howell et al. (1999); Majmudar and Behringer (2005). Today, it remains the most well-developed method for quantifying stresses Amon et al. (2017), and methods for particle making Cox et al. (2016); Barés et al. (2017a) and image post-processing Daniels et al. (2017); Kollmer (2018); Lantsoght and Docquier (2018) are under active development.

After two decades of quantitative efforts, the scientific successes of the photoelastimetry method are numerous. In Robert Behringer’s group alone, it was responsible for identifying the erratic stress fluctuations in sheared granular matter Howell et al. (1999); Barés et al. (2017b); Zadeh et al. (2018), Green’s function response Geng et al. (2001), particle-scale anisotropy of the contact force networks Majmudar and Behringer (2005); Zhang et al. (2010), shear jamming Majmudar et al. (2007); Bi et al. (2011); Ren et al. (2013); Zheng et al. (2014); Wang et al. (2018), the dynamics of granular matter under impact Clark et al. (2012); Lim et al. (2017); Zheng et al. (2018), the Reynolds pressure, the Reynolds coefficient Ren et al. (2013), and more.

Far beyond the bounds of his laboratory at Duke University, the method has been used to examine particle shape dependence Zuriguel and Mullin (2008), identify interparticle contacts Lherminier et al. (2014), observe sound propagation Shukla (1991); Owens and Daniels (2011); Huillard et al. (2011), test the validity of statistical ensembles Puckett and Daniels (2013); Bililign et al. (2019), examine sensitivity to initial conditions Kollmer and Daniels (2018), identify dilatancy softening Coulais et al. (2014), measure force chain order parameters Iikawa et al. (2016), and observe the effects of fluid flow Mahabadi and Jang (2017). In interdisciplinary efforts, photoelastimetry permits scientists to evaluate the grain-scale stresses caused by growing plant roots Wendell et al. (2012); Kolb et al. (2012); Barés et al. (2017a), and examine situations relevant to faulting and earthquakes Daniels and Hayman (2008); Hayman et al. (2011); Geller et al. (2015); Lherminier et al. (2019)

The structure of the paper is as follows. First, we briefly review the physics of photoelasticity in §2; this section can be skipped for those only interested in qualitative uses of the method. In §3, we present various ways of fabricating photoelastic particles by cutting, casting or printing, followed by imaging-techniques in §4; these two section can stand alone for the creation of a demonstration apparatus. Finally, we present quantitative methods in §5. In all cases, additional information and technical specifications are provided on a wiki to which many of the paper authors have contributed Zadeh et al. (2019).

2 Photoelasticimetry theory

Photoelasticity arises from the birefringent properties of most transparent materials, in which the speed of light (via the index of refraction) depends on the polarization of the incident light wave. In some cases, such as glass and polymeric materials, birefringence arises only when the material is subject to anisotropic stress, with the refractive indexes depending on the eigenvalues of the local stress tensor. As such, photoelasticity can provide measurements of the internal stress in the material.

Measurements are best taken using circularly (rather than linearly) polarized light, in order to provide isotropic measurements. Circularly polarized light is created using a linear polarizer followed by a quarter-wave ($\pi/2$) phase shift between the two orthogonal components, as shown in Fig. 2. On the other side of the birefringent material, a circular polarizer with opposite polarity (the “analyzer”) blocks any light that doesn’t match its polarization. If the material is unstressed, there is no transmitted light and a dark image results. However, anyplace in the material where there is anisotropic stress, the wave components which are polarized along the two principle axes of the local stress tensor will travel with different speeds. This speed difference results in a relative phase shift for these two components of the wave, converting circularly polarized light to elliptically polarized light. As a consequence, a portion of the wave is not completely blocked by analyzing polarizer, and is therefore recorded as a bright region of the image.

This property allows photoelasticity to be used to quantitatively measure local stress, via an inverse method. We begin by assuming that the relation between the local stress and the refractive index is linear. We consider the difference in transmission between the two principle axes:

[TABLE]

where $\sigma_{1},\sigma_{2}$ are the two eigenvalues of the local stress tensor, and $n_{1},n_{2}$ are the two refractive indices in the corresponding directions. The material constant $C$ is known as the stress-optical coefficient. The relative phase shift of wave components in the eigendirections of the local stress tensor is determined from

[TABLE]

with $\lambda$ the wavelength of the incident light and $d$ the distance traveled inside the material (its thickness). For this phase shift, the intensity of the wave that emerges from the analyzer is given by

[TABLE]

More details about these relationships and how to calibrate the material parameters for quantitative measurements can be found in Trushant S Majmudar (2006); Daniels et al. (2017) and §5.

Note that while Eq. 3 relates internal stress and image intensity for a single wavelength of light, the effect is also preset for a superposition of wavelengths (e.g. white light), as shown in Fig. 1. The difficulty is that inverting Eq. 3, to infer local stress from light intensity, requires the use of a single wavelength in order to be tractable. Even so, the non-uniqueness of the solution due to the $\sin^{2}$ term makes the inversion problem challenging. §5 of this paper provides techniques for performing this task.

3 Fabricating photoelastic particles

While many transparent materials have photoelastic properties, only some of them have a large enough value of the stress-optical coefficient $C$ (see Eq. 1) to produce a measurable effect under reasonable loads. Furthermore, it is necessary to pay special attention during the fabrication of a material in order not to create an object containing significant residual stresses. Therefore, the creation of a photoelastic granular material requires considering all of the following: (1) the load you will to apply; (2) whether the shape of the particles matters; (3) how precisely you expect to make quantitative measurements, if at all; (4) what imaging method you will use; and (5) budget available. In this section, we explain several ways to meet these goals, ranging from simply cutting particles from pre-existing sheets up to casting bespoke particles.

Before choosing the material for the particles, it is important to consider that the photoelastic signal is a periodic function of the stress (see Eq. 3). The larger the deformation, the higher the stress, and the larger the number of fringes will be observed (see example in Fig. 2). Therefore, increased material stiffness is required for experiments with larger load, so that the fringes do not become denser than the resolution of the imaging system. Conversely, if the material is too stiff (or the loading too weak, or the thickness too small), then insufficient photoelastic signal will be observed. It is therefore advisable to perform some preliminary trials before committing to a large batch of particles.

3.1 Cutting sheets

Historically, it was simplest to make photoelastic particles by simply cutting shapes out of a pre-existing sheet of photoelastic material. This could be a flat sheet of Plexiglas*™* or a rubber, which is then cut with a bandsaw (the original Behringer particles), milling machine, spinning cookie cutter, or waterjet. The choice of material needs to match the experimenter’s dual requirements of deformability and sensitivity ($C$ in Eq. 1), to be compatible with the intended applied load. Most commonly, it is convenient to simply purchase sheets of transparent Vishay PhotoStress*™* vis or polyurethane pre , but for high loads Plexiglas*™* or polycarbonate are also suitable choices. The stiffer materials can be machined using any appropriate tool, with residual stresses annealed out by heating them to just below their glass transition temperature. In all cases, a diversity of thickness, stiffness, and color are available. Sheets of thickness $1/4$” and hardness $60$A are quite broadly-applicable, such that some particles of this type remain in use after two decades of use Howell et al. (1999); Barés et al. (2017b).

If circular-shaped particles (discs) are desired, a custom-built “cookie cutter” tool can make numerous, identical particles from a single sheet with little waste Daniels et al. (2017). Both the rotation and downward cutting speeds have to be properly chosen as a function of the material not to induce residual stresses when cutting. An example of particle obtained with this method in given in Fig. 3a, with slight horizontal marks made by the cutter. For more complicated shapes (see Fig. 3b), a computer-controlled milling machine outfitted with a narrow-diameter mill can trace arbitrary outlines, again with care taken to minimize residual stresses. In this case, there is also some roughness at the outer surface. Finally, it is possible to create arbitrary shapes with waterjet cutting, first used by Wendell et al. (2012) and further developed by Wang (2018). Depending on the skill of the operator, the resulting particles can have straight and smooth edges with few residual stresses, but can also sometimes leave a narrow channel at the start/end point of the cut. The narrow cutting width of the waterjet ($\sim 0.1$ mm) has the benefit of permitting very complex edges, as shown in Fig. 3c.

3.2 Casting particles

A second method to make photoelastic particles is to mold them directly, from such materials as polyurethane, gelatin, or nearly any other castable polymer or water gel. This allows an even larger diversity of complex, 3D shapes Barés et al. (2017a) so long as a mold can itself be fabricated. Using this method, it becomes possible to tune the stiffness (via the controlled addition of crosslinking molecules) or to add inclusions Cox et al. (2016) as shown in Fig. 3d.

The first step of the casting method requires creating a backing mold: a positive relief of the desired particle shape. The backing mold can be machined or 3D printed of nearly any material stiff enough to maintain a shape. As shown in Fig. 3f, this backing mold is then used to cast the final (reusable) silicone mold that makes the actual particles. A commercial silicone mold formulation such as MoldStar*™* provides easy-to-use formulations mol ; we have found that 15 Slow fits most needs.

Urethane is one popular choice of particle material Barés et al. (2017a). The commercial product ClearFlex*™* cle is available in different stiffnesses ($50$A fits most needs), and it can be custom-tuned by varying the crosslinker ratio. For the benefit of particle-tracking (see §5.1), it can be helpful to dye the clear urethane (see Fig. 5). The product SoStrong*™* so_ provides suitable dyes. Casting urethane requires some care to avoid the production of bubbles, compensate for material shrinkage, and develop fast enough work-flows; a number of helpful tricks from the community are shared at the online wiki Zadeh et al. (2019).

Another popular choice is to use biological gels which are cheap and easy to use, but non-permanent. As shown in Fig. 3e, gelatin has excellent photoelastic properties Kilcast et al. (1984); Lim et al. (2017); Workamp et al. (2017), as does agar or konjac Tomlinson and Taylor (2015) (with the later being less transparent). The stiffness of these materials is easy to tune by varying the ratio of gelling agent, and it is possible to achieve arbitrarily low elastic moduli. However, because biological gels are composed of water and a food source, they are vulnerable to drying, swelling, and bacteria; so they are not stable over long times. Some groups have had success by crosslinking the gel Damink et al. (1995); Workamp et al. (2016) when making particles Workamp et al. (2016). In this case, glutaraldehyde is directly added to the liquid gelatin preparation before molding or it can be diffused into the gelatin particles once they have gelled Workamp et al. (2016) to increase the particle stiffness and stabilize the material.

3.3 3D printing particles

Finally, when the particle shape is sufficiently complex that a backing mold cannot be made, it is possible to simply 3D print individual photoelastic particles. This choice is not as straightforward as it might seem, since most of the available printing materials are non-transparent, porous, and contain residual stresses. However, optical fabrication processes (stereolytography) now allow for the production of semi-transparent photoelastic particles. Two available options are Durus White*™* used by Lherminier et al. (2016), and VeroClear*™* ver used by Wang et al. (2017). The VeroClear has better clarity and optical response, but both are quite stiff and are therefore only appropriate for experiments with large loading stresses.

4 Imaging methods

In designing a photoelastic experiment, it is also necessary to consider the imaging method. The most appropriate choice will depend on considering the magnitude of the photoelastic response (§2) for the applied load, as well as desired degree of precision and the speed with which the images need to be collected to capture the dynamics. In this section, we elucidate several different methods that allow for adapted to a variety of constraints.

4.1 Transmission vs. reflection imaging

In §4 and Fig. 2, we considered a geometry in which the light is transmitted directly through the particles from a polarized light source, and then through a second filter (analyzer) before reaching the camera. In cases where the experiment will be optically-accessible from both sides of the particles, this is the easiest polariscope to construct, and has seen the most usage over the past several decades Liu et al. (1995); Majmudar and Behringer (2005); Daniels et al. (2017). Care must be taken when constructing the optical pathway of filters, that the polarizers be perfectly crossed and aligned exactly as shown in Fig. 2. This can be done either by purchasing two pairs of linear polarizers and quarter-wave plates pol and doing the alignment by hand, or by purchasing a pair of left and right circular polarizers which are pre-aligned sandwhiches containing both components. In that case, care must be taken to orient the quarter-wave plate side of the polarizer towards the granular material so that light passes through the filters in the correct sequence. If the circular polarizers are pre-mounted inside a camera filter, their default configuration will be backwards from what is required for the polariscope shown in Fig. 2.

In some cases, the granular material may not be optically-accessible from both sides, for instance due to the loading mechanism or because the particles are resting on an opaque surface Puckett and Daniels (2013); Zhao et al. (2017). It can therefore be desirable to create an optical setup in which the the apparatus is lit and imaged from the same side: this is a reflective polariscope. As in transmission polariscope, both the polarizer and the analyzer are circular, but now a single polarizer serves in both roles (see Fig. 4). Two successful options for reflecting the polarized light back through the granular sample are to rest the particles on a mirrored surface, or to coat all particles with a mirror-effect paint. Details about how to construct such an apparatus are provided in Daniels et al. (2017).

4.2 Multi-wavelength imaging

As described by Eq. 3, quantitative stress measurements require monochromatic light measurements. In order to minimize the overlapping of photoelastic fringes at high stresses (as can seen in Fig. 2), it is important to work with monochromatic light. Furthermore, polarizers and quarter-wave plates are optimized for a given wavelength, usually green light. Therefore, a quantitative apparatus should be designed so that green light is used for photoelastic measurements, and other wavelengths are used for monitoring quantities such as particle positions and orientations.

One method, suitable for quasi-static dynamics, is to perform sequential imaging after each loading step Ren et al. (2013); Barés et al. (2017b); Zhao et al. (2017). For example, the camera takes an image (1) between crossed polarizers, to observe the photoelastic response; (2) without at least one of the polarizers, to detect particle positions; and (3) additionally illuminated by UV light, to monitor an additional characteristic such as particle orientation. An example of this type of imaging is shown in Fig. 1 and Fig. 6d, in which a bar drawn in UV-sensitive ink is visible.

Another option is to measure these same fields simultaneously, by choosing to illuminate different properties using different wavelengths of light. Conveniently, color cameras already perform color-separation into red, green, and blue (RGB) channels. For example, the particle positions can be illuminated in unpolarized red light, while the photoelastic measurements are recorded by polarized green light, and a subset of tagged particles identified by UV-ink tags that emit blue light Puckett and Daniels (2013); Kollmer and Daniels (2018). Such a setup is shown schematically in Fig. 4b, with a sample image in Fig. 4c.

A third approach is to use dyed particles, which can then be illuminated with white light. In this case, the particles effectively act as a bandpass filter, and the passed wavelength can be used to make photoelastic measurements while blocked wavelengths provide position data. A particle dyed dark blue is shown to demonstrate this effect in Fig. 5. These techniques allow a single image to provide both position detection and stress measurements at the same instant in time, making it possible to monitor the dynamics of a system.

5 Analyzing images

We have seen numerous examples of photoelastic images taken from a variety of geometries and lighting conditions; next, we examine several key methods of extracting quantitative information from such images. The level of precision available to the user – whether particle position, orientation, or force – depends on the focus and resolution of the images, the photoelastic sensitivity of the particles (parameter $C$ in Eq. 1), and the magnitude of the applied load. In what follows, we will focus on the use of cylindrical (disk) particles.

5.1 Tracking particle positions and orientations

For particle detection, the first step is to produce a high-contrast (even binary) image to distinguish the pixels occupied by particles from those that are not. This is best done on whichever color-channel (red, green, or blue) has the highest contrast between the particles and their background. It can be helpful to increase the contrast by filtering single high/low outlier pixel values (see Fig. 6a), and also to apply a low-pass Fourier filter (see Fig. 6). These steps work for particles of any shape.

Starting from a filtered image improves particle-finding methods, of which we will consider the two most convenient: convolutions Shattuck (2015) and Hough transforms. If the particles are not circular, then the convolution method is preferred over the Hough transform; it can also achieve higher precision but is more computationally-intensive. Image convolutions are performed between a pre-set image (“kernel”) of a single ideal particle (approximately $3/4$ the diameter of the smallest particle works well), and the original or binarized image. After performing the convolutions, binarizing the resulting image using a threshold of $99.5\%$ of the peak convolution value will result in a field of well-isolated peaks corresponding to the particle centers, with the area of each peak indicating that particle’s diameter. For circular (non-binarized) images, it is possible to use a circular Hough transform to locate the particle centers and sizes Peng et al. (2007); Puckett and Daniels (2013). This transform works best on an image that has been prepared, via edge-detection, to isolate the circumference of particle. The algorithm uses a computationally-efficient voting process to identify the centers of such circles of a specified radius. A sample result of particle-detection via the Hough transform in shown in Fig. 6c.

From the the list of particle centers present in each image of a series, it is possible to track Lagrangian trajectories of those particles through space and time. Two efficient, accurate, and open-source algorithms are those of Crocker and Grier (1996) and Blair and Dufresne . The main idea for both is that, given the positions for all detected particles at a given step, the aim is to select the closest possible new positions from the list of all detected particles in the next step. To this end, the algorithm considers all possible pairings of (old-new) positions, and selects the pairing that would result in the minimal total squared displacement. If a particle is missed over a small interval, the two halves of the trajectory later be re-attached. In cases where the image resolution was insufficient to reliably detect particle positions, it is still possible to use particle image velocimetry (PIV) to measure Eulerian flow fields for a series of images Willert and Gharib (1991).

Finally, it is possible to tag the particles with a contrasting stripe, in order to measure particle-rotation during dynamics. One common method, shown in Fig. 6d-f, is to mark the top of each particle with a UV ink bar. This ink is not visible under white light, and therefore doesn’t interfere with ordinary particle-tracking as described in the first part of this section. However, once illuminated by UV light alone, the ink bars become visible with high enough contrast to be recorded by a camera. The contrast is sufficient to binarize the image, and the position of each particle center can be used to conduct a least squares of the bar pixels in order to determine its orientation (see Fig. 6f).

5.2 Scalar measurements of force

For many experiments, it is desirable to make quick, low-resolution measurements of the forces on each particle, rather than solving the full inverse problem to obtain vector contact forces as will be described in §5.3. Scalar measurements can be the best choice for a variety of reasons: poor image resolution, dim lighting, fast dynamics, or a very large number of images. It is the easiest analysis when beginning a new project, as a tool for preliminary investigations, and remains popular because if its computational efficiency.

The simplest form of scalar analysis is to use the mean light intensity within particles. For small forces, the mean image intensity increases linearly with the applied force, as shown in Fig. 7a for a sample experiment under a known imposed load. By analyzing a small system such as a linear chain, it is possible to perform a calibration which provides a quantitative measure of the force on particle. As shown in Fig. 7b for a numerical version of this calibration, there is a threshold force above which the mean light intensity plateaus and this method is no longer quantitative. This threshold corresponds approximately to the development of the first additional set of fringes. For applied forces within the linear regime, this method is valid.

Once this threshold is crossed, the variation in intensity due to the fringes becomes a significant factor, allowing for an additional quantitative technique which accounts for the presence of the fringes. For images at high enough resolution to observe fringes, the force on the particle can be measured using the gradient of the image intensity field $I(i,j)$, known commonly as the $G^{2}$ technique. It is defined as follows:

[TABLE]

where $i$ and $j$ refers to the row and column number of pixels and $\braket{\cdot}$ means averaging over pixels in the region of interest Majmudar (2006). Figs. 8 shows the distribution of light intensity $I$ (a) and its squared gradient $|\nabla I|^{2}$ (b) for a photoelastic disk under normal compression.

This method has been used to measure both the average pressure for the whole image Howell et al. (1999); Geng et al. (2001); Farhadi and Behringer (2014); Zheng et al. (2014); Farhadi et al. (2015); Cox et al. (2016); Daniels et al. (2017); Zheng et al. (2018), and at the particle scale Clark et al. (2012); Coulais et al. (2014); Iikawa et al. (2015); Barés et al. (2017b); Wang et al. (2018). For bulk measurements, we observe that $G^{2}$ is a monotonic function of pressure, with some nonlinearity in the relationship Howell et al. (1999). For particle-scale measurements on disks, if applied forces are small enough such that photoelastic fringes can be clearly resolved, $G^{2}$ is proportional to the sum the individual vector contact forces $\sum_{i}|\vec{F}_{i}|$, where ${\vec{F}}_{i}$ are the contact forces on the disk Zhao et al. (2019) (see Fig. 8a). Thus, where tangential forces are not large, $G^{2}$ approximately measures the particle scale pressure.

Importantly, $G^{2}$ will differ based on the material choice, the lighting condition and the image resolution. At particle scale, for an individual disk, the coefficient of proportionality $k=G^{2}/\sum_{i}|\vec{F_{i}}|$ is shown to be proportional to $I_{0}^{2}Cd/\lambda R^{2}N$ Zhao et al. (2019), with $\lambda$ the light wavelength, $C$ the stress-optical coefficient, $d$ the disk height, $R$ its radius, $N$ the camera resolution defined as number of pixels per meter and $I_{0}$ the background light intensity. As an example, Fig. 8c plots the collapse of the linear part of $G^{2}$ rescaled using $N$ and $R$ as functions of $\sum_{i}|\bf{F_{i}}|$ for diametrically loaded disks (see Fig. 8a) with other parameters kept the same. The dependence of $G^{2}$ on $R$ is important for analyzing poly-disperse systems. Fig. 8c also compares the dependence of $G^{2}$ and intensity $I$ on the contact forces for a given test, showing that the range of force that $G^{2}$ can measure is $4$ to $5$ times larger than that of $I$. Note that the linear dependence of $G^{2}$ value for individual disk on contact forces will differ for other particle shapes Zhao et al. (2019).

5.3 Vector measurements of force

The detailed pattern of light and dark fringes in $I(i,j)$ is set by the local stress values through Eq. 3. Therefore, it is possible to use a known stress field (computed from the set of vector contact forces ${\vec{F}}_{i}$) to predict a light intensity field. Examples of such calculations are shown in Fig. 7b. In order to determine the values for each ${\vec{F}}_{i}$ from the fringe pattern requires performing the inverse of this process, as illustrated in Fig. 9. Because Eq. 3 cannot be directly inverted due to the $\sin^{2}$ non-monotonic function, the numerical estimation of the set of ${\vec{F}}_{i}$ must be done via an optimization process. This is achieved in a sequential process:

Detect the particle positions (see §5.1) 2. 2.

Estimate the total and individual contact force magnitudes on each particle (see §5.2) 3. 3.

Optimize ${\vec{F}}_{i}$ on each particle according to a set of constraints

Note that the second step contains two parts. It is usually possible to use the $G^{2}$ method to estimate the total force, but this may not be possible for the individual contacts, depending on the image resolution. If the quality of the photoelastic image is good enough, each contact force can be similarly estimated from calculating $G^{2}$ in the vicinity of each contact Daniels et al. (2017); Kollmer (2018). This method allows for the removal of contacts that do not transmit loads, and for the inverse problem to be performed on individual particles. However, the contact-scale step is not always possible for lower-resolution images, for example in the case of stirred granular media Lantsoght (2019). Instead, the optimization step (below) must also include a process for testing multiple different contact-configurations in which the total force magnitude is taken to be equally-distributed in the initiation of the optimizer Lantsoght and Docquier (2018).

The third step proceeds to optimize the vector contact forces ${\vec{F}}_{i}$, starting from the estimated values in steps 1 and 2, in order to achieve a photoelastic response as close as possible to the image intensity $I(i,j)$. The input (initial guess) and output (optimized) parameters are therefore the set of all contact force magnitudes and orientations. For cylindrical, linearly-elastic particles subjected to multiple contact forces, the stress-field can be approximated by an analytical expression Timoshenko and Goodier (1970); Daniels et al. (2017).

Combining this stress field with Eq. 3 creates a reconstructed image of each particle which can be evaluated against the measured image. The chosen optimization algorithm works by evaluating the agreement between these two images – the measured $I(i,j)$ and the numerically-generated $I_{n}(i,j)$ – and modifying the values of ${\vec{F}}_{i}$ to bring them into closer agreement. The similarity between $I,I_{n}$ can be computed with such quantities as the mean squared error which is quick to calculate, or the structural similarity index Wang et al. (2004) which better-accounts for the the image structure. One commonly-used optimization protocol is Levenberg-Marquardt optimization Daniels et al. (2017).

No matter the choice of similarity index or optimization protocol, it is beneficial to have the initial guesses for forces to be of highest quality possible in order to have the optimizer convergence to a valid result. When the initial guesses are too far from the correct values, the optimizer may land in a local minimum and not escape. Ideally, a single run of optimization is sufficient to determine ${\vec{F}}_{i}$ for all particles, but in practice it can be beneficial to run the process multiple times. When the optimizer is having difficulty converging, for instance due to poor-quality initial guesses for the contact forces, it is possible to propagate already-identified values of ${\vec{F}}_{i}$ to their reciprocal force on the adjacent particle. This improves the initial guesses for the next particle, and this process can be repeated multiple times, sequentially propagating information through the packing. At the end of the optimization process, all contact forces (magnitude and orientation) have been determined for each of the particles in the system. At this stage, it is still possible to iteratively make additional improvements, for instance by checking for consistency with Newton’s Third Law, or the equilibrium condition on each particle Daniels et al. (2017).

6 Outlook

We close with a summary of the newest developments in photoelastimetry: taking these techniques and applying them to faster dynamics, non-circular particles, and three-dimensional systems.

Photoelastic measurements in fully three-dimensional systems is challenging, since a curved particle simultaneously acts as a lens. Nonetheless, it is possible to obtain semi-quantitative information about the stress state of the system using birefringent spheres Liu et al. (1995); Yu et al. (2014). A promising route for 3D studies is therefore to use terahertz photoelasticity, where the wavelength of light is hundreds of microns. Although early measurements confirmed the translucence of photoelastic materials at terahertz frequencies Mahon et al. (2006); Heimbeck et al. (2011), the first attempts at measuring strain birefringence failed because the microscopic displacements of atoms were too small compared to the wavelength of the radiation. However, it is possible to develop metamaterials whose meta-“atoms” exhibit greater displacements. An example of such a material is shown in Fig. 10ab. The basic principle, recently demonstrated Everitt et al. (2019), is that an applied stress reversibly distorts the shape of the meta-atoms, and this anisotropy is detected when the material is placed between crossed polarizers and illuminated from above by a terahertz generator. To make truly 3D measurements, this system can be combined with terahertz holographic imaging Mahon et al. (2006); Heimbeck et al. (2011).

Another challenge in applying the techniques described in §5 to more realistic situations, is to allow their application to particles that are non-circular or cohesive, such as commonly arise in geophysical and industrial contexts. For circular particles with cohesion (i.e. electrostatic interactions, liquid bridges, chemical bonds), the inverse methods described in §5.3 likely still apply. However, the key technical challenge for non-circular particles is that the inverse problem is only well-specified for the case of circular particles, where the trivial surface geometry (all tangential forces are torques, all normal are central) dramatically simplify the formulation of the inverse solution to Eq. (3). New mathematical solutions are required to created a more general formulation to apply to particles of arbitrary shape.

Another persistent challenge is to move from quasi-static experiments to fully dynamic. Previously, this was accomplished by performing semi-quantitative force analysis, as described by the $G^{2}$ method in §5.2. Recently, advances in high-speed imaging have made it possible to record images with sufficient brightness and resolution to perform photoelastic analysis on continuously-avalanching flows, and thereby obtain vector contact forces via adaptations of the methods in §5.3 Thomas and Vriend (2019). This process is enabled by measuring the particle positions directly from a single movie of photoelastic images, in order to capture as much light as possible. Furthermore, because the particles are accelerating, it is not possible to use particle-scale force-balance as a constraint in determining the vector contact forces.

Finally, computational advances in inverting the photoelastic images would improve the applicability of photoelastimetry to all of the problems presented here. Better algorithms for avoiding unintended-minima and the ability to simultaneously solve for all vector contact forces in the system (a massive set of constraints) would improve the reliability of these methods.

It is our hope that the recent proliferation of publications of photoelastimetric techniques in open formats Zadeh et al. (2019); Lantsoght and Docquier (2018); Kollmer (2018) will spur development in all of these new directions.

Acknowledgements

We would like to thank Rémy Mozul for his technical support with the wiki Zadeh et al. (2019). Several conversations and collaborations have led to sharing the techniques described in this paper. We are grateful to Bernie Jelinek and Richard Nappi for sharing their technical knowledge about photoelastic material cutting. The outlook section contains insights gained from Chris M. Bingham, Willie J. Padilla, Anthony Llopis and Nan M. Jokerst (recent work on terahertz photoelasticity), and from Nathalie Vriend and Amalia Thomas (fast-imaging photoelasticity).

Finally, we thank the late Robert Behringer for his kindness, his depth of knowledge gained from developing photoelastic techniques for two decades, and his stimulating attitude towards every new generation of scientists passing through his laboratory. This review article is a product of his excellent mentorship.