X-ray 3D imaging–based microunderstanding of granular mixtures: Stiffness enhancement by adding small fractions of soft particles

To obtain the P-wave modulus of the biphasic sample, M, the ultrasound velocity (V_p = L_z/t_p) is measured knowing the travel time of the wave into the sample (t_p) and the sample height after deformation (L_z). The bulk density of the sample is determined for different rubber fractions as ρ⁰ = ν ρ_g⁰ + (1 − ν) ρ_r⁰, where ρ_g⁰ and ρ_r⁰ are the glass and rubber densities (ρ_g⁰ = 1,540 kg m⁻³ and ρ_r⁰ = 860 kg m⁻³ (40)) at preparation (p = 0), and ν is the rubber fraction of samples ν = V_r⁰/V⁰, where V_r⁰ and V⁰ are the volume occupied by the rubber beads and the volume occupied by all particles (V⁰ = V_r⁰ + V_g⁰), respectively. Thanks to the calculated mass density, ρ⁰, of the sample, and the experimentally obtained wave velocity, the longitudinal P-wave modulus is being determined by M = ρV_p² = ρ⁰(L_z⁰/L_z)V_p².

Fig. 1 illustrates the relationship between the small-strain stiffness, denoted as M, and the five vertical stresses, p, which are the resultant of the applied axial stress for samples with varying fractions of soft/stiff particles. Each experiment was repeated five times for every rubber content to obtain new configurations (i.e., new particle arrangements). The SD of the modulus is included in Fig. 1 for p= 80 kPa and 160 kPa with vertical bars.

As the applied vertical stress on the sample increases, its porosity decreases, making it denser and resulting in an increase in stiffness. Fig. 1 shows three different regimes highlighted in different colors as we move along the x-axis: i) Up to ν ≈ 0.3, the response is mainly controlled by stiff particles (red color in the background of Fig. 1). The tests conducted under uniaxial compression using an oedometer cell show a slight increase in the P-wave modulus, upon adding a small fraction of rubber particles; ii) adding more rubber beads leads to a drop of modulus for 0.4 ≲ ν ≲ 0.6. This indicates that the samples are weakened, as shown by the cyan color background in Fig. 1; iii) rubber occupies most of the total volume for ν ≳ 0.6 where the sample is fully controlled by soft particles. Thus, the stiffness remains constant even with the addition of more rubber particles (highlighted by the yellow color background in Fig. 1). This is due to the softness of rubber particles, which can deform significantly, leading to surface contacts instead of point-to-point contacts between particles. These samples are more like a porous rubber media with glass inclusions, and the packing seemingly resembles that of an inhomogeneous rubber block (64, 65). It must be noted that the maximum deformation applied to the samples (especially soft samples) is 10% of the initial sample height (L_z⁰= 80 mm), which is within the elastic limit of the rubber particles, so they do not deform permanently (30).

The oedometric cell experiment under uniaxial compression with zero lateral strains supports earlier research conducted by the authors using triaxial cells with isotropic stress control (40). However, it is still unclear how rubber strengthens samples for ν = 0.1 and 0.2. The microstructural mechanisms responsible for the transition from stiff- to soft-dominated networks have yet to be resolved. Therefore, we employ in situ µXRCT to provide a microinvestigation at the particle scale. Another representation of the modulus plotted against the applied stress is shown in SI Appendix, Fig. S1. In addition to the wave propagation tests performed under the stress-control condition, we performed complementary tests following a strain-control protocol for samples ν = 0, 0.1, 0.2, and 0.3, which is given in SI Appendix, Fig. S2.

Mechanical response investigations of samples using (ultrasound) wave propagation have allowed identifying the most interesting regime (0.1 ≲ ν ≲ 0.6) for in situ µXRCT experiments. The particles’ contact networks of some samples are studied next using the extracted information such as particles’ centroid and their contact locations with one another. To begin with, we classify particles according to the number of particles they are connected to. In Fig. 2A, we can see the distribution of particles based on the number of contacts they have. The y-axis (N_p) represents the number of particles, while the x-axis (c) shows the number of contacts within the respective subnetworks. These results are based on images obtained at p = 160 kPa and different rubber fractions ranging from 0.1 to 0.6. SI Appendix, Fig. S3 displays the analysis carried out on images collected at p = 80 kPa.

