Photoelastic force measurements in granular materials

Photoelasticity has long been used by engineers to quantify the internal stresses within solid bodies.¹ The technique is based on the understanding, dating to Maxwell,² that light can be polarized (the electric and magnetic fields have a well-defined orientation) and its speed depends on the medium’s index of refraction. Many crystalline materials are birefringent: their crystalline axes specify a fast and slow direction for the propagation of light.

Non-crystalline materials such as polymers and glasses can nonetheless be photoelastic, with the degree of birefringence at each point in the material depending on the local stress. When a photoelastic material is placed between two polarizers and subjected to stress, each region of the material rotates the polarization of light according to the amount of local stress and the stress-optic coefficient of the material. This creates a visual pattern of alternating bright and dark fringes within the material, and this pattern depends in a detailed way on the orientation of the polarizers, the shape of the material, and how it is stressed.

Before the creation of finite-element computing methods, photoelastic analysis was a key technique for predicting regions of high compressive or tensile stress during the design of manufactured parts. As such, it has been an important part of engineering practice in both its qualitative and quantitative forms.^3–6 It was therefore inevitable that photoelastic methods eventually be applied to the study of idealized conceptions of granular materials, as was first done by Wakabayashi using epoxy tracer particles in 1950s,⁷ followed by Dantu (using ordered/disordered arrays of glass cylinders⁸) and later by Drescher and de Josselin de Jong⁹ (in polymer disks) and Liu et al.¹⁰ (in glass spheres). The network-like structures visible in Figs. 1 and 2 have come to be known as force chains.

These pioneering experiments played a key role in the physics community’s realization that heterogeneous force transmission identifiable at the particle scale is of great importance for understanding the mechanics of granular materials.^10,11 A series of papers by Behringer and co-workers established the photoelastic method as the key technique for making quantitative advances.^12–14 In the past two decades, many different experimental teams have used photoelastic force visualization to quantify myriad phenomena. For example, it has been possible to reveal the particle-scale anisotropy of the contact force networks,¹⁵ identify interparticle contacts,^16,17 measure fabric and stress tensors,¹⁸ examine particle shape dependence,¹⁹ test the validity of statistical ensembles,²⁰ identify dilatancy-softening,²¹ measure force chain order parameters,²² evaluate the grain-scale stresses caused by growing plant roots,^23,24 observe the effects of fluid flow,²⁵ follow slow dynamics under shear,^13,26–28 and quantify fast dynamics.^29–32 Because the photoelastic response is known analytically for circles/ellipses, several of these studies have measured the vector contact forces at every interparticle contact within the granular material.^33–35 

This paper first provides a review of photoelasticity and how it can be used to create high-quality images of stress within the granular material (Section II). These images can be used either semi-quantitatively (Section III) or quantitatively (Sections IV and V), with the latter providing a method to determine the vector contact forces for each of many circular particles in the system. We provide a description of software resources available to perform this task. Finally, we describe key techniques for building an experimental apparatus (Section VI), including resources for the purchase of key components.

This section provides a brief overview of polarized light and photoelasticity, to aid the reader in building a qualitative understanding of the method. A standard optics textbook (such as Hecht³⁶) can provide quantitative details about polarized light, and the classic two-volume monograph by Frocht¹ provides a comprehensive treatment of photoelasticity as a method. Several recent book chapters also provide an introduction to the method’s use in the context of granular materials.^35,37

Electromagnetic waves travel in a direction co-perpendicular to the sinusoidal electric and magnetic fields of which they are comprised (see Fig. 2(a)). Linearly polarized light consists of waves that all have the same orientation of their electric fields. Two common ways to generate linearly polarized light are reflection from a metallic surface and transmission through a polarizing filter (also known by a brand name, Polaroid^TM). These filters consist of a thin sheet of iodine-impregnated polyvinyl alcohol, with polymer chains preferentially aligned along one axis. The mechanism through which they polarize the light is that the iodine provides conduction electrons which polarize along the polymer chains and thereby cancel the electric field of the incident light polarized in the parallel direction but not in the perpendicular direction. Typically efficiencies are to absorb >95% of light in the parallel direction and transmit <40% of light in the perpendicular direction.³⁸ 

Circularly polarized light consists of electric and magnetic fields that remain mutually perpendicular but rotate with constant magnitude with left-handed (counter-clockwise) or right-handed chirality (clockwise) with respect to the direction of propagation. To create circularly polarized light, light is first linearly polarized and then passes through a quarter wave plate that advances the phase of one polarization with respect to the other. A quarter wave plate is a birefringent material of chosen orientation (fast/slow axes) and thickness matched to the wavelength of interest. The orientation of the fast axis at $± 45 °$ (see Fig. 2) with respect to the linear polarization determines whether the output is left-circularly polarized (LCP) or right-circularly polarized (RCP). Typically, the two optical elements are purchased as a single, fused unit. Note that, unlike linear polarizers, a circular polarizer can be arbitrarily rotated by any angle around the transmission axis with no effect on the filtered light but cannot be flipped over.

