Analysis of particle contact using frustrated total internal reflection

Within the field of soil mechanics a continuum assumption is generally adopted in order to avoid the complications of modelling micro-mechanical behaviour. However, certain constitutive behaviour can only be explained by investigating particle level interactions. Numerical investigations, such as those using the Discrete Element Method (DEM) to model soil particles as clusters of spheres, have delivered a greater understanding of the micro-mechanical behaviour. One of the limiting factors in current DEM approaches is modelling of the particle–particle or particle–surface contact behaviour. Hence, an experimental methodology has been developed and used to study particle–surface contact behaviour. The experimental methodology involves loading particles onto a piece of sapphire glass and observing the resulting contact area. In order to distinguish between the contacted area and the rest of the particle, the principle of frustrated total internal reflection and evanescent waves was used which results in only objects in very close proximity to the glass being illuminated and visible. This methodology hence allows the number of contacts and the area of those contacts to be tracked during loading and over time. This paper presents the validation of the experimental methodology by comparing the observed contact behaviour of plastic beads against Hertzian contact theory. In addition, the results from tests on sand samples are presented which show a density of 0.40 and 0.80 contacts per $$D_{50}^2$$D502 for coarse and fine grained sand respectively at an isotropic stress state which subsequently increases to 0.90 to 1.00 contacts per $$D_{50}^2$$D502 at peak deviatoric stress. It was also found that the fine sand particle contacts carried a maximum load of approximately 0.27 N per contact whereas the coarser sand was able to carry substantially higher loads.

Abstract Within the field of soil mechanics a continuum assumption is generally adopted in order to avoid the complications of modelling micromechanical behaviour. However, certain constitutive behaviour can only be explained by investigating particle level interactions. Numerical investigations, such as those using the Discrete Element Method (DEM) to model soil particles as clusters of spheres, have delivered a greater understanding of the micromechanical behaviour. One of the limiting factors in current DEM approaches is modelling of the particleparticle or particle-surface contact behaviour. Hence, an experimental methodology has been developed and used to study particle-surface contact behaviour. The experimental methodology involves loading particles onto a piece of sapphire glass and observing the resulting contact area. In order to distinguish between the contacted area and the rest of the particle, the principle of frustrated total internal reflection and evanescent waves was used which results in only objects in very close proximity to the glass being illuminated and visible. This methodology hence allows the number of contacts and the area of those contacts to be tracked during loading and over time. This paper presents the validation of the experimental methodology by comparing the observed contact behaviour of plastic beads against Hertzian contact theory. In addition, the results from tests on sand samples are presented which show a density of 0.40 and 0.80 contacts per D 2 50 for coarse and fine grained sand respectively at an isotropic stress state which subsequently increases to 0.90 to 1.00 contacts per D 2 50 at peak deviatoric stress. It was also found that the fine sand particle contacts carried a maximum load of approximately 0.27 N per contact whereas the coarser sand was able to carry substantially higher loads.
Keywords Discrete element modelling Á Contact behaviour Á Optics Á Image processing