At first glance, all-network data appear to follow a Gaussian distribution function, as shown in Eq. 1. The function provides histograms that represent the number of particles based on the number of contacts they have:

In Eq. 1, c is the number of contacts carried per particle, a is a fitting parameter, σ is the SD with respect to the mean, and $\overline{C}$ is the coordination number of a sample, which is defined by

Dashed and solid black lines are the fitted plots using Eq. 1 for samples with ν = 0.1 and 0.6 scanned at p = 160 kPa, respectively. The distribution of particles and their contacts indicates that the samples almost coincide on the same function, regardless of particle characteristics. However, when more rubber particles are added (ν ≳ 0.4), the plots shift slightly to the right, due to the nature of rubber particles, which can deform and create more contacts with neighboring particles.

Next, the networks are split into glass (G–G, Fig. 2B) and rubber (R–R, Fig. 2C) subnetworks for each sample, where only particles belonging to a single species are taken into account. For completeness, we plot in SI Appendix, Fig. S4, the connectivity of glass and rubber where only the contacts between rubber and glass particles are considered and not the ones from the glass–glass and rubber–rubber subnetworks. Both N_p and c decrease as ν increases. On the contrary, the opposite trend is observed in the rubber subnetwork plots. Comparing the data of ν = 0.1 (red) and ν = 0.2 (blue), the number of particles with a high number of contacts is quite similar in the glass networks, i.e., the right tail of the glass network is very similar for both cases. In contrast, this is not observed for rubber networks where almost all rubber particles carry just less than three contacts for ν = 0.1.

Earlier studies have shown that particles carrying less than three contacts, so-called rattlers, are not mechanically stable and do not contribute to force networks since their few contacts are momentary (66–69). The networks that are considered in this study exist in a dense regime under gravity, and as a result, there are very few rattlers. This can be observed in Fig. 2A, where the plot shows that only a small number of particles (N_p) have less than four contacts (c). However, when considering only a rubber network with ν ≈ 0.1, the isolated clusters of rubber particles exhibit similar characteristics to rattlers in a loose packing. This is due to the fact that rubber has a stiffness that is three orders of magnitude smaller than that of glass, and as a result, rubber particles surrounded by glass, similar to an isolated “void,” do not contribute significantly to the force network in comparison to glass–glass contacts. Therefore, the overall response of the entire network is controlled by the glass particles. The main difference between a sample with ν = 0.1 and a pure glass sample is the deformability of the bulk due to the rubber islands, which reduces the length of the sample.

In contrast, samples with ν ≳ 0.2 contain more rubber particles with coordination above three. Some of these particles participate in the overall force networks by forming chains, loops, or clusters between themselves and glass particles, but the amount of rubber particles is yet not sufficient to influence the overall network. The sample with ν = 0.3 shows a qualitatively similar response as ν = 0.2. Nevertheless, the overall response of the network is weaker since stiff particles are replaced by soft ones, which leads to a decrease in M-modulus when ν ≳ 0.3.

Finally, to illustrate the evolution of networks, glass and rubber subnetworks are shown in Fig. 3 for two samples, ν = 0.1 (A and B) and ν = 0.5 (C and D), compressed under p = 160 kPa. Comparing (B) and (D) subfigures (rubber network), we observe many isolated rubber particles with almost zero or one contact which cannot form a chain/path for wave propagation. In contrast, a highly dense network has formed for the rubber network of ν = 0.5 in subfigure (D), indicating an increase in the coordination number and demonstrating the transition from a stiff- to soft-regime.