The basic principle of a darkfield polariscope is to place a photoelastic sample between a polarizer and an analyzer of the opposite polarization (e.g., vertical and horizontal linear polarizers or LCP and RCP circular polarizers), viewed in transmission. Sample configurations are shown in Fig. 4. Note that care must be taken during the installation of circular polarizers, to ensure that the polarizer and analyzer are both placed with their quarter-wave plates facing the sample.

In each case, initially unpolarized light becomes polarized after passing through the polarizer. As the light passes through the sample, the amount of birefringence at each point depends on the local stress state (photoelasticity). This birefringence causes the light to rotate its polarization as it passes through the sample, and when it exits the sample the polarization has been rotated by a total of θ degrees. The light now passes through a final polarizer called the analyzer. The final intensity of light depends on the local stress in the birefringent sample as well as the orientation of the analyzer. In a stressed sample, this leads to an image of light and dark fringes (see Fig. 3).

While both linear and circular polariscopes are possible, there is a practical difficulty with linear polariscopes that limits their utility: In a linear polariscope, there are two possible reasons for fringes to appear: (1) the photoelastic response of the sample rotates the polarization to match/mismatch the analyzer; and (2) if the principal stress of the sample happens to be aligned with the polarizers, it will not be affected by the photoelastic response of the material. Therefore, for linear polariscopes, the fringe pattern depends on both the relative angles between the polarizer (fixed) and the principal stress (changes at different locations and times during experiment). For this reason, quantitative work with polariscopes should always be conducted with circularly polarized light, which provides azimuthal symmetry.

While the resulting image of light and dark fringes depends uniquely on the stress in the sample, there is an ambiguity in the phase. Various groups have proposed various methods to unwrap the phase map of the resulting image to extract the stress in the sample. However, many images of the sample are needed with different wavelengths of light³⁹ and different analyzer orientations.⁴⁰ These methods must resolve both isochromatics (associated with the magnitude of local stress) and the isoclinics (the angle of stress). In investigations of granular materials, we are solely interested in the contact forces of each sample rather than the local stress state itself. Therefore, it is sufficient for contact-force measurements to use an analyzer which has been arbitrarily rotated about its azimuthal axis (with respect to the polarizer); this results in an image of the isochromatics only. Throughout the remainder of the text, we assume images are acquired using the circular polarizers set up as shown in Fig. 4.

Generally speaking, the samples in Figs. 1 and 2 illustrate that particles with more force are brighter (darker) against a dark (bright) background because the light that passed through them had its polarization rotated. Often, this contrasting effect is sufficient to identify general features of stress-transmission within a granular sample. Methods for quantifying these features are described in Section III. Importantly, there is a quantitative, one-to-one mapping between the stress field in the sample and the resulting pattern of fringes.¹ Therefore, it is possible to solve the inverse problem and use polariscope images to determine the vector contact forces on each particle in the sample. Quantitative methods for performing this inversion are detailed in Sections IV and V.

Originally, images of the type shown in Figs. 1 and 2 were used for semi-quantitative investigations, without the knowledge of the vector contact forces between the particles. Even in modern experiments, there is much that can be learned from comparing the spatial patterns of the force chains. For example, images of the force chains can reveal history dependence via spatial correlations,^41,42 the ensemble-averaged response to point forces,⁴³ and the principal axis of the response.⁴⁴ By subtracting images taken before/after an event, it is possible to examine spatial patterns of failure^26,45 and the speed of sound.³¹ 

To quantify how the number of dark/light fringes increases with the stress on each particle, Behringer and co-workers^11,13,14,46 calculated the squared local gradient, averaged over all pixels in an image

They determined (see Fig. 5) that this quantity provides an empirical measure of the 2D stress on the system. More recently, this has come to be known as the gradient-squared (G²) method and has been used on the system-scale, chain-scale, particle-scale, and contact-scale to provide a semi-quantitative way to track changes in granular forces.^32,47–49

The G² method has several important advantages over the photoelastic inverse methods to be described in Section V. First, it works on low-resolution images (10 pixels/particle is enough) as long as there is sufficient intensity contrast. Second, image processing times are very fast: at the time of writing, it takes 0.05 s to calculate the pressure within a 1 megapixel image on a single processor. Finally, there are no parameters in Eq. (1). To translate G² into pressure, only a single calibration experiment needs to be done under the same lighting conditions as the experiment, in order to properly account for the image intensity and contrast. For these reasons, the G² method has remained popular, even as the speed and accuracy of vector contact force measurement codes have improved.