Introduction
Traditionally geotechnical design has focused on capacity (failure) calculations with the application of large factors of safety and hence overly conservative designs. However, with the ever growing demands for infrastructure, there is increasing motivation to undertake performance based design [1] which uses advanced constitutive models to make accurate predictions of what ground movements will actually occur during and after construction. This approach allows for more efficient design and the opportunity to construct in areas which, with capacity based design, would have previously been uneconomical.
The development of advanced constitutive models relies on a fundamental understanding of soil behaviour, both at the macro and particle scale. Particle scale investigations have mostly focused on the application of the Discrete Element Method (DEM), which typically uses spheres to develop a simulated model of a real soil matrix. One of the challenges associated with such simulations, which has drawn attention from researchers [2][3][4][5], is accurately modelling the contact and breakage behaviour at the particle-particle and particle-surface contact points. In order to validate the implemented contact models, the simulated response of the overall soil matrix, for example the stress-strain behaviour, can be compared with data from experiments. However, it is difficult to determine the source of deviations between the simulated and measured responses as the particleparticle interaction behaviour, which controls the overall soil matrix response, cannot be extracted from experimental investigations. Therefore, an ability to compare particle scale behaviour between experimental and DEM simulations would allow for significant improvements in the accuracy of DEM simulations. Particle scale features which could be compared include the coordination number, the statistical distribution of normal and tangential forces between particles and the contact areas.
Taking a different approach, finite element modelling has recently been applied to virtualized single grain particles subjected to loading [6] in an attempt to better understand the contact mechanics. However, as with all numerical simulations, such modelling is only as accurate as the input parameters used and hence does not negate the need for experimental based insights.
In addition to understanding the response of a soil matrix to an increasing load, the response under a constant load, due to creep, is also of interest. An experimental means of examining how contacts change over time when subjected to a constant load would provide an insight into the mechanisms behind creep in granular soils.
This paper will present a technique which is capable of providing the experimental insights described previously, including providing information about the contacts per unit area, the area of each contact and the average force per contact.
Contact tracking has previously been the topic of research in areas such as tribology however the developed experimental methodologies have not yet been applied to the study of granular matter within the field of soil mechanics. A common experimental approach to examine contact behaviour has been to make use of the principle of frustrated total internal reflection and evanescent waves [7][8][9]. This approach results in only the parts of an object in very close proximity to a contact plane to be illuminated by the light source.
In the case of the research presented in this paper, soil grains are loaded against a pane of sapphire glass. Light is passed into the glass from the side and hotspots are produced at all locations where the grains and glass are in contact. By imaging the variation in size of these hot-spots with increasing load, the load-contact area relationship can be determined.
In this paper, the experimental device will be described and its functionality validated by comparing the force-contact area results from single beads of two different plastic materials with predictions using Hertzian contact theory. In addition, a cylindrical sample of beads will be tested resulting in multiple contact points across the glass window. Each of these contacts will be analysed and the contact force determined. These individual forces will then be summed and compared to the global externally applied load. Finally, real sand grains will be investigated, with the contact area and contact density per unit surface area being measured as a function of the grain size, applied load and packing density.

Frustrated total internal reflection (FTIR)
Prior to describing FTIR, the effect of refraction and total internal reflection (TIR) must be explained. Refraction is an optical effect in which light passes with a non-perpendicular angle from one medium to another if the two materials have a different density or refraction index. The angle at which the light exits the interface is higher than the entering angle if the second medium has a lower refraction index than the first medium (as shown in Fig. 1). Total internal reflection occurs if the entering angle becomes greater than a critical angle at which the exiting angle would be theoretically higher than 90 degrees. Therefore, light cannot pass into the second medium and is reflected back through the first medium with the same angle. This critical angle is given by Eq. 1: where n 1 and n 2 are the refraction indices of the first and the second media respectively. During total internal reflection (TIR), a phenomenon involving what are termed 'evanescent waves' occurs. These waves are produced at the boundary of the medium at which the total internal reflection is occurring. If an object is located near the interface or touching the interface, the evanescent waves will interact with the object and light intensity is transmitted. This intensity can be either absorbed, reflected, refracted or diffused by the object depending on its own optical properties.
The properties of these evanescent waves depend on factors such as: the wavelength of the light; the distance to the object; the incident angle; and the refraction indices of both media [10]. The magnitude of the interaction of these evanescent waves with the object in close proximity decreases exponentially with distance between the interface and the adjacent object. The intensity decreases to 36% (exponential of À1) of the intensity at the interface at what is termed the 'penetration depth'. In this research, the light used is a blue LED strip giving a peak of wavelength of 450 nm and the high density medium is a sapphire window with a refractive index of 1.76. Therefore the penetration depth ranges from 100 to 50 nm for incident angles ranging from 45 to 90 respectively.
The interaction between the evanescent waves and the adjacent object results in only the part of the object in close proximity to the interface being visible while the remaining part of the object remains unilluminated (Fig. 2). With the size of the particles investigated in this study being significantly larger than the penetration depth magnitude, the part illuminated by the evanescent wave effects can be considered to represent the contact area.