The next microparameter investigated here is the coordination number, Eq. 2, which reveals more quantitative information about the microstructure. In Fig. 4, we demonstrate the evolution of the coordination number, $\overline{C}$, of glass subnetwork, rubber subnetwork, and overall networks for samples with different rubber fractions at p = 160 kPa. If we take into account all the particles in the networks, we can observe that the coordination number (indicated by the green line) increases slightly as the rubber content increases. This is because the rubber particles are more deformable, which leads to the formation of more contacts between all particles due to the lower bulk volume. As rubber content increases, the glass subnetwork deteriorates further while the rubber–rubber network develops.

Earlier studies have shown that if the coordination number of a system under compression exceeds a certain value, it enters a regime known as a “jammed” state, where contacts are permanently established (3, 68–71). At jamming, where the system transits from liquid-like to solid-like state, the system becomes mechanically stable with finite bulk and shear moduli. Different parameters influence the coordination number. Higher friction lowers the coordination number (${\overline{C}}_J\approx$ 5.5 for friction equal to one; dashed line in Fig. 4), while lowering friction increases both the coordination number (${\overline{C}}_J\approx$ 6 for frictionless particles) and the volume fraction of jamming (72, 73).

One key observation one can make from comparing the coordination numbers of subnetworks of samples with rubber fractions ν = 0.1, 0.2, and 0.3 is whether each of the subnetworks is in a jammed state or not. The coordination numbers of the glass–glass subnetworks of samples with rubber fractions ν = 0.1 and 0.2 are above 5.5, confirming that these networks are surely jammed. This indicates that the glass subnetworks of ν = 0.1 and 0.2 are in a stable state with a finite stiffness not requiring rubber particles. However, the coordination number of the glass–glass subnetwork of the sample with a rubber fraction of ν = 0.3 is close to five, slightly below the jammed state coordination. Therefore, to reach a mechanical equilibrium with a finite stiffness, the glass–glass subnetwork relies on the rubber particles for ν ≳ 0.3.

Incorporating rubber particles—although it increases the coordination number and eventually offers a jammed packing—the overall stiffness of the packing starts to decay due to the softness of rubber. This explains the nonlinear, nonmonotonic behavior of the M-modulus observed in Fig. 1 for low rubber content samples. The presence of more and more rubber particles for ν = 0.1 and 0.2 makes the samples more deformable than an only glass sample, ν = 0, and the sample with ν = 0.2 is also more deformed in height than the sample with ν = 0.1, as we can see from SI Appendix, Fig. S2A.

Moving along the x-axis of Fig. 1, a further decline in the stiffness occurs at rubber fractions greater than 0.3. This decrease is associated with the change of contact networks from a stiff- to soft-dominated network (40). The regime 0.3 ≲ ν ≲ 0.5 in Fig. 4 is where the phase change starts since both glass and rubber networks carry almost the same coordination number. For rubber fractions greater than or equal to 0.5, the overall network behavior, such as stiffness, is fully controlled by soft particles, as their networks have gained a higher coordination number than that of stiff networks.

Previous observations (74) indicated that acoustical waves traveling through the granular medium are transmitted by strong force chains. Thanks to µXRCT imaging (61, 75–77), we are able to extract detailed information about the force chains of mixtures to explain the percolation chains of granular mixtures by postulating that the wave travels along the stiffest (fastest) path.

As earlier studies have shown, see refs. 39 and 40 and Fig. 1, increasing pressure leads to an increase in wave velocity and modulus. When a soft particle comes in between hard particles, due to its softness and viscoelastic features, it dissipates more energy, which leads to slower wave travel. If two identified chains have a similar length, the one which is only made of stiff particles is preferred by the wave to travel faster.