As illustrated in Fig. 6, a set of z contact forces on a disk determines the pattern of photoelastic fringes. Because the analytical solutions for the stress created by point contact forces on a disk are known,¹ it is possible to determine the intensity field $I ( x , y )$ from the set of z vector forces. In analyzing images from experiments (to be shown in Section V), we will need to do the inverse: for a known $I ( x , y )$ recorded by a camera, determine the vector contact forces which created that image. In this section, we present the theory required to determine the mapping from $f → → I$ and calibrate the measurements so that they can be used for the inverse problem in Sec. V.

The special case of diametric loading provides both an instructive example for understanding quantitative photoelasticity and the means of calibration necessary for any quantitative experiment. As illustrated in Fig. 7, the number of dark/light fringes increases with the diametric load of magnitude f (equal and opposite forces along the diameter of the disk). We can understand this behavior starting from the known solutions for the stress tensor at the center of the disk,^1,50

The intensity of the fringe pattern at any point is given by

where $( σ 1 − σ 2 )$ is the difference in the principal stress at the point (x, y), h is the material thickness, λ is the wavelength of light, and C is the stress-optic coefficient (a λ-dependent material property).

The spatial variations in I, visible in Fig. 8, are present because $σ 1 − σ 2$ is spatially varying. We will consider the simpler case of monochromatic light, so that $F σ = λ / C h$ (the stress-optic coefficient) is a constant. In the case of isochromatic fringes, Eq. (5) reduces to

Using Eq. (4), we have $σ y y = − 3 σ x x$⁠, where the off-diagonal components $σ x y$ and $σ y x$ are zero and the shear stress is zero at the center of the disc. The difference in the principal stress is the difference in the eigenvalues, which then is given by $σ 1 − σ 2 = | σ x x − σ y y |$⁠. Then, at the center of the disc, the number of fringes N_fr is given by

For an unstressed disk, $σ 1 = σ 2 = 0$ everywhere in the disk, the difference in the phase is zero and thus the isochromatic image is dark. The first bright fringe occurs when the difference in the principal stress is $σ 1 − σ 2 = F σ / 2$⁠. The bright and dark fringes on the disk arise from the periodic function $( sin 2 x )$ in the phase retardation. As shown in Fig. 7, the total number of fringes increases as the sin function cycles through an increasing number of bright and dark fringes for increasing values of $σ 1 − σ 2$⁠.

Thus, $F σ$ can be empirically determined by diametrically compressing the disk by a known force and recording the fringe number as follows. For diametric loading, the fringe number is defined as the number of fringes counted between the edge of the disk and the center of the disk, as shown in Fig. 8(a). The fringe number at the boundary is zero. Moving toward the center, each bright/dark band increments the fringe number by 0.5. In order to measure C (the material parameter), a typical experiment proceeds by placing a disk of known radius and thickness into a load cell and illuminating it with monochromatic, polarized light. As the diametric load on the particle is slowly increased, a measurement is taken each time a new fringe appears (N_fr increases). A scatter plot of the diametric force and number of fringes allows for a fit to the left and right sides of Eq. (7). Since all other parameters are known, this provides a measurement of C for that particular material, as shown in Fig. 8(b). Most of the commonly used materials (see Section VI A) have $C ≈ 3000$ Br.

Ultimately, we will determine all contact forces acting on a particle based on the pattern of isochromatic fringes recorded in images of photoelastic disks. This goal is illustrated for a single particle in Fig. 6. To achieve this will require a formula for the solution I(x,y), as a function of an arbitrary set of z vector contact forces, acting at arbitrary contact locations around the perimeter of a disk. In this section, we derive the solution to the stress field $σ ( x , y )$ inside an elastic disk in mechanical equilibrium with z point forces on the boundary, along with a short review of the necessary theory of elasticity. The solution presented here largely follows the work by Refs. 1, 50, and 51.

For any linear elastic solid with boundary forces, there are three basic assumptions: (a) the condition of equilibrium (force and torque balance); (b) a constitutive relation (here, linear stress-strain); (c) Saint-Venant’s condition of compatibility (no gaps or overlaps). In two dimensions and in the absence of body forces, the condition of equilibrium for a static solid is given by

Equivalently, $∇ ⋅ σ = 0$⁠. The second assumption is a linear relation between the strain and stress tensors,

The third assumption is Saint-Venant’s condition: the solid deforms continuously leaving no gaps and creating no overlaps.⁵² For plane strain in the absence of body forces, the compatibility of infinitesimal strains is given by

The solution to Eq. (10) in two dimensions is solved using Airy stress functions, φ, where the components of the stress tensor are derivatives of φ,

Then the compatibility equation of the stress function is biharmonic