Experimental device and data processing
At the core of the experimental setup is a 50 mm diameter, 5 mm thick sapphire window. The particles are located above the window and loaded onto it. This loading is applied by a load frame in the case of a single particle. However, for the larger cylindrical samples (50 mm diameter, 95 mm height), which are contained within a flexible latex membrane, in addition to the application of axial load, the air pressure is reduced within the sample hence replicating more traditional triaxial loading conditions. These two configurations are shown schematically in Fig. 3.
The experimental device includes the following components (parts labelled in Fig. 3): 1. A sapphire window on which the samples are placed; 2. An LED strip fitted around the circumferential edge of the sapphire window in order to provide the light required for the FTIR; 3. A digital camera, with an exchangeable lens, placed beneath the window to provide images of the hot-spots produced from the FTIR (i.e. the contact areas); 4. Exchangeable lens with a 10Â microscope magnification lens for examining single particles and a macroscopic lens for examining larger samples. The interface used is a sapphire window from SP3 Plus. The sapphire material was selected due to its relatively higher hardness than quartz. Experiments have shown that, although a standard quartz glass window is suitable for softer particles such as the acrylic and polypropylene beads used during the validation process, it cannot be used for sand grains due to permanent damage to its surface. When damaged, the surface of the glass will scatter light where damage is present and hence affect the ability to accurately decipher contact points.
In order to extract the contact area from the obtained images, an intensity threshold cut-off is applied to the image within MATLAB. The reaming pixels provide an un-calibrated measure of the contact area. This is subsequently calibrated using millimetre paper which is placed on top of the sapphire window where objects come into contact as illustrated in Fig. 4. The lens distortion is corrected by applying a 4th order polynomial function map to straighten and calibrate the image. By measuring the force corresponding to each image, a relationship can be established between the contact surface area and the applied force. Subsequently the contact area can be used to back-calculate the applied force.
For single particle experiments, the loading was carried out with a loading rate of 1 N/min. For larger

Validation with Hertzian model
Prior to investigating the behaviour of real sand grains, a process was undertaken to confirm the validity of both the equipment and analysis process. This involved the loading of 3/16 00 (4.8 mm) diameter spherical black gloss acrylic and white matt polypropylene beads. The applied load was recorded and images taken. The obtained results were then compared to a Hertzian model prediction. The Hertzian model [11] predicts the load-contact area relationship between a sphere and a planar surface.
The relation between the applied force (F) and the contact area (S) is given by where K, the stiffness, is defined as where R is the radius of the spherical bead, and E Ã is defined by E and m are the Young's modulus and the Poisson's ratio respectively and the subscripts represent each body (the particle and the plane). In this case, the sapphire window has a significantly larger Young's modulus than the plastic beads. Consequently, where the subscript p represents the plastic beads' elastic properties.

Acrylic bead
An example image of the contact between an acrylic bead and the sapphire window is shown in Fig. 5 (as a negative for clarity). The force applied in this image was 3.3 N. Three black acrylic beads were tested separately and the images processed as described previously. The results, and a comparison with the Hertzian prediction, are shown in Fig. 6. The Hertzian model is fitted to the three results, by optimising the choice of the stiffness parameter (K), and as shown in Fig. 6 the model provides a good correlation with the experimental data. The optimised value for the stiffness parameter, K, was 265 N=mm 3 and therefore, E Ã ¼ 2:64 GPa which is in accordance with the order of magnitude quoted for this material (E ¼ 2:76 GPa to 3:30 GPa). A typical value for the poisson ratio was adopted for calculation of the Hertzian curve (m ¼ 0:37). The consistency of results between the three separate acrylic beads, combined with the agreement with the Hertzian model validates both the experimental apparatus and the subsequent image processing methodology.

