Local mobility and microstructure in periodically sheared soft particle glasses and their connection to macroscopic rheology

Oscillatory shear rheology is widely used to characterize the rheological behavior of materials as diverse as colloidal suspensions, polymer melts and solutions, surfactants mesophases, liquid crystals, and soft glassy materials [Hyun et al. (2011)]. A typical dynamic oscillatory test is performed by subjecting the material to a periodic strain or stress of sinusoidal shape at angular frequency ω and measuring its mechanical response as a function of time. The linear viscoelastic (LVE) regime of the material is probed when small amplitude oscillatory strains are applied. The LVE properties are characterized by the (elastic) storage modulus $G′$ and (viscous) loss modulus $G″$⁠. Materials can be probed at a large range of frequencies or time scales, often inaccessible by steady shear experiments, which yield the characteristic relaxation times and moduli. $G′(ω)$ and $G″(ω)$ fully describe the material response near equilibrium and provide insights into the relation between microstructure and rheology. Predictions of moduli based on molecular or microscopic theories are now available for a wide range of materials.

At large amplitude oscillatory shear (LAOS), the equilibrium structure is distorted far from equilibrium and the response is no longer proportional to the excitation. The analysis used for small strain amplitudes is no longer meaningful and is insufficient to characterize the nonlinear rheology of materials. The complexity and richness of LAOS tests reside in two main features. First, the amplitude and the frequency can be varied independently. Second, LAOS tests are usually nonequilibrium processes in that the material periodically returns to the same state of deformation without necessarily having time to reach steady state during the cycle. This makes LAOS experiments useful for studying highly nonlinear situations where materials are subjected to large deformations for a short period of time, as in fatigue tests.

The analysis of LAOS measurements constitutes a very active area of research [Hyun et al. (2011)]. The oscillatory stress responses are often visualized in the form of closed curves in the strain and stress versus shear-rate planes, and are generally referred to as elastic and viscous Lissajous–Bowditch plots, respectively [Dealy and Wissbrun (1990); Philippoff (1966)], or as three-dimensional (3D) space with stress, strain and strain rate coordinates [Cho et al. (2005)]. General methods of analyzing nonlinearities consist of decomposing the response waveforms into various basis functions, yielding a series expansion with an arbitrary number of terms. Fourier-transform rheology uses trigonometric functions to represent the response into a series of harmonics in the frequency domain [Wilhelm et al. (1998); Wilhelm (2002)]. Cho et al. (2005) decompose the response into a sum of so-called elastic and viscous contributions. The idea was later generalized by Ewoldt et al. (2008) who used a set of orthogonal Chebyshev polynomials. Another method of decomposition uses combinations of characteristic basis functions, which are associated with specific physical events [Klein et al. (2007)].

Despite these important achievements, the interpretation of LAOS experiments in terms of physical mechanisms remains challenging. The question of how data extracted from LAOS experiments correlate to constitutive models has motivated the investigations of various materials, such as polymer solutions and melts with various molecular architectures [Pearson and Rochefort (1982); Hyun and Wilhem (2009); Wagner et al. (2011)], wormlike micelles [Gurnon and Wagner (2012)], and yield stress fluids [Ewoldt et al. (2010)]. A conceptually useful framework analyzes the response waveforms as temporal sequences of physical processes that repeat over time [Rogers et al. (2011a, 2011b); Rogers and Lettinga (2012); Rogers (2012)]. With this technique, it is possible to acquire multiple parameters that would otherwise require several independent rheological tests. In some instances, it has been possible to relate the LAOS response to the structural changes occurring during cycles [Rogers et al. (2012); Gurnon and Wagner (2012); López-Barrón et al. (2012)]. Nonetheless, fundamental descriptions where the nonlinear shear response is derived explicitly from microscopic theories are absent. Here, we determine the microstructural changes that occur during LAOS experiments for soft particle glasses with particle scale simulations and use that to establish the microscopic events that occur during the strain cycle. We then connect these microscopic events to the macroscopic behavior of soft particle glasses during oscillatory shear.

Soft particle glasses form a broad family of materials made of deformable particles, which are jammed at volume fractions exceeding close-packing of equivalent hard spheres (⁠$ϕ=0.64$ for monodisperse particles). These materials are as diverse as microgels [Cloitre et al. (2003)], emulsion droplets [Mason et al. (1995)], star polymers [Likos (2006)], and block copolymer micelles [Buitenhuis and Forster (1997)]. Their structure, which is amorphous, shares common features with hard sphere glasses. Each particle is trapped in a “cage” formed by its neighbors, which restrict and even arrest macroscopic motion. However, whereas hard sphere colloids only experience excluded volume interactions, soft particles at high volume fractions are compressed against each other by bulk osmotic forces and particles which are in contact interact via soft elastic repulsions. The sources of particle elasticity can be different for different materials [Bonnecaze and Cloitre (2010)], like in microgels where the particle elasticity is due to the electrostatic repulsions of counterions, in concentrated emulsions where it is the surface tension between the two phases and in star polymers where the elasticity is due to the foldable arms. The rheological properties of soft jammed materials are dominated by the existence of the so-called yield stress below which they behave like weak elastic solids whereas they yield and flow like viscous liquids above it.