In general, the solution to Eq. (12) is guessed by using symmetries or adjusting a known solution to the boundary conditions.

Consider a stress on a disk, as shown in Fig. 9, with radius R and thickness h. We define the origin at point O₁ where the direction of pressure is O₁O₂, with angles $θ 1$ and r₁ as shown in Fig. 9. Points on the circumference of the circle have a radial stress $σ r r = − 2 f π h cos θ$ in the direction r₁. The angle between $σ r r ′$ and the tangent is $θ 2$⁠. The components of the stress tensor are

Substituting $σ r r$ and using $θ = θ 1$ and $r 1 = 2 R sin θ 2$ everywhere on the circumference, then

These two tensors are the superposition of three stresses at point A on the circumference as follows:

• normal stress: $− f π R h sin ( θ 1 + θ 2 )$⁠.

• tangential stress: $− f π R h cos ( θ 1 + θ 2 )$⁠.

• uniform stress in the direction of r₁ with normal and tangential components: $− f π R h sin ( θ 2 − θ 1 )$ and $− f π R h cos ( θ 2 − θ 1 )$⁠, respectively.

If the disk is force and torque balanced, torque balance implies $∑ i z f i cos ( θ i , 1 + θ i , 2 ) = 0$⁠, and the sum of the uniform tension must also be zero. Therefore, the remaining term is the sum of the normal tensions at the boundary $∑ i z f i π R h sin ( θ i , 1 + θ i , 2 )$ along the circumference, which is simplified using the inscribed angle theorem $θ 1 + θ 2 = π 2 ± α$ and the identity $sin ( π 2 ± α ) = cos ( ± α ) = cos α$⁠. To ensure that the boundary of the disk is stress free, a uniform tension of magnitude $f π R h cos α$ is added to the radial stress distribution. The solution for z concentrated forces on a disk where $θ i$ is the angle between fi and r_i and α as shown in Fig. 9 is

All other components of σ are zero.

The solution is in z-polar coordinates. To write the solution in a single coordinate system, we use a Cartesian coordinate system with the origin at the center of the disk. Let β be the angle between the positive Y-axis and the line segment $C O 1 ¯$⁠, as shown in Fig. 9. We compute $σ r r$ as in Eq. (15) and transforming coordinate systems with a rotation by $θ i = − ( β i + α i )$ as

where $T ( θ i )$ is the rotation matrix.

The eigenvalues of σ are

where the subscripts 1 and 2 correspond to the + and − signs, respectively.

We have presented the solution to z-forces on a disk and given a calibration of $F σ$ for the photoelastic disk. Section V details the numerical methods involved in taking this solution and finding the contact forces in an array of disks.

The ultimate goal of the force-inversion process is to start from images such as Fig. 10, find all of the particles, and determine the vector contact forces between them using the theoretical framework described in Section IV. While the equations for the image intensity due to a known set of particle positions and contact forces are known analytically, the inverse of the problem is not. Several research groups^15,20,28,35 have created computational solutions to this inverse problem, including in C++⁵³ and MATLAB,⁵⁴ a technique pioneered more than a decade ago³³ by Majmudar.

Here, we describe our most recent implementation, PEGS, inspired by the earlier work of Majmudar³³ and Puckett.³⁴ PEGS is written for the MATLAB platform and is available for download and development on the GitHub open source software platform.⁵⁴ It makes use of Matlab’s built-in parallelization toolbox, making it suitable for running on high-performance computers to reduce computational time.

In general, it is necessary to have two images for each experimental configuration: one taken in unpolarized light for detecting particle locations, and one taken in polarized light for measuring the forces on those particles.¹⁶ Several possible ways to do this are discussed in Section VI D.

Throughout this section, we illustrate the process using the sample set of images shown in Fig. 10. The top image is presented as taken in RGB color by the camera. This comes from a setup illuminated by two monochromatic light sources: red unpolarized light and green circularly polarized light (see Fig. 16). This configuration allows the camera to automatically collect the two necessary images within a single color frame.

The PEGS inversion process comprises three programs. The pre-processor program detects particles, identifies neighbors, and validates which of those neighbors are contacts. This program operates on both the polarized and unpolarized images. The positions of the particles and contacts are then used by the fringe-inverter program, which performs the image-inversion process by which the photoelastic fringe patterns from a circular polariscope are used to determine all vector contact forces. This process operates on polarized-light images only. Finally, the post-processor program contains a suite of tools for error-correction and visualization. This provides a means to check the output for quality and then use this feedback to adjust the program’s parameters before submitting batch jobs to analyze a full dataset.