Polypropylene bead
As a comparison, white polypropylene beads were also tested. Figure 7 shows an example of an image, in negative, of the contact point created. The image shows that the contact is not as clear as for the acrylic bead. Therefore, the surface area measurement is more sensitive to the adopted grey level threshold. This optical behaviour for this bead type makes it less suitable for conducting this type of experimental analysis compared to the acrylic beads.
Measuring the contact area as a function of the force applied gives the results shown in Fig. 8. The results show that the different polypropylene beads do not give perfectly consistent results. The Hertzian model can still be fitted to the experimental data, giving a stiffness of K ¼ 90 N=mm 3 . Therefore, E Ã ¼ 0:9 GPa which is close to the range provided by the manufacturer (E ¼ 1:1 GPa to 1.3 GPa). A poisson ratio of 0.42 was adopted for the polypropylene beads. It is noteworthy that the load-contact area relationship differs between the loading and unloading stages due to polypropylene demonstrating viscoelastic behaviour. During the testing of Sample 3, the loading was paused for 3 min before unloading and an increase in contact area was observed, further demonstrating the potential impact of viscoelastic behaviour on such testing. The experiment with polypropylene beads was only carried out to illustrate the fact that some plastic materials demonstrate viscoelasticity which can impact the results from these types of experiments. As this stage of the experimental research is focused on validation of the FTIR methodology, such viscoelastic behaviour reduces the usefulness of the data. Little to no viscoelasticity was observed for the acrylic beads as shown in Fig. 6 which makes acrylic a better candidate for this validation process.

Sample with multiple contacts
Following the successful validation of the apparatus and analysis method with single particles, cylindrical samples of 50 mm diameter and 95 mm height were tested to visualize the behaviour of a group of beads under load. Across all samples a consistent voids ratio of 0.65 was maintained.

Acrylic beads
For the acrylic bead samples, only a quarter of the sample height from the sapphire window was filled with acrylic beads with the remaining volume being made up by polypropylene beads. Given the sole measurement of load comes from the area measurements at the interface, the presence of these polypropylene beads will not have an impact on the obtained results. Indeed, the load transferred from the load cell as an external action is reacted solely by the acrylic beads on the interface. Figure 9 shows an example negative image of the contact points formed between the acrylic beads and the sapphire window. It should also be noted that only The presence of dust particles, despite careful cleaning of the beads during sample preparation, and reflections around the edges of the beads results in the contact areas not being as clearly defined as when individual beads were tested. Nonetheless, it is still possible to detect the contact points manually and subsequently automatically track them through the image processing methodology adopted.
The area of each contact point is measured within the field of view while the load cell force and vacuum pressure are recorded. The contact areas are converted to contact forces using the Hertzian contact model and the calibrated stiffness from single beads. Figure 10 shows the cumulative force provided by the contacts working from the centre to the outside edge of the measurement area, i.e. the lowest line (blue) indicates the force transmitted through the single central-most particle and the uppermost (green) line indicates the cumulative force transmitted through all particles within the measurement area. As was explained earlier, the field of view for the camera is smaller than the total sapphire glass area and hence, so that direct comparisons can be made, the total force, as measured by the load cell, in Fig. 10 has been reduced to represent the force applied on the measurement area (red line). This result shows how accurate this measurement technique is and how sensitive it is by clearly tracking the rapid changes in force which were imposed during the loading sequence.
It is also clear from Fig. 10 that there is no consistent pattern to the variation in contact force across the measurement area, for example it is not apparent there are larger contact forces near the centre and smaller forces at the outer edges. In fact the variation in contact force between particles across the measurement area would appear completely random. This indicates the samples were homogenous in their preparation and proximity to the boundaries of the sample did not impact the behaviour of the contacts. Figure 11 shows the distribution of the contact forces at different stages of loading which further illustrates clearly this random distribution.
The sum of all contact point loads divided by the measurement area compared to the applied pressure, as measured by the vacuum regulator and load cell, is shown in Fig. 12 for two identical tests on the acrylic beads. It is clear from Fig. 12 that the results demonstrate consistency and repeatability, further validating the experimental approach.