From top to bottom of the samples, a number of chains can be identified. Among possible paths for every sample, depending on selection criteria, the shortest length (the most linear-like) is selected as the fastest propagating path. Fig. 5 demonstrates the packings including one of the identified short chains, where red and blue particles represent glass and rubber particles. To identify the shortest possible chains, the following steps are taken: i) The first layer of particles on top of the samples is identified; ii) for each individual particle in this layer, the neighboring particles are determined. Two particles are considered neighbors if the distance between their centers is equal or less than the diameter of the particles (4 mm); iii) among the neighboring particles, the one with the smallest vertical coordinate (z-direction) is selected; iv) steps (ii) and (iii) are repeated until the last particle of the chain is found at the bottom of the sample. To have a comparison among the identified chains, their lengths are divided/scaled by the sample heights (almost identical in all cases, L_z⁰ = 80 mm). Table 1 reports the normalized average length of the first ten identified short chains ($\overline{\mathcal{L}}=\mathcal{L}/L_z$), their SD, number of particles involved in the chains ($\overline{\mathcal{N}}$), and the shortest path among all the ten identified chains (ℒ). As mentioned earlier, we have repeated µXRCT two more times for samples prepared with ν = 0.1, 0.2, 0.3 with new configurations by removing and adding the particles for every test. The statistical characteristics of these tests are also included in Table 1 for completeness.

From Fig. 4, we learned that the sudden drop of modulus for samples with ν ≳ 0.4 is due to a change of networks. The chains contain both glass and rubber particles, whereas in the case of lower rubber content, ν ≲ 0.3, many chains with only stiff particles are being identified. Although the chain length $\overline{\mathcal{L}}$ remains approximately the same among the identified chains, the chains of ν ≳ 0.4 are much softer, which leads to a dramatic change in the wave speed (modulus).

Comparing the length of the chains (glass only) at lower rubber content (ν ≲ 0.3) for three different realizations confirmed that the particles’ arrangements for samples with ν ≈ 0.1 and 0.2 provide shorter (i.e., straight/stiff) chains. This can be correlated with the fact that the glass subnetworks of these two samples are in a jammed state, therefore stable. The subnetwork of the sample with ν = 0.3 is (possibly) under jamming, which makes the system less stable with weaker force chains (78). Most of the glass networks for ν ≈ 0.1 and 0.2 offer relatively short chains in comparison to ν ≈ 0.3; therefore, it is not surprising to obtain a higher modulus.

Low-frequency ultrasound propagation using broadband piezoelectric transducers is a nondestructive way to investigate the small-strain (compression) stiffness of granular materials (14, 83). In the present contribution, cylindrical samples consisting of monodisperse soft (rubber) and stiff (glass) particles with diameters of 4 mm are poured into a custom-made X-ray transparent oedometer cell made of polymethyl methacrylate (PMMA). Glass and rubber beads used in this study are similar to those employed in an earlier study (30, 40), with Young’s moduli of 65 GPa and 0.185 GPa, and Poisson’s ratios of 0.2 and 0.5 for glass and rubber, respectively. The cell holds cylindrical samples with a length L_z⁰ = 80 mm and a diameter D = 80 mm with initial volume $V^0=\frac{\pi }{4}{L}_z^{\circ }{D}^2$ that are then placed under uniaxial static compression. A pair of P-wave broadband piezoelectric transducers with a frequency of 100 kHz (Olympus-Panametrics Videoscan V1011) are mounted on top and bottom of the cell and connected to an ultrasonic square wave pulser/receiver unit (Olympus-Panametrics 5077PR). The amplified signals were recorded with a resolution of 15 bit and a sampling rate of 125 Ms/s using a digital oscilloscope (PicoScope 5444B). To improve the signal-to-noise ratio, 32 signals were averaged. Samples composed of various rubber fractions ranging from ν = 0 to 1 in 0.1 increments were prepared, where the total volume of the cylinder is occupied by glass particles when ν= 0 and by rubber particles when ν = 1. The samples were prestressed in a gradual manner using static mechanical force. At different quasistatic axial stresses, p = 40, 80, 120, 160, 200 kPa, acoustic P-waves are propagated from the top transducer (source) into the prepared samples, and the bottom transducer (receiver) collects the signals. To ensure reproducibility, each mixture was repeated five times, with particles being removed and added to the cell each time to obtain a new particle arrangement and eliminate configuration dependency.