In a previous work, we showed that steady shear properties result from a subtle interplay between elastic interactions and structural rearrangements between interlocked particles. We proposed a micromechanical model which provides physical insights into yielding and flow mechanisms and leads to quantitative predictions of the shear stress and the normal stress differences without any adjustable parameters [Seth et al. (2011)]. Here, we extend this approach to soft particles glasses subjected to oscillatory shear flow. In Sec. II, we describe the model, the simulation technique, and the experimental protocol used to conduct oscillatory experiments on well-characterized microgel suspensions. In Sec. III, we present the simulation results at small and large amplitudes for the dynamics at the particle scale, the microstructure at the mesoscopic scale, and the macroscopic rheology. The viscoelastic moduli at low strain amplitudes, the strain sweeps, and the Bowditch–Lissajous (BL) plots are found in good agreement with real experiments, which validates the technique. In the discussion of Sec. IV, we connect the evolution of the particle mobility and the microstructure to the macroscopic rheology and demonstrate that a cycle can be decomposed into a sequence of physical processes that repeat periodically. In Sec. V, we conclude and identify some interesting analogies between the behavior of soft particles glasses and other disordered suspensions.

We model soft particle glasses as 3D packings of N periodically replicated non-Brownian elastic spheres dispersed in a solvent with viscosity η at volume fractions exceeding the random close-packing of hard spheres [Lacasse et al. (1996); Seth et al. (2006)] [Fig. 1(a)]. The particles have a Gaussian size distribution with an average radius R and a standard deviation of 10% in order to avoid crystallization. The 3D packings are built as follows. We first generate periodically replicated random close packed configurations of hard spheres in a simulation box using the compression algorithm introduced by Lubachevsky and Stillinger (1990). After forming the close-packed microstructure, the spheres are treated as deformable particles and compressed by decreasing the box size in small decrements until the desired concentration is achieved. After each decrement, the system is allowed to relax using a conjugate gradient algorithm to minimize the system energy and ensure the net contact forces and torques on each particle vanishes. The particles form flat facets at contact resulting in the relative deformation $εαβ=(Rα+Rβ−rαβ)/R$ where $rαβ$ is the center-to-center distance between particles α and β, and $Rc=RαRβ/(Rα+Rβ)$ is the contact radius [Fig. 1(c)]. The suspension is subjected in the x-direction to an oscillatory shear strain of amplitude and frequency, $γ0$ and $ω$⁠, respectively: $γ=γ0sin(ωt)$⁠. The resulting velocity is in the x-direction with y axis being the gradient direction.

To study the dynamics of the suspension, we extend the micromechanical model initially proposed by Seth et al. (2011) for steady shear flow. The elastic repulsion between two particles $α$ and $β$ is modeled with a modified Hertz potential