Polypropylene beads
Two samples of pure polypropylene beads were also tested. Figure 13 shows an example of an image in negative. The image is clearly better than for the sample of acrylic beads. One reason is that the beads are white and hence the camera exposure time is smaller than for the black acrylic beads leading to the dust particles being less visible. In addition, the white polypropylene beads provide a much brighter contact point than their black acrylic counterparts leading to Sum of contact forces/measured area (kPa) Test 1 Test 2 Fig. 12 Sum of all contact forces divided by the observing area as a function of to the total applied pressure improved contrast, which assists in the thresholding process. Another positive aspect of the polypropylene beads is that they are matt and not reflective-there is no reflection of the bead edges observed in the images, therefore facilitating easier contact point identification and post processing. The results show that the calibrated contact area cumulative force is significantly smaller than the measured total force using the load cell, even when accounting for the difference in the measurement areas. This is due to the viscoelastic behaviour, discussed earlier, which affects the load-contact area relationship and hence impacts on the ability to use polypropylene bead for such an experiment.
The samples are prepared using the same methodology as the plastic beads, with the sands being washed to remove any dust particles in advance of sample preparation before being attached to a piece of rubber which is subsequently attached the bottom edge of the low stiffness spring shown Fig. 3.
With sand particles having more complicated geometry than simple spheres, Hertzian contact theory may not necessarily apply. It is therefore difficult to provide an accurate measurement of the contact forces of each particle in a whole sample. Instead, it is of interest to use the now validated method of measuring contact area to study other characteristic behaviour of the single sand grain. Several grains of fraction A and fraction B sand have been tested and a typical image is shown in Fig. 14. It was observed that although three contact points are needed for a single particle to be stable, during loading (due to support from the top of the particle) the number of contact points may vary. Figure 15 shows six different sand particles tested and analysed for their contact surface area. The curves plotted represent the surface of each contact point (i.e. not the cumulative contact area for the particles). The results show that the contact areas are in the order of magnitude of 10 À3 mm 2 . The results also show a balancing of forces between contact points. As load increases, some contact areas may decrease but others increase at the same time.
By considering the data from the first Fraction A grain at a load of 1 N (Fig. 15, top left), it is possible to evaluate the radius of curvature of contact points.  From the elastic properties of the sapphire window and the quartz silica sand included in Eq. 4, a value of E Ã ¼ 50 GPa can be obtained. The two contact areas at 1 N force are 0:5 Â 10 À3 and 1:8 Â 10 À3 mm 2 . If both contact points behave as Hertzian contacts and the elastic properties of the sand grain is homogeneous, the total force can be obtained from Eq. 6 (derived from Eq. 2 for multiple contact points).
where i is the index of the contact points, R is the radius of curvature, and S is the contact area. If we assume that both contacts points have the same radius of curvature, R 1 ¼ R 2 , then a radius of curvature of 1 mm can be back-calculated from Eq. 6. This value is close to the dimensions of the grain which has a diameter of approximately 1.5 mm, giving confidence on the reliability of the experiment and processing methodology.
In future developments it will also be possible to directly determine the radius of curvature by measuring the geometry of the grain around the contact points. This is made possible by using focus-variation microscopy: a device that can measure live topography whilst loading the grain sample.
Further investigations are required to understand the variability of the surface areas from one sand grain to another. The values depend more on the geometry of the contact region than on the actual size of the grain, with sand angularity and sphericity being the most relevant parameters to take into account. However, this investigation is beyond the scope of this paper but will form the focus of future research. Firstly, the transparency of the sands particles affects their optical behaviour-light passes through the sand grains and reflects back to the camera. In addition, the resolution of the camera is not high enough in the macroscopic scale to identify the contact surfaces of sand grains. Consequently, a thin plastic film was placed between the sapphire window and the sand sample facilitating both removal of reflections and provision of a more visible contact area on the camera. Additionally, its presence removes the possibility of dust gathering on the surface of the sapphire window during sand pouring. The thin plastic film is made from polyethylene with a thickness of 22:5 lm. Sand samples were prepared by pluviation at two different heights providing two different void ratios (e) for each sand fraction. The samples are isotropically consolidated by vacuum at 95 kPa and then sheared (axially loaded) until their peak deviatoric stress was reached. Typical images of both sand fractions are shown in Fig. 16.
The number of contact points per unit surface area were counted and associated with the total pressure applied on the sapphire plate. The results, given in Fig. 17, show that the dense samples have a higher contact point density than their loose counterpart, as would be anticipated. At pure isotropic consolidation of 95 kPa, Fraction A samples show a contact point density of 0.75 to 0.8 contact points per D 2 50 which is nearly double than for Fraction B samples (0.42 contact points per D 2 50 ). However, the Fraction B samples under increasing deviatoric stress gain contact points at a faster rate than Fraction A and the two different density samples reach similar values at peak deviatoric stress (i.e. 0.90 1.00).
By dividing the total vertical force on the sample by the total number of contact points (extrapolated from the contact point densities quoted above), it is possible to evaluate the average force for each contact point (Fig. 18). The Fraction A samples show a steady increase in contact force. However, Fraction B samples show the average contact force increases until it reaches a cap value of approximately 0.27 N. The effect of capping can be seen in Fig. 17 where a change in slope occurs for both Fraction B curves at 175 kPa correlating to the point at which the rate of change in number of contacts increased. This cap in contact force is a new observation in the field of soil mechanics and requires further research to explore fully, especially given it is observed only in the Fraction B sand. A notable difference between the two sands is that Fraction B is a more rounded sand grain and thus will have reduced interlocking between grains. Consequently, as the particle to particle contact loads in the soil increase during shearing, there is an increased chance of slippage occurring at particle contacts in the Fraction B sand. This slippage will result in redistribution the loads carried by the 'just slipped' contact to other contacts which are under lower load. To explain in more detail, future research will be needed, involving both experimental and DEM studies of simplified samples of spherical and nonspherical shaped particles of the same size.