µXRCT imaging (84, 85) is arguably the most promising approach to attain the 3D microstructure of particle packings in high resolution by reconstructing the internal structure from a high number of 2D projection images (radiograms) acquired from different directions. The main advantage of µXRCT imaging in comparison to other nondestructive imaging methods (e.g., magnetic resonance imaging) is its ability to resolve the internal microstructure with high accuracy (86, 87). It is generally possible to visualize features down to the single-digit micrometer range (84, 85). Thus, it has been widely used to attain microscopic features of particulate systems. To further enhance a micromechanical insight into the mixture packings, the phase transition from stiff- to soft-dominated (0.1 ≤ ν ≤ 0.6), and the stiffening behavior of the packings with a low number of soft inclusions (0.1 ≤ ν ≤ 0.2), at stresses of 80 kPa and 160 kPa, approximately, were visualized by in situ µXRCT. For this, a lab-based modular buildup µXRCT system further described in ref. 88 was used. The underlying used experimental setup is illustrated in Fig. 6A showing the combination of conventional imaging setup and wave propagation measurements into an oedometer cell under uniaxial preload. It is worth mentioning that due to the existence of rubber particles, it is not possible to apply diffraction measurements of contact fabric, i.e., interparticle force, on the samples.

Unlike wave propagation tests, µXRCT scans are expensive both in energy consumption and time. Thus, we performed µXRCT scans on two vertical stresses p= 80 kPa and 160 kPa. Nevertheless, three different configurations (networks) were prepared for samples with rubber fractions ν = 0.1, 0.2, and 0.3, whereas for samples with rubber fractions ν = 0.4, 0.5, and 0.6, only one test was performed. All 24 (two vertical stresses × six rubber contents, 3× ν : {0.1, 0.2, 03} + 1× ν : {0.4, 0.5, 06}) scans were performed with identical image acquisition settings. To capture the whole cell content (D = 80 mm; L_z⁰ = 80 mm), a geometric magnification of 1.36 was used, given by the ratio of source–detector distance and source–object distance. Using 2 × 2 detector binning, the final resulting voxel edge length is 149.6 µm and sufficiently high to resolve the beads having a 4 mm diameter. The X-ray tube voltage and flux were set to 110 kV and 110 µA, and the detector exposure time to 1,000 ms. For the scan, the loaded sample was rotated 180° clockwise and counterclockwise in 0.25° increments to cover the necessary complete turn. To enhance the final image quality with regard to ring artifacts, two slightly shifted projections for each angle were acquired and combined. Further information on the experimental setup and the used parameters can be found in ref. 89. For each scan, stress relaxation in the order of 10 kPa occurred quickly after switching from force- to a displacement-controlled mode for the imaging. However, we observed a quick recovery of the load upon reloading. Each scan took about 120 min and was performed displacement-controlled to avoid potential movement errors. The 3D volume reconstruction based on the radiograms of each scan was done with the FDK algorithm (90) using the commercial software Octopus Reconstruction (v.8.9.4-64 bit) (91). Image processing of the reconstructed gray value image stacks was performed with Dragonfly software (v.2020.2) (92); see Fig. 6B. First, a mathematical smoothing function (median filter) was applied to reduce unwanted noise in the tomograms. Accurate segmentation of images was achieved due to the high difference in attenuation coefficients between glass, rubber, and air. Finally, glass and rubber particles were identified using a mathematical morphology technique, watershed transform, on defined markers of individual particles to determine grain volumes, grain centroids, and all contact locations and orientations (93–95). To identify contacts between particles, immediate neighbor voxels containing two distinct grain IDs were considered (96, 97). To select a so-called true contact, a minimum of 10 voxels between two particles was set as a threshold for a contact candidate (87, 95).

For further analysis, we demonstrate the samples only at p = 160 kPa since the topological response of p = 80 kPa is almost unchanged. Results of the analysis of p = 80 kPa are reported in supplementary information file. The in situ acquired data of all scanned samples, including metadata, are open-access published (98–101); see SI Appendix for more details.

K.T., M.R., and H.S. designed research; K.T. and M.R. performed research; K.T. contributed new reagents/analytic tools; K.T. analyzed data; S.L. contributed to discussions; H.S. coordinated the research; S.L. and H.S. revised the paper; K.T. and M.R. wrote the paper.

The authors declare no competing interest.