In the above expression, $E*$ is the particle contact modulus $E*=E/2(1−ν2)$ (⁠$E$⁠: Young's modulus, ν: Poisson ratio), and $n⊥$ is the direction perpendicular to the particle–particle facet at contact. Parameters n and C depend on the level of compression between particles: n = 1.5 and C = 1 for ε_αβ < 0.1; n = 3 and C = 32 for 0.1 ≤ ε_αβ < 0.2; n = 5 and C = 790 for 0.2 ≤ ε_αβ < 0.6. This form of interaction was initially derived for elastomeric spherical particles subjected to large compression [Liu et al. (1998)].

When the suspension is subjected to flow, there is a sheared, thin film of solvent at the facets between the particles which creates an elastohydrodynamic (EHD) drag force of the form [Seth et al. (2011)]

where $uαβ,//$ and n_// are the relative velocity and the direction parallel to the flat facet developed at particle–particle contact, respectively.

Each particle experiences a net elastic repulsion and an EHD drag due to all the particles in contact with it and these forces are assumed to be pairwise additive. The particles also experience an effective drag force due to relative motion with the solvent. Neglecting particle and fluid inertia, the sum of all forces on each particle is zero. The resulting equation of motion can be made dimensionless by scaling lengths, time and velocity by R, $η/E*$ and $RE*/η$⁠, respectively. It has the form

x_α is the particle position; $fr(ϕ)$ is a coefficient that accounts for the hindered mobility of particles due to the high concentration; it was equal to 0.01 for these simulations, but the results did not vary significantly for values of $fr(ϕ)$ less than 0.05; $uα∞=γ̇ηE*yex$ represents the velocity field due applied strain, $γ̇=γ0ωcos(ωt)$ being the instantaneous shear rate and e_x the unit vector in the x-direction.

These N coupled equations of motion were integrated numerically to determine the evolution of the spatial position and velocity of each particle. Periodic boundary conditions were applied in the x- and z-directions and Lee–Edwards boundary conditions were implemented in y-direction. At each step, the positions of particles leaving the box were shifted by a cumulative strain $γ=∫0tγ̇dτ$⁠. The open source code LAMMPS [Plimpton (1995)] was used to perform the simulations. Packings of 1000 and 10 000 particles were generated. Configurations with 1000 particles were used for the large amplitude simulations, $γ0$ = 2.0, 1.0, 0.5, and 0.1. Each simulation was run up to 50 cycles and averaged over five different initial configurations. Configurations with 10 000 particles were used for the small amplitude simulations, $γ0$ = 0.05, 0.01, and 0.003. These simulations were run for up to 20 cycles and averaged over five different initial particle configurations. The cycles for averaging were chosen after reaching steady state. The approach to steady state was defined as the cycle beyond which all other cycles were the same, i.e., the stress–strain data for all cycles overlapped on top of each other. The average stress at different points during the cycle was computed from the Kirkwood formula $σ=−1V∑α∑β>α(fαβEHD+fαβe)(xα−xβ)$⁠, where V is the box volume [Larson (1997)]. The viscoelastic moduli were computed by taking the Fourier transform of the stress–strain data computed over a cycle.

The microstructure of suspensions was characterized through the dynamic pair distribution function g(r) [Brady and Morris (1997); Morris and Katyal (2002); Sierou and Brady (2002); Seth et al. (2011)]. The dynamic pair distribution function g(r) was computed at different points in the cycle. It can be decomposed into a series of orthogonal spherical harmonic functions $g(r)=gs(r)+∑l=1∞∑m=−llglm(r)Ylm(θ,ϕ)$ [Hanley et al. (1987)].

The suspensions used in this study are prepared from polyelectrolyte microgels in water. The microgels were synthesized by standard emulsion polymerization at low pH (≅2) using the two monomers ethyl acrylate (64 wt. %) and methacrylic acid (35 wt. %), and dicyclopentenyloxyethyl methacrylate (1 wt. %) as a crosslinker. In order to remove unreacted monomers, surfactants and other impurities present in the final emulsion, the polymer latexes obtained from the synthesis were cleaned by ultrafiltration. The solid content of the stock suspension was determined by thermogravimetry and samples were subsequently prepared by dilution with ultrapure water (η_s = 9.8 × 10⁻⁴ Pa s). At low pH, the microgels are insoluble in water, and hence they are essentially hard particles. When sodium hydroxide (1M) is added to the microgel suspensions, the acidic units become ionized and the osmotic pressure of the counterions provokes the swelling of the microgels. In a previous study we have shown that most of the counterions are trapped in the polymer network, so that the net charge carried by the particle is extremely low and electrostatic interactions are negligible [Cloitre et al. (2003)]. In dilute suspensions, the swollen particles have a spherical shape with a hydrodynamic radius R (≅295 nm). Above C_m = 1.5 × 10⁻² g/g, the suspensions exhibit solid-like properties with a yield stress. The concentration can be increased much further because microgels are able to deform and deswell osmotically. The concentration of the suspension investigated in the following is C = 2 × 10⁻² g/g.

Rheological measurements were carried out using an Anton Paar MCR 501 rheometer mounted with a cone and Peltier plate geometry with a diameter of 50 mm, a 2° angle, and a truncation of 48 μm. The shearing surfaces were sandblasted to provide a surface roughness of 2–4 μm, which prevented the occurrence of slip. A solvent trap was placed around the sample to minimize water evaporation, the interior atmosphere of the trap being saturated using a few droplets of distilled water. All measurements were made at 20.0 ± 0.1 °C. Prior to any measurement, the suspensions were presheared at a shear rate 500 s⁻¹ for about 30 s in order to erase their mechanical history. Then they were kept at rest for a waiting time of about 4 h, which is sufficiently long to make aging effects negligible.

We performed several types of rheological tests. Steady shear experiments, where the stress σ was measured as a function of the applied shear rate $γ̇$⁠, were performed by applying constant shear rates varying from 10³ to 10⁻⁴ s⁻¹ and recording the stress until steady state was reached. Oscillatory frequency sweeps were used to measure the storage modulus G′ and loss modulus G″ as function of the angular frequency ω (10⁻² < ω < 10² rad/s) at small strain amplitudes in the LVE regime (γ₀ = 6 × 10⁻³). The storage and loss moduli show the characteristic variations exhibited by many soft materials, i.e., a nearly constant plateau in G′(ω) and a much lower G″(ω) with a small minimum around a frequency ω_m. We define the elastic modulus of the suspension G₀ as the value of G′(ω) at ω_m. LAOS measurements were performed in strain-controlled mode using the expert mode package provided by the Anton Paar rheometer software. The strain and stress signals were recorded at a sampling rate of 256 data points per oscillation cycle. The strain amplitude was varied from 3 × 10⁻³ to 10 and the frequency was in the range between 0.3 and 15 rad/s. The results are presented in the form of oscillatory strain sweep plots and Lissajous–Bowditch curves.

In order to compare simulations and experimental results, we need to characterize the material properties of our suspensions. The storage modulus G₀ was determined from oscillatory frequency sweep tests as explained in Sec. ???. The yield stress σ_y was obtained from the variations of the stress in strain amplitude sweep tests at low frequencies. The yield strain γ_y is obtained from the relation σ_y = G₀ γ_y. An important parameter involved in the micromechanical model presented above is the particle contact modulus $E*$⁠, which is unknown a priori and cannot be easily measured for submicron particles. To circumvent this difficulty, we used the following method which was first proposed by Seth et al. (2011). We performed steady shear experiments to measure the flow curve $σ(γ̇)$ of the suspension. It reaches a constant value at low shear rates, which corresponds to the yield stress measured in oscillatory tests. To determine $E*$⁠, we fitted the experimental flow curve to the theoretical expression predicted by Seth et al. (2011) 

where k ≅ 80 is a numerical coefficient deduced from simulations. Since σ_y, γ_y, and η are known independently, $E*$ is the only free parameter. The experimental and fitted flow curves are shown in Fig. 2 and the resulting value of $E*$ is reported in Table I. In the inset of Fig. 2, by comparing the experimental value of G₀/$E*$ to the theoretical variations of G₀/$E*$ with volume fraction [Seth et al. (2006); Seth et al. (2011)], we estimate the effective volume fraction ϕ, which is of the order of 0.80. This justifies that the simulations reported in the flowing were performed at ϕ = 0.80. The values of the elastic modulus, yield strain, and yield stresses measured experimentally and obtained from simulations are summarized in Table I.

In this section, we study particle dynamics in 3D jammed packings of soft particles subjected to oscillatory shear flow. We track individual particle motions in the simulations and we analyze the mean square displacements of the particles.

Figure 3 represents the variations of the mean square displacements in the x-, y-, and z-directions versus the number of oscillations n_osc (=ωt/2π) for a fixed frequency and various strain amplitudes. In the x-direction, the mean square displacement is modulated by oscillations due to the periodic forcing unlike along y- and z-directions. The top panel shows data computed at low strain amplitudes γ₀ close or smaller than the macroscopic yield strain, γ_y; in the bottom panel, γ₀ is much larger than γ_y. At the smallest strain amplitudes investigated (γ₀/γ_y ≤ 1.5), the mean square displacements 〈Δx²〉, 〈Δy²〉, and 〈Δz²〉 slowly increase in time before reaching a plateau value which is much less than a particle diameter. This indicates that the particles explore restricted regions of space around their mean positions. In other words, they can be seen as trapped in cages formed by their neighbors. At larger strain amplitudes (γ₀/γ_y ≥ 3), the particles move over distances larger than their radius, i.e., they escape their local environments and jump to other cages where they experience different neighbors. These rearrangements give rise to large mean square displacements of individual particles even though the macroscopic strain is purely oscillatory.

Figure 4 presents the mean square displacements in the x-, y-, and z-directions at small and large strain amplitudes when the frequency is varied. In the top panel, the strain amplitude is smaller than the yield strain (γ₀/γ_y = 0.09) and the nondimensional frequency ηω/$E*$ varies between 2 × 10⁻⁸ and 10⁻³. In this low strain amplitude regime, the particles get trapped at all frequencies after a few oscillations. The steady state, where the mean square displacements reach their plateau values, occurs at longer times for lower frequencies. It is interesting to note that the plateau values of 〈Δx²〉, 〈Δy²〉, and 〈Δz²〉 are similar, as already noted in Fig. 3. The bottom panel shows data computed for a large strain amplitude (γ₀/γ_y = 30) and nondimensional frequencies varying again between 2 × 10⁻⁸ and 10⁻⁶. The mean square displacements steadily increase at all frequencies showing that the particles are not permanently trapped at these large amplitudes. It is interesting to note that the mean square displacement per cycle remains of the same order even though the time taken for completing one oscillation is orders of magnitude different for the different frequencies.

For strain amplitudes larger than the yield strain, the particle mean square displacements 〈Δx²〉, 〈Δy²〉, and 〈Δz²〉 increase with time. After steady state is reached, the variations are linear with periodic modulations in the oscillatory flow direction (Figs. 3 and 4). This allows us to define effective diffusion coefficients of the particles in the x-, y-, and z-directions (the oscillatory shear is in the x-direction) by fitting the mean square displacement curves as follows:

In practice, the mean square displacements were computed over the four last cycles. In Fig. 5(a), we present the variations of the nondimensional effective diffusion coefficients D_x, D_y, and D_z normalized by the scaling factor D₀ = R²Ε^∗/η, for different strain amplitudes at $ηω/E*=2×10−8$⁠. Again, the amplitude of the maximum strain γ₀ is normalized by γ_y. The form of D₀ expresses the fact that the only natural time and length scales are the characteristic time η/$E*$ and the particle radius. All three diffusion coefficients are equal, except at the largest strain amplitude where the diffusion coefficient in the x-direction is slightly larger. This indicates that the shear-induced diffusivity associated with particle rearrangements is isotropic under these conditions, even though the shearing motion is applied along the x axis. In the following, we will only consider the averaged diffusion coefficient, D = (D_x + D_y + D_z)/3. The diffusion coefficients increase linearly with strain amplitude beyond γ₀/γ_y ∼ 1. This further confirms that the large scale rearrangements responsible for the increase of the shear-induced diffusivity start taking over around the yield strain.

We also studied the effect of frequency on the particle scale diffusivity at large amplitudes. In Fig. 5(b), we plot the variations of the averaged diffusion coefficient D with the nondimensional frequency ηω/$E*$ for two values of γ₀/γ_y > 1. D increases linearly with the applied frequency. In Fig. 5(c), the averaged diffusion coefficients for different strain amplitudes greater than the yield strain and frequencies are plotted against the nondimensional shear-rate amplitude $γ̇0η/E*$ where $γ̇0=γ0ω$⁠. The data obtained for different strain amplitudes and frequencies collapse onto a single master curve. They are well described by a linear variation. It is interesting to note the similarities between the nondimensional variable $γ̇0η/E*$ and that involved in the constitutive equation (4). This suggests that the same underlying mechanisms, i.e., disorder and EHD interactions lay at the heart of the behavior of soft particle glasses at the local and macroscopic scale.

From Fig. 5(c), one may also conclude that the diffusivity $D≅0.1γ̇R2$⁠. Remarkably, this is practically the same scaling for shear-induced hydrodynamic diffusion for hard spheres in suspensions with volume fractions greater than 30% [Nott & Brady (1994)], despite the fact that the physics driving the behavior (long-ranged hydrodynamic interactions versus near-field elastohydrodynamics) are quite different.

Let us first revisit the microstructure of the suspensions at rest following Seth et al. (2006). Figure 6(a) shows the radial pair distribution function computed at rest at a volume fraction $ϕ$ = 0.8. The constituent particles are soft, and therefore at concentrations greater than the hard sphere random close-packing limit (⁠$ϕ$ = 0.64), they are pressed against one another and the radial separation at which the pair probability is maximum is less than twice the particle radius. This contrasts with the case of hard spheres where the radial pair distribution function shows a sharp rise at the hard sphere diameter. Figure 6(b) shows the two-dimensional pair distribution function at rest in the x–y plane. Since the range of compression between particles is small compared to the particle radius, it is beneficial to represent the microstructure in the r–$θ$ or azimuthal plane, so that we can focus on the range of separations where particle contacts occur and study the angular distortion more closely. In Fig. 6(c), we clearly see that the microstructure at rest is radially symmetric without any preferential orientation or angular distortion. As a reference, the particle separation where the pair distribution function is maximum is indicated by a white line. In the following, we will use the azimuthal representation to present the two-dimensional pair distribution functions computed within the cycle at small, medium, and large strain amplitudes.

We first study the evolution of the microstructure within a cycle when the strain amplitude γ₀ is smaller than the yield strain (γ₀/γ_y < 1). Figure 7(a) shows the variations of the applied strain γ and of the resulting stress σ. We observe that γ and σ are nearly in phase and linearly proportional, which corresponds to the LVE regime. We have computed the pair correlation function at five particular locations of the cycle, denoted by t_i (i = 1–5). t₁ and t₅ are the points where the stress amplitude has positive and negative extrema, respectively; t₃ is where the stress is zero; at t₂ and t₄, the stress takes intermediate values of opposite signs. Note that in this small strain limit where γ and σ are nearly in phase, the strain has also its positive and negative extrema approximately at t₁ and t₅ and is zero in t₃. The five pair distribution functions which are computed are shown in Fig. 7(b). At all times, the pair distribution function appears to be essentially uniform. The particle separation r/R at which the pair correlation function is maximum is about the same as at rest. To better analyze the two-dimensional pair correlation function, we have calculated the spherical harmonics g_2,−2(r) at t₁, t₃, and t₅, which are plotted in Fig. 7(c). At t₃ where the stress is zero, g_2,−2(r) vanishes within the statistical fluctuations. At t₁ and t₅ where the stress has its extrema, g_2,−2(r) exhibits a small negative minimum at small center-to-center distances r/R and a small positive maximum at large distances during the positive part of the cycle (t₁) and vice versa during the negative part (t₅). This means that there is a small accumulation of particles center along the compression axis or in the upstream quadrant and depletion along the extension axis or in the downstream quadrant. In conclusion, a small amplitude oscillatory shear flow induces only very slight distortions of the pair distribution function with respect to the static situation.

Figure 8 analyzes the change of microstructure when the strain amplitude is slightly larger than the yield strain (γ_0/γ_y ∼ 1), the maximum stress being of the order of the yield stress. We refer to this situation as the medium amplitude case. Figure 8(a) shows the variations of the applied strain and resulting stress within a cycle. The stress waveform is no longer sinusoidal indicating that we are outside the LVE regime. We have computed the two-dimensional pair correlation function at six characteristic intracycle positions: t₁ and t₄ where the stress has its positive and negative extrema, respectively; t₃ and t₆ where it is zero; and t₂ and t₅ where the strain has its positive and negative extrema. The corresponding azimuthal plots are presented in Fig. 8(b). We observe that the pair distribution functions are no longer radially symmetric and reveal significant distortions of the microstructure. These distortions are the most visible at the positions in the cycle where the stress is the largest (t₁ and t₄, t₂ and t₅). The positions of the maxima of the particle pair distribution functions vary with the angle and the distribution of particles centers along the radial direction is no longer symmetric. In particular, there are particles at smaller separations than found at rest or in the LVE regime. This indicates that neighboring particles tend to accumulate in the upstream quadrant (π/2 < θ < π), where they are more compressed, and to deplete along the extensional axis (0 < θ < π/2), where they are less distorted. It is interesting to note that the angular position of the accumulation and depletion effects in the negative part of the cycle is shifted by π/2 radians with respect to their positions in the positive part due the reversal of the flow. At t₃ and t₆, where the stress vanishes, the accumulation–depletion effect is significantly reduced, the particle distribution becomes more uniform and resembles to that in the static state.

The accumulation–depletion mechanism that we have identified is also supported by the variations of the spherical harmonics shown in Fig. 8(c). g_2,−2(r) has a negative peak at low separations (r/R ≅ 1.83), i.e., large degrees of compression, which corresponds to an accumulation of particles along the compression axis. There is a positive peak of smaller amplitude at large separations (r/R ≅ 1.9), i.e., low degrees of compression, which corresponds to the depletion of particles around the downstream extension axis of the flow. From the depth of the minimum and the height of the maximum, we deduce that they are fewer particles in the extensional region. The asymmetry between the peaks associated with upstream accumulation and downstream depletion is more pronounced when the stress is large (t₁ and t₄, t₂ and t₅). On the contrary, at t₃ and t₆, where the stress vanishes, the asymmetry is weak and g_2,−2(r) approaches that computed for the low amplitude case. We also note that the spherical harmonics g_2,−2(r) computed at locations where the stress takes opposite values (t₁ and t₂, t₃ and t₆) are equal but with opposite signs, which again express the fact that the particle distribution simply reverses upon flow reversal.

Figure 9 shows the evolution of the microstructure during strain amplitudes that are large relative to the yield strain (γ₀/γ_y ≫ 1); here, the maximum stress exceeds the yield stress. In Fig. 9(a) the stress waveform is considerably distorted with respect to the strain, indicating that the response is highly nonlinear. In Fig. 9(b), we show the pair distribution functions for the same characteristic times as in the medium amplitude case. The evolution of the microstructure within the cycles is somewhat similar to that for the medium amplitude case but the amplitudes of the distortions are much larger. At times t₁ and t₄, where the stress reaches its maximum positive and negative values, which largely exceeds the yield stress, we observe important asymmetries in the pair distribution function both along the azimuthal and radial directions. There is a significant accumulation of particles upstream and depletion downstream indicating an important redistribution of particles. At t₂ and t₅, where the strain reaches its positive and negative values, the stress has significantly decreased and the accumulation and depletion regions are much less pronounced. Finally, at t₃ and t₆, where the stress goes to zero, the particle distribution is nearly isotropic, the maximum compression and the average separation between the particles being comparable to that observed at rest or at small strain amplitudes.

The spherical harmonics g_2,−2(r) computed at the characteristic times t_i are presented in Fig. 9(c). As already discussed for the medium amplitude case, the negative minimum and positive maximum during the positive part of the strain cycle are associated to the accumulation of particles in the compressive quadrant and the depletion in the extension quadrant. The peak heights are the largest at t₁ and t₄, when the stress is the largest. They decrease rapidly at t₂ and t₅, when the stress amplitude drops and finally they nearly vanish at t₃ and t₆, when the stress goes to zero. Again we observe symmetry between upstream and downstream upon flow reversal.

Oscillatory experiments at low strain amplitudes probe the LVE response of the unperturbed microstructure. We have performed systematic simulations in this regime, varying the frequency of oscillations over several decades, and compared the results with experimental observations. The variations of the storage and loss moduli with frequency are shown in Fig. 10. At low frequencies, we probe the elasticity and dissipation associated with the deformation of the particles through their contacting facets. The storage modulus exhibits a plateau that is much larger than the loss modulus. The value of the plateau modulus, G₀, is found in good agreement with independent calculations based on uniaxial deformation of small amplitude [Seth et al. (2006)]. For convenience, all the storage and loss moduli presented are scaled by the plateau modulus G₀. At high frequencies, both the storage and the loss moduli increase. For reference, we show the power law variation which is often followed by the loss modulus in experiments on compressed emulsions and foams (⁠$G″$ ∼ ω^1/2) [Liu et al. (1996); Cohen-Addad et al. (1998)]. Although more simulations at higher frequencies would be necessary to draw a definite conclusion, the data calculated in our simulations fit reasonably well to this expectation. In Fig. 10, we have also plotted the storage and loss moduli which were measured for the microgel suspension described in Sec. II B, using a conventional rheometer. For the sake of comparison, the experimental moduli are normalized to the plateau modulus measured experimentally and the reduced frequency is calculated using the value of the contact modulus which was determined experimentally and the solvent viscosity. The correspondence between the simulation and the experimental data validates the ability of the 3D EHD model to reproduce quantitatively the LVE behavior of real soft particle glasses.

Similarly, we can calculate the effective storage and loss moduli as well as the stress amplitude for any strain of arbitrary amplitude and construct simulated strain sweep plots. The results are presented in Fig. 11 together with the corresponding experimental data measured for the $ϕ=0.8$ microgel suspension. The moduli are again normalized by the plateau modulus G₀, the stress by the yield stress, and the strain amplitude by the yield strain. Both sets of data agree quantitatively, which further confirms the capacity of the 3D model and our simulations to reproduce nonlinear responses in oscillatory shear rheology. Figure 10 gives evidence for a change of mechanical behavior around γ₀/γ_y ≅ 1, where the maximum stress exceeds the yield stress. The small strain amplitude regime where γ₀/γ_y < 1, intermediate strain amplitude regime where $γ0/γy≈1$⁠, and large strain amplitude regime where γ₀/γ_y > 1 can be mapped on the regimes defined in the analysis of the microstructure in Sec. III B. In the low strain amplitude regime, where the microstructure is essentially nondisturbed by the mechanical excitation, both moduli are constant with $G′>G″$ indicating linear elastic behavior. The intermediate strain amplitude regime, where large scale rearrangements start to occur, is associated with a break in the stress curve, a decrease of the storage modulus and a bump in the loss modulus. At large amplitudes, where particles continuously rearrange and move over long distances, both $G′$ and $G″$ vary as power laws according to $G′∝γ0-μ$ and $G″∝γ0-ν$ with μ ≅ 1.45 and ν ≅ 0.80 (μ/ν ≅ 1.8). Similar variations are ubiquitous in soft glassy materials [Hyun et al. (2002); Miyazaki et al. (2006); Erwin et al. (2010)]. One major drawback of this analysis, however, is that these effective moduli which are calculated or measured experimentally in the nonlinear regime provide only qualitative information about the actual rheological properties of the materials. In Sec. ???, we go further by presenting a full nonlinear analysis of the rheological response using Bowditch-Lissajous (BL) plots.

Stress–strain plots or elastic BL for different strain amplitudes are shown in Fig. 12. Figures 12(a)–12(c) present results from simulation. In Fig. 12(a), the strain amplitude is smaller than the yield strain (γ₀/γ_y ≅ 0.09). The elastic BL plot is an ellipse, which is the typical response expected for an LVE material. The slope of the long axis is the storage modulus while the short axis, which characterizes the openness of the ellipse, depends on the loss modulus. The elongated shape of the ellipse indicates that $G′>G″$⁠, which is found in agreement with the previous results shown in Fig. 10. Ellipses with similar features are obtained for different strain amplitudes as long as γ₀/γ_y < 1, since we remain in the LVE regime. In Fig. 12(b), we present data at intermediate strain amplitudes. The BL curve for γ₀/γ_y = 1 is still an ellipse but that for $γ0/γy=3.0$ it begins to deform. The distortions are the largest when the strain becomes of the order of the yield strain and the stress comparable to the yield stress. The BL plot is then stretched along the strain axis. Figure 12(c) presents data for large strain amplitudes. The BL plots then adopt a parallelogram shape. The lateral sides of the parallelogram correspond to the regime where the stress is smaller than the yield stress. The top and bottom horizontal sides of the BL curves correspond to the portion of the cycles where the stress is larger than the yield stress. It is interesting to note that the point where yielding starts is signaled by a small overshot of the stress, which is reminiscent of the static yield stress. For the sake of comparison, Figs. 12(d)–12(f) present experimental data for the concentrated microgel suspensions. They agree reasonably well with the simulated data both in shape and in amplitude. There is a small difference near the yield point around σ = σ_y, which is less apparent for the microgel suspension. Note that a well-defined static yield stress has been observed star polymer glasses [Rogers et al. (2011a)]. We have also investigated the importance of frequency both in simulations and in experiments. The BL plots essentially keep the same characteristic shapes in the low, intermediate, and large amplitude regimes. In the large amplitude regime, the slope of the lateral sides of the BL plots increases with frequency, in accordance with the increase of the linear storage modulus; the absolute value of the stress in the top and bottom sections of the BL plots, which corresponds to macroscopic flow, increases.

The particle scale simulations and experiments described in Sec. III provide a deep understanding of the microstructural changes occurring during the oscillatory shear of soft particle glasses. At small strain amplitudes, the particles remain trapped in the cages formed by their neighbors as they periodically move back and forth to follow the applied strain. Interestingly, the particles execute some in-cage motion of small amplitude around the mean position. The mean square displacement of the particles slowly increases in time and tends to a constant plateau after many cycles. The spherical harmonic g_2,−2(r) reveals that there is a small accumulation of particles in the compressed quadrant and a depletion in the extension quadrant, the effect reversing upon flow reversal. Since the particles remain caged, this change of microstructure simply corresponds to a local redistribution of contacts and a mild distortion of the static cages, the overall microstructure being unchanged and reversible after a few oscillations.

At strain amplitudes larger than the yield strain, the particles can escape their position and rearrange over large distances. The cages are continuously advected and renewed by the oscillatory shear flow. The mean square displacements of the particles in the three x-, y-, and z-directions increase linearly in time, which indicates that the shearing oscillatory motion along the x-direction induces a spatially isotropic shear-induced diffusive motion of the particles. The evolution of the microstructure along a cycle is controlled by the value of the shear stress σ relatively to the yield stress σ_y. Along the portions of the cycle where |σ| < σ_y, i.e., where the instantaneous shear rate is small, both the pair correlation function g(r) and the associated harmonic g_2,−2(r) exhibit little changes with respect to the low amplitude case or the static situation. This suggests that the particles are locally trapped in their cages so that the material instantaneously responds like an elastic solid. However, contrary to what happens to the low amplitude case, the particles start to yield and rearrange once the yield stress is exceeded. The evolution of the microstructure is then characterized by an accumulation–depletion mechanism where particles accumulate in the compressive directions of the flow, where they are more compressed, and deplete in the extension directions where they are less compressed. To slide on top of one another and rearrange, particles have to overcome the barrier force at the point of maximum accumulation, where the particles are highly compressed. This mechanism is quite similar to that we previously described and discussed for steady shear flows [Seth et al. (2011)]. This analogy indicates that the portions of the cycles where |σ| > σ_y, i.e., where the instantaneous shear rate is large, correspond to macroscopic shear flow.

In conclusion, we have demonstrated that soft jammed suspensions sheared periodically at large amplitudes undergo a sequence of well-defined physical events: Elastic caging in the low stress portions of the cycle, cage breaking near the yield stress, and flow rearrangements in the portions of the cycle where the stress exceeds the yield stress. This approach was recently proposed from the perspective of rheology by Rogers et al. (2011a). Here, we show that this framework is also relevant with respect to the evolution of the microstructure. To finish, it is interesting to discuss the respective significance of the strain and stress amplitudes. The ratio of the strain amplitude to the yield strain, γ₀/γ_y, discriminates the low and medium/large amplitude regimes. For γ₀/γ_y < 1, particle diffusivity is very low, the microstructure is similar to that at rest, and the mechanical response is in the linear regime. For γ₀/γ_y > 1, particles diffuse through the suspension, the microstructure is distorted with respect to quiescent states, and the mechanical response is the nonlinear regime. It is worth noting that a single quantity, γ₀/γ_y, thus suffices to characterize the onset of yielding at the particle scale (local mobility), mesoscopic scale (microstructure), and macroscopic scale (rheology). Now when γ₀/γ_y > 1, the instantaneous value of the stress relative to the yield stress characterizes the state of the material, i.e., caged (σ/σ_y < 1), yielding (σ/σ_y ∼ 1) or flowing (σ/σ_y > 1), along the cycle.

The elastic component of the shear stress can be computed from the g_2,−2(r) spherical harmonic of the pair correlation function according to

where ρ is the average number density of particles and f^e is the magnitude of the elastic force between two particles given by Eq. (1) [Morris and Katyal (2002); Seth et al. (2011)]. It is interesting to note that the integral in this expression is dominated by the point of maximum accumulation, where they are many particles which are highly compressed [f^e and g_2,−2(r) take their extreme values]. The elastic contribution to the stress can be calculated easily using the expression of the elastic force f^e and the spherical harmonics g_2,−2(r). The results are plotted in Fig. 12, where we observe the data fall onto the BL plots obtained numerically. This result is interesting since it shows that the total stress is dominated by the elastic contribution coming from the microstructure and the elastic repulsive potential, the viscous contribution being negligible. We previously drew a similar conclusion for steady shear flows [Seth et al. (2011)].

We now propose an alternative technique based on the previous result that the material behavior along cycles can be decomposed in a sequence of processes, namely caging, yielding, and flow. We start from the regions of the LAOS cycle where the strain amplitude is near its extrema or the shear rate is zero. From Sec. IV A above, we know that the particles are instantaneously trapped in cages, the cages being only mildly distorted from their static configuration, so that the suspensions essentially respond elastically. Following Rogers et al. (2011a), we define the cage modulus as the local slope of the BL plots at σ = 0

The results are plotted in Fig. 13 for both the simulations and the experiments. Again we note the good agreement between the two sets of data. The other interesting result is that the cage modulus for increasing strain amplitudes remains constant and comparable to the storage modulus at low strain amplitudes. For comparison, the effective storage modulus G calculated over one cycle decreases when γ₀/γ_y > 1. The data in Fig. 13 are obtained for a frequency of $ηω/E*=2 × 10−8$⁠. Similar results are obtained at higher frequencies, the cage modulus being the low strain amplitude storage modulus for that frequency. This means that the lateral segments of the BL plots are the signatures of the cage elasticity and they can be used to determine the magnitude of the linear storage modulus in LAOS experiments.

Let us now turn our attention to the portions of the cycles associated with postyielding behavior (σ > σ_y). From Sec. IV A above, we know that the material flow in a way similar to that in steady shear flows. We can thus use the relation between the shear rate and the stress as an indicator of the flow curve. In Fig. 14, we plot the data obtained for LAOS simulations and experiments at different strain amplitudes and frequencies altogether with the macroscopic flow curves computed and measured independently for steady shear flows [Seth et al. (2011)]. We observe that the data obtained from LAOS rheology are found in good agreement with the steady shear flow curve as shown in Fig. 14. Thus, in the large stress region of the cycle the flow mechanism is similar to that in steady shear where the material exhibits a shear thinning viscoelastic behavior. The unifying effect of the shear-rate amplitude on the particle scale diffusivities exemplified in Fig. 5 must lie at the root of this collapse of LAOS cycles onto to the flow curve. Figure 14 also show that the experimental and simulation data are found in good agreement, which supports the validity of the micromechanical model.

Our results demonstrate the capacity of our micromechanical model to reproduce quantitatively the oscillatory shear rheology of concentrated microgel suspensions. Beyond the good agreement between simulations and experiments for microgel suspensions, our results account for the generic linear and nonlinear rheological behavior of many materials made of soft jammed particles [Ewoldt et al. (2010); Hyun et al. (2011)]. Our work is based on a multiscale study which simultaneously analyzes the dynamics at the particle scale, the statistical properties of the microstructure at the mesoscopic scale, and the rheology at a macroscopic level.

The results concerning the local dynamics of particles echo a former investigation of yielding in concentrated emulsions [Hebraud et al. (1997)]. Using diffusive wave spectroscopy to measure the motion of droplets in concentrated emulsions subjected to oscillatory shear, Hebraud et al. found that the particles reversibly retrace their trajectories at low strain amplitudes but undergo irreversible rearrangements above the yield strain. In our simulations, the plateauing of the mean square displacements at low strain amplitudes indeed indicates that particles do not undergo large scale rearrangements, while the linear increase at large scale amplitudes shows that irreversible rearrangements take place. Hebraud et al. were also able to detect that a finite fraction droplets systematically exhibit chaotic motion while the rest of the sample follows reversible trajectories over many periods. This surprising feature would deserve more attention in our simulations, which will help to understand the physical mechanisms associated with yielding.

Our findings concerning the particle scale dynamics also exhibit interesting analogies with the behavior of non-Brownian hard sphere suspensions in a viscous fluid, which are subjected to oscillatory shear motion [Pine et al. (2005); Corté et al. (2008)]. The physical origin of interparticle interactions is a priori significantly different in both systems: In non-Brownian suspensions, particles interact through long-range hydrodynamic interactions; in dense suspensions of soft particles, particles experience short-range EHD forces at particle contacts. It was found that irreversibility in non-Brownian viscous suspensions occurs above a well-defined critical strain amplitude above which particles diffuse over long distances and below which they organize after a few cycles into a configurations that remain undisturbed. This critical strain in non-Brownian hydrodynamic suspensions thus plays the role of the yield strain in our dense suspensions. This analogy suggests some common underlying physics originating from the existence of many contact interactions, independently of the detailed nature of interparticle forces. In addition, the threshold between the reversible and the irreversible states in hydrodynamic suspensions was recently described as a phase transition. It would be interesting to analyze and push forward this concept in dense suspensions of soft particles, where it may provide a new description of yielding as the onset of irreversibility and unpredictability. The micromechanical model and the simulations presented here can give access to crucial dynamical quantities like the fraction of active particles at any time, the equilibration time, and spatial correlations between particles and help to analyze our conjecture in more detail. Preliminary results show that the fraction of active particles, i.e., those that move irreversibly, is indeed zero at small strain amplitudes and goes to unity beyond the yield strain.

At the mesoscopic level, the statistical properties of the particle distribution provide unambiguous signatures of the onset of yielding. At very low strain amplitudes, the mesoscopic structure adjusts itself to the periodic shearing motion through localized relaxation of contacts, expressing that the deformation of soft particle suspensions is essentially nonaffine [Lacasse et al. (1996)]. At larger strain amplitudes, the occurrence of rearrangements is associated with a periodic modification of the pair distribution function, where particles accumulate upstream where they are more compressed and deplete downstream where they are less compressed. At high shear rates, concentrated hard sphere suspensions also exhibit an accumulation and depletion of particle density along the compression and extension axes, but the stress is linearly proportional to the shear rate and not quadratic as described here [Brady and Morris (1997); Crassous et al. (2008)]. These results indicate that the soft EHD interactions are central to the rheology of deformable particles. The accumulation–depletion mechanism is dominant along the portions of the oscillatory cycles where the stress exceeds the yield stress but negligible elsewhere. In particular the pair distribution function when the stress is zero resembles that at rest, indicating that particles are instantaneously trapped in cages. This shows that the behavior of the material along one oscillation can be analyzed in terms of a sequence of microstructural events that repeat periodically.

At the macroscopic level, the interpretation in terms of a succession of physical processes provides a way of mapping the information obtained from LAOS onto the results of several independent rheological tests like linear viscoelasticity measurements, strain sweep tests, and steady shear flow for both simulations and experiments.

The authors gratefully acknowledge support from the National Science Foundation (Grant No. 0854420) and the computational resources of the Texas Advanced Computing Center (TACC) at The University of Texas at Austin. M.C. thanks Professor Dimitris Vlassopoulos for enlightening discussions about LAOS. L.M. acknowledges support of an Eiffel Excellence scholarship from the French Ministry of Foreign Affairs.