Conclusions
The optical principle of frustrated total internal reflection was used in this work to measure contact areas of different types of objects under applied loads. The first type of object examined was spherical beads of two different plastic materials: acrylic and polypropylene. These tests were carried out to assess the feasibility of the experimental methodology by comparing the results with Hertzian theory. Results conclude that acrylic is a better candidate than polypropylene for this validation process because of its lower viscoelasticity characteristic. The elasticity modulus derived from the experiment was in accordance with the range of the acrylic materials and the power law described by the Hertzian theory is closely followed by the experimental data. Subsequently, a group of acrylic beads were arranged as a conventional triaxial sample and tested under consolidated drained conditions with vacuum confinement. Each contact point at the bottom edge of the sample has been tracked and its area measured during the shearing process. These areas were then used to determine each contact force and it was found the sum of these individual contact forces equated to the externally measured applied load. The results from these tests showed that the distribution of local forces is not dependent of the location of the contact points.
The second type of object examined was sand particles. Single sand particles were loaded onto the sapphire window and results showed that the contact area is of the order of magnitude of 10 À3 mm 2 . Two different grain sizes were tested: Leighton Buzzard Fractions A and B. The smaller grains display a lower contact area but the difference between Fractions A and B is not significant. This suggests that grain size is not the only factor affecting the contact area. Other factors such as sphericity, angularity and radius of curvature at the contact should also be considered.
The two fractions of sand were also tested within a triaxial configuration at two different packing densities and the number of contacts per unit area was measured. Results show that for a higher packing density, the contact density is indeed higher, in line with expectations. Also Fraction A showed a higher contact density (0.72 contacts/D 2 50 for loose and 0.80 contacts/D 2 50 for dense) than for fraction B (0.40 contacts/D 2 50 for loose and 0.42 contacts/D 2 50 for dense) at pure isotropic stress conditions. This is due, in part, to the higher angularity of Fraction A. Contact densities increase with increasing deviatoric stress and approached to 0.9 contacts/D 2 50 at peak deviatoric stress. Average contact forces for each particle showed a steady increase at both densities of Fraction A whereas for Fraction B the loads increased until a cap of 0.27 N for both densities. This cap is not present for Fraction A and is hypothesised to be caused by the higher angularity of Fraction A which provides better interlocking. This novel observation is interesting for the study of grain crushing and requires further exploration.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.