The first step is to obtain the locations and radii of all particles in the packing, from an image such as Fig. 10(b). From experience, we find that this needs to be done with an accuracy of $≲ 0.05 d$ in order to successfully invert the fringe pattern; the quality of the fringe-inversion fits will improve significantly with an increased centroid accuracy. We have found that the MATLAB imfindcircles() function provided in the Image Processing Toolbox is able to achieve the necessary level of precision; other Hough transforms provide similar performance. Several algorithms are available within that function, and for the results shown below we used the Hough transform (see Ref. 55 for a review), applied separately for each particle size. The results shown in Fig. 11(a) use red/blue circles to mark the locations of all particles with large/small radii, respectively. Many other particle-finding routines, including convolution methods, could alternatively be used. For a review, see the work of Shattuck.³⁵ 

In order for the fringe-inversion process to work efficiently and effectively, it is important to have a high-quality list of force-bearing contacts. The initial list of candidate-contacts is created by first identifying all neighbors by a threshold distance between the centers of a pair of particles of known radii, located at $x → i$ and $x → j$⁠. A naïve definition of a contact would be to call two particle contacts if their separation was within some tolerance d_tol of the sum of their respective radii,

However, even with sub-pixel resolution, this method is inadequate for determining whether the contacts are force-bearing. For each disk, the number of neighbors n satisfying Eq. (18) is greater than or equal to the actual number of contacts z, and d_tol must be set generously to overcome uncertainty in particle positions.

Therefore, we perform a secondary screening step to determine which of the neighbors should be included in the subset of z contacts. The criteria for inclusion are that the photoelastic response in an area adjacent to the point of contact must be above a threshold value.^16,33 In the analysis that follows, we use a threshold based on the $⟨ G 2 ⟩$ response (see Eq. (1)), rather than the intensity. If both particles show sufficient response near the contact point, then the contact is declared to be load-bearing.

A sample comparison is shown in Fig. 11, where all of the small white and green circles indicate possible contacts, but only the green circles meet the threshold criteria for being a contact. This significantly reduces the possibility, shown in Fig. 11(b), that a fringe located near the outer edge of one particle cause a false positive. This thresholding method is sensitive to having sufficient resolution and contrast. As a practical measure, it is necessary to adjust the contact threshold and the size of the evaluated area around the contact point for each specific dataset. In the PEGS software,⁵⁴ we provide a verbose option in the program allowing the user to get a visualization of the detected contact network, as shown in Fig. 11(c).

In addition, uncertainty in which contacts should be included in z can be overcome after the fact: the force-inversion code can later determine that $f i ≈ 0$ (a null contact). This simply increases the parameter space over which the inversion process must search, from 2z − 3 to 2n − 3. Therefore, it is better to set d_tol and the G²-threshold to be generous in including likely contacts on the list. While this will result in an increased computational load, the force solver is free to set the force of any contact to zero, but it is not allowed to create new contacts.

Using the list of particle positions, radii, and contacts, we can iterate over all particles in the system to determine which set of $f → = ( f , α )$ forces from the general solution (Section IV B, Refs. 1, 50, and 51) would produce the observed pattern of isochromatic fringes. In the text below, we use the same naming of the location of the contacts β, forces f, and impact angles α as used in Fig. 9.

The inversion process proceeds one disk at a time, for which we already know the locations of the z contacts from Section V A. Importantly, the location of all contacts must be along the line between the two particle centers and is specified by the polar angle β. This is satisfied by default in the algorithm described in Section V A.

We generate an initial guess for each force magnitude f by first applying the gradient squared method from Eq. (1) in Section III to the entire disk. We then distribute this total force among the z contacts in proportion to the value of G² within a small area around each contact point. To make an initial guess for the angle α of each force, we select a value for which the given contact locations β and guessed values of f would collectively be in force balance (see Section V C). This is included as an optional step, since force balance will not be applicable for systems for which there is significant dynamics.

Using the set of z values of $( f , α , β )$⁠, we solve for the stress field $σ ( x , y )$ on a grid chosen to match the resolution of the original image. Optionally, it is possible to downscale the grid for faster computing time at a reduced resolution. Together with Eq. (6), this creates a pseudo-image $I ps ( f , α , β )$ of the isochromatic fringes. We compare the guess to the original image I_obs by the residual function

where x and y are the spatial coordinates of the images. This residual function is suitable for optimizing via any standard algorithm, with the values of $f 1 … z , α 1 … z$ are varied in order to find the minimum value for e, as was first done by Refs. 16 and 33. We have had good success using the Levenberg-Marquardt^56,57 implementation contained within the MATLAB lsqnonlin() function. Note, however, that this will only find a local minimum in the $( f , α )$ landscape and is therefore sensitive to the quality of the initial guesses. Because the stress field calculation is computationally expensive, we calculate the σ solutions for each grid-point in parallel. This option, available natively in MATLAB, enables the code to run efficiently on multi-core/threaded computers or high performance computing clusters.

A sample outcome for the image in Fig. 10 is shown in Fig. 12. The left column contains the original image, and the right column contains the pseudo-image resulting from the optimization procedure. Visually, the agreement is excellent; Section V D will show how to use Newton’s third law to provide a quantitative evaluation. Note that the number of fringes rapidly increases in the vicinity of strong contacts, as shown in the single-disk pseudo-images in Fig. 10. Because this information is not available for comparison in the observed images due to the optical quality of the particles, we find that the accuracy and efficiency of the results are improved by fitting only the innermost 0.95r of each disk rather than using all pixels up to the outermost edge.

For both the initial guesses (mentioned above), but also for all subsequent optimization steps, it is possible to improve the quality of the guesses by making use of force and torque balances. Because these additional constraints reduce the number of unknowns in $( f , α )$ from 2z to 2z − 3, we have found that this optional step greatly improves both the accuracy of the fit and the time to convergence.

For z contacts, the force and torque balance constraints are given by

For z = 2, the above equations reduce to

By examining Fig. 9, α can be determined geometrically by the internal angles in the triangle O₁O₂C. The magnitude of the contact force f is the only free parameter.

The initial test of the algorithm for minimizing $e ( f , α )$ is used on a set of diametrically compressed disks, where I_obs is the image of the disk and f_applied is known. In Fig. 8, the result of the best fit $I ( f , α )$ is compared with I_obs. For z = 2, initial guesses quite far from the true $( f , α )$ will still converge; however higher fringe numbers require better initial guesses. The error in the magnitude of contact force was typically $≲ 5 %$ and did not depend on the magnitude of the force (see Fig. 13).

The solution for $z > 2$ requires a bit more work. The force and torque constraints in Eq. (20) reduce the number of free parameters by 3. To solve for these constrained parameters, one could do a coordinate transformation from $( f , α )$ to the Cartesian (f_x, f_y); however this is not necessary. The solution requires a coordinate transform of the contact locations β by rotating each contact so that $β z = 0$ (i.e., taking $β i = β i − β z$⁠). After, $( f , α )$ are found, we simply rotate the coordinates back. Solving Eq. (20) for f_z−1, f_z, and $α z$⁠, we have

This will ensure a set of $f 1 … z$ and $α 1 … z$ that is both force and torque balanced.

The force-inversion code provides the full set of $f → = ( f , α )$ values for each contact, and the residual score e (Eq. (19)) for each disk. Using Newton’s third law, it is possible to quantitatively evaluate the reliability of the results shown in Fig. 12. To illustrate this on the small number of example particles, Fig. 13 compares the force fit on one side of a contact to the (presumably identical) force fit on the other side. Ideally, all contacts would have $f A → B = f B → A$⁠; therefore deviations from this line indicate errors in the inversion process. We observe that these are of a consistent relative magnitude throughout the range of forces observed.

In addition, since each contact is subject to two independent optimizations, it is possible to flag contacts with unreliable results. Individual values, or entire disks for which the residual score e is too high, can be automatically or manually replaced by the value from the opposite-contact. This feature also allows the data to be transformed so that all contacts (rather than all particles) experience force-balance, by simply assigning each contact a force which is the vector average of the two independent measurements.

Visualizations such as Fig. 12 also provide feedback about improving the resolution and lighting. For example, the force-inversion technique is limited by resolution and contrast in I_obs and also governed by the number (and quality) of initial guesses provided to the algorithm. In the low force limit, I_obs is dark in the center of the disk and contacts appear as small dim spots near the perimeter. Greater image resolution and contrast can improve the accuracy of the contact forces found during optimization. Similarly, the resolution of the image also limits how well large forces are calculated, where fringes can become pixelated. Large forces (and large z) are more computationally expensive, and more/better initial guesses are required to find the correct local minimum for $e ( f , α )$⁠.

Several commercially available materials are available with a uniform photoelastic response. These generally come in various stiffness and are sold in one of the two main formats: sheets, or castable liquids. Methods suitable for both formats are discussed below. Because of the expense associated with all types of particles, it is advisable to make only a dozen particles in the first batch and check them between polarizers to determine whether the quality is good enough. For instance, frozen-in stresses often take the form of a bright ring around the outer edge of the particle when viewed in a polariscope. It is sometimes possible to eliminate/reduce these stresses by slow, gentle heating over a few days or weeks.

In designing the particles, several geometrical considerations are important. The granular materials will be most amendable to experimentation when all particles (disks/cylinders) stay in a single plane. This can be achieved by making particles which are wider than they are tall so that they do not tip over when sheared. When calculating the amount of material needed, it is helpful to cut at least 25% more particles than seem necessary; particles inevitably become lost or damaged over time, and it is difficult to later make identical particles. These are then best replaced by particles made in the same batch as the others to avoid the necessity of particle-specific calibrations.

In humid climates, photoelastic particles will experience some degree of adhesion to themselves and their substrate due to liquid bridging. One method of reducing this effect is to dust the experiment with a very small amount of baking powder, a technique pioneered by the Behringer Lab.

The following are common sources of sheet stock. Vishay PhotoStress Plus⁵⁸ is a proprietary material with excellent uniformity. Nearly clear urethane sheets are available from various industrial sources⁵⁹ and are a cheaper alternative. Both of these options come in sheets of various stiffness and thickness and have optical clarity and uniformity suitable for quantitative measurements. Sheets are available with a small amount of dye, which can improve particle-tracking. Although clear cast acrylic (Plexiglas^TM) is widely available, it is generally too stiff (low $F σ$⁠) for most granular experiments. In addition, it is weakly dichroic (it exhibits polarization-dependent light absorption) and is therefore not suitable for quantitative photoelastic analysis.

For any sheet material, several cutting methods are possible depending on the technical capabilities of the machine shop. For circular particles, a skilled machinist can create a custom tool (see Fig. 14) which acts as a spinning disk cutter to create circular particles of a specified diameter and with little waste. Alternatively, a small-sized computer-controlled mill can trace the outside of individual particles and make repeated custom (non-circular) shapes. However, the sidemill tool will waste the material in proportion to its diameter. In either case, the machinist will need to lubricate the cutter with dish soap rather than the usual machine oil, which degrades the polymer. With good lubrication, it is possible to create particles with smooth edges all the way down, and no frozen-in stresses. It is also possible to use a water jet cutter to make custom shapes with little waste if care is taken to minimize the presence of a notch or flap at the starting/ending location.

Two simple cutting options, which initially seem promising, are less practical than they might appear. The use of a punch is ill-advised since these leave a lip at the bottom of each particle. Unfortunately, both PhotoStress and urethane break down due to the heat of a laser cutter and release toxic gases, so these do not work. For a low-cost option, it is instead advisable to simply use a band saw to cut strips and then particles as was done in early experiments.⁴³ 

Both Vishay PhotoStress Plus⁵⁸ and urethane^60,61 are also available as a two-part castable liquid. Both formulations are available in a variety of stiffness and can be clear or lightly dyed. The companies provide protocols for making and releasing molds. Note that in order to achieve bubble-free particles, it is necessary to cure the material under vacuum. Castable liquids have allowed for the creation of particles of unusual shapes and also with magnetic inclusions.⁶² 

Two main polariscope configurations are possible: the light can pass through the granular material once on its way to the camera (transmission) or twice (reflection). The most common configuration is the transmission polariscope, shown schematically in Fig. 15. If the experiment is horizontal, the upper plate could be omitted but is often included to prevent out-of-plane buckling of the granular layer. The pair of polarizers may be either of opposite chirality (darkfield, more common) or the same chirality (brightfield).

In principle, it does not matter whether the polarizers are inside or outside of the supporting layers. In practice, however, if the supporting layers are themselves made of a birefringent material (e.g., acrylic), then the polarizers must be placed directly adjacent to the particles to avoid measuring photoelastic effects present in the supporting layer. The camera and lightbox will be placed on either side of this sandwich. Such a configuration has the additional difficulty, however, that the polarizers will become scratched over time.

There are several advantages in arranging the system so that the camera-side polarizer is attached directly to the camera. First, this avoids the expensive purchase of a second, large polarizing filter. Second, easy-to-use polarizing filters are available for many commercial camera lenses (see Section VI C for details).

In some cases, a transmission polariscope is undesirable when the supporting layer must be opaque. This is the case, for instance, for a porous membrane (frit) providing air-levitation as in Ref. 20. To accommodate such situations, it is possible to build a reflection polariscope which requires optical access from only one side of the granular material. This is shown schematically in Fig. 16.

To create reflecting particles for use in a reflection polariscope, a simple method is to coat one side of the particles with silver spray paint (e.g., Rust-Oleum Mirror Effect, used by the Behringer Lab). These particles will reflect the incident light and simultaneously flip its polarization. Therefore, a darkfield polariscope results from configuring both the illumination source and the camera with the same chirality polarizer.

Two formats of polarizers are typically available, large sheets and camera-based filters. Care must be taken to identify which polarization is specified and to select left/right circularly polarized products as desired. For large sheets, current sources include API Optics, http://www.apioptics.com/,³⁸ Edmund Optics, http://www.edmundoptics.com/,⁶³ and Alight Polarizers, http://www.polarization.com/.⁶⁴ Photography shops sell camera-based filters that are sized and mounted for use on standard camera lenses.

Care must be taken when installing the filters in the polariscope, so that their orientation is correct. A simple check for a darkfield (lightfield) polariscope is to place them face-to-face and rotate one polarizer with respect to the other. When the quarter-wave plates are both facing inwards, rotating either polarizer will always maintain a dark (light) view. If either polarizer is flipped inside out with respect to the other, then the view will go alternately light and dark as one polarizer is rotated with respect to the other. Note that circular polarizers built for standard photographic purposes will have their quarter-wave plate facing the camera sensor; this is the wrong orientation for a polariscope analyzer (see Fig. 2). To remedy this, it is usually possible to unscrew the mounting frame and flip the filter over within its holder so that the quarter-wave plate faces outwards from the camera. Recall, also, that for a reflecting polariscope, having a polarizer and analyzer with matched chirality (LCP vs. RCP) gives a darkfield image. A final difficulty is that the specific chirality of photographic filters may or may not be provided at the time of purchase.

Since the fringe pattern for birefringent particles depends on the wavelength of light through both the stress-optic coefficient and rotation of the light phase (see Eq. (5)), it is best to record monochromatic light for quantitative results. There are two main ways to do this, either by using a monochromatic light source or by filtering polychromatic light at the source or at the camera itself. A very simple method, suitable for non-quantitative studies, is to simply visualize only one of the three RGB (red-green-blue) channels of a full-color image, a technique pioneered by Ref. 16, by applying multiple filters in succession, and now possible to achieve using single-images. Because these channels have a bandwidth of ∼100 nm, quantitative studies using this type of color separation are wise to still use monochromatic light. This could also be achieved via a monochromatic filter which screws directly onto the photographic lenses, available at any camera supply store.

Different lighting solutions are appropriate depending on whether the system is a transmission or reflection polariscope. In either case, it is better to achieve uniform lighting before collecting data rather than requiring a post-processing step to equalize the brightness levels across all images. In all cases, it is preferable to use low-temperature lights such as light-emitting diodes (LEDs) or fluorescent bulbs rather than halogens or incandescents which will heat the experiment and change the material properties of the particles and/or their supporting layers.

For a transmission polariscope, several commercial solutions are available. Common types of lightboxes include artists’ copy-boards, sign and menu boards, and medical x-ray lightboxes. LED-based boards are typically quite thin and have no AC flicker when running on rechargeable batteries. A doctor’s lightbox, while commonly available used at low cost, uses fluorescent lightbulbs which have significant 60 Hz flicker unless they have been specially designed to eliminate this or use an electronic ballast that upconverts the frequency to 20 kHz. It is also straightforward to build a custom lightbox out of low-cost strips of LEDs which can be cut to length to suit any shape experiment. For a reflection polariscope, many other options are possible, including colored LED spotlights (either decorative or theatrical) which provide monochromatic light.

To perform particle-tracking and photoelastic measurements at the same time, it will be necessary to have at least two different sources of light. Typically, a polarized, monochromatic light source provides the photoelastic image via a polarizing filter placed on the camera. A non-polarized light source, which will still be predominantly transmitted through a camera-mounted polarizing filter, provides other measurements such as the particle positions. There are numerous methods for achieving this goal: multiple cameras with different filters and post-processing to align the images, electrically or mechanically triggering different lights in sequence; placing/removing one of the polarizing filters; or using RGB color-separation of simultaneous sources within a single camera. This last option is illustrated in Fig. 10, for which the two images came from the red, green LED configuration shown schematically in Fig. 16.

This configuration leaves the third color channel (blue) available for supplemental use. For example, it is possible to draw features with an ultraviolet marker that will glow blue when exposed to black light. This technique has been used to monitor the rotation of particles,^18,28 to identify and track a sub-population of particles,²⁰ and to outline the edges of very clear particles to improve their detection.

K.E.D. is indebted to Bob Behringer, in whose lab she was first introduced to the use of photoelastic materials and who has championed the technique for many years. The photoelastic experiments by the three authors have been supported by the National Science Foundation (Nos. DMR-0644743 and DMR-1206808) and the James S. McDonnell Foundation. This review article began as a resource sheet and talk for the NSF-funded (No. EAR-1537902) workshop “Analog Modeling of Tectonic Processes” held in Amherst, MA in May 2015 and was further developed for the Spring School “Imaging Particles” held in Erlangen, Germany in April 2016, funded by the German Science Foundation (DFG) through the Cluster of Excellence “Engineering of Advanced Materials” No. EXC-315. The authors thank Arne te Nijenhuis for his work on Fig. 2, Amalia Thomas and Zhu Tang for their work on Fig. 8(b), and Amalia Thomas for her work on Fig. 14(b), and to Amalia Thomas and Nathalie Vriend for describing their techniques for casting urethane particles.