Brownian dynamics simulations of bidisperse hard discs moving in two dimensions in a given steady and homogeneous shear flow are presented close to and above the glass transition density. The stationary structure functions and stresses of shear-melted glass are compared quantitatively to parameter-free numerical calculations of monodisperse hard discs using mode coupling theory within the integration through transients framework. Theory qualitatively explains the properties of the yielding glass but quantitatively overestimates the shear-driven stresses and structural anisotropies.
Shear flow can drive dense colloidal dispersions into states far from equilibrium. Especially of interest is the possibility to shear melt colloidal solids, in particular metastable colloidal glasses and gels, and to investigate shear-melted (yielding) colloidal glasses. Does a yield stress and/or yield strain exist (Petekidis et al. 2002)? Are shear-melted states necessarily heterogeneous (e.g. shear banded)? Does ageing prevent stationary states under steady shearing? Do hydrodynamic interactions dominate the properties in flow? These are just some of the questions whose answers will provide insights into the still murky glassy state.
Recently, the integration through transients (ITT) approach has been used to generalize mode coupling theory (MCT) of structural glass transition to the case of steady shearing (Fuchs & Cates 2002, 2009), and to arbitrary homogeneous flows (Brader et al. 2008), albeit based on a number of approximations. First, homogeneous flow profiles are assumed from the outset; second, hydrodynamic interactions are neglected; and third, the mode coupling (mean-field-like)–decoupling approximation, splitting a four-point density fluctuation function into the product of two-point density correlators, has been applied. While simplified (schematic) models of MCT-ITT (Fuchs & Cates 2003) can well be fitted to linear rheology (frequency-dependent storage and loss moduli) and nonlinear rheology (flow curves, viz. stress as a function of shear rate) of model dispersions (Siebenbürger et al. 2009), no fully quantitative comparison between the theory and data without additional approximations has been presented up to now.
We provide this first quantitative comparison, solving the MCT-ITT equations for hard discs confined to move in a plane, and performing Brownian dynamics simulations for a binary hard disc mixture, also in two dimensions. Computational (memory) constraints force us to work in two dimensions, which is still of some experimental relevance, as glasses have been observed in quasi-two-dimensional dispersions (Bayer et al. 2007). Use of a binary mixture suppresses the nucleation rate sufficiently for even our long simulation runs to remain free of crystal nuclei. As hydrodynamic interactions are absent also in the simulation, and as Lees–Edwards (LE) boundary conditions in combination with the thermostat impose a homogeneous shear rate, the integration of the equations of motion in the simulation ensure a stringent test of the theory. Results for stationary shear and normal stresses, and for the shear-distorted microstructure are presented and discussed. The input quantities entering the theory, and transient density correlation functions, which (crucially) enter during intermediate MCT-ITT steps, are also characterized.
Our work bears some similarity to the study by Miyazaki and Reichmann, who, however, concentrated on fluid states under shear, and on their time-dependent fluctuations (Miyazaki & Reichmann 2002). We focus here on the yielding glass state, and its stationary, time-independent averages, which are not accessible to the theory in Miyazaki & Reichmann (2002). In order to compare quantitatively with theory, we aimed for better statistics than in comparable previous two-dimensional simulation studies of sheared glasses (Yamamoto & Onuki 1998; Furukawa et al. 2009).
The basic concept of the algorithm has been described in detail in three dimensions by Scala et al. (2007) and can easily be adapted to two dimensions x and y. In order to prevent crystallization at high densities, we consider a binary mixture with a diameter ratio of D2/D1=1.4 at equal number concentrations (Nsmall=Nbig). N=1000 hard discs move in a simulation box of volume V with periodic LE boundary conditions at packing fraction . The mass of the particles is set as equal to unity, thermal energy kBT sets the energy scale, and D1≡D is used as the unit of length in the following. We choose our coordinate axes such that the flow is in the x-direction and the shear gradient is in the y-direction. After putting the particles on their initial positions, we provide Gaussian distributed velocities that will be overlain by the linear shear flow. To propagate the system forwards in time, we employ a semi-event-driven algorithm. For every particle at time t0, the algorithm determines the possible collision time τ with any other particle. This is easily achieved by solving the equation 2.1 where Di and Dj denote the diameters of the particles i and j, rij denote the relative vectors between both particles and vij denote the relative velocities. The smallest solution for all particle pairs determines the next event in the algorithm. All particles can then be propagated according to . With conservation of energy and momentum, binary collision laws impose new velocities for particles j and i, 2.2 Owing to the boundary conditions, any particle in the vicinity of the box boundary can collide with an image particle coming from the other end of the box. The boxes are constructed in such a way that they are translated with velocity , where L denotes the size of the box. The boundary conditions thus have the consequence that the velocities tend to increase. Therefore, a thermostat has to be introduced. After a time τB=0.01, a so-called Brownian step sets in, which assures that the particles move diffusively for longer times. In the Brownian step at time τB, all particle velocities are freshly drawn from a Gaussian distribution with 〈|v|2〉=2. After that, linear shear flow is imposed so that holds.
As the system starts from a cubic lattice, it is necessary to wait for the system to relax before meaningful stationary averages can be taken. The quantity of interest in the work presented here is the potential part of the stress, , with the relative force components (in the α-direction) of particles i and j (F(t)ij)α, and the particles’ relative distance components (r(t)ij)β. The kinetic part plays no role, and thus has already been omitted. As we consider hard particles, the forces must be calculated from the collision events. By observing the collisions within a certain time window Δτc, forces may be extracted using the change of momentum that occurs during the observation time. This leads to the evaluation algorithm (Foss & Brady 2000; Lange et al. 2009), 2.3 where the summation is over all collisions of particles i and j at time tc within the time window Δτc. The procedure effectively sums the momentum changes (Δvij)α in the α-direction multiplied by the relative distance of the particles (rij)β in the β-direction. Here and below, the brackets 〈… 〉s denote the average over different simulation runs.
Additionally, the shear stress can be computed via the contact value g(r), 2.4 where g11(r,θ),g12(r,θ),… denote the θ-dependent partial contact values of the two components and d11,d12,… denote the minimal distance between two particles. θ is the polar angle.
The Green–Kubo relation holds for the non-sheared system. Thus, the shear viscosity can be extracted from the simulation via (Alder et al. 1970) 2.5 where the sum runs over all collisions. The second pivot in this analysis is the equal-time structure factor and its deviation from the quiescent system. Exploiting the fact that is a real quantity, we can extract it via 2.6 where the double sum runs over all pairs of particles i and j (NI). The pairs ij are determined by the LE boundary conditions and the constraint of having the lowest distance among all possible image particles in the surrounding boxes.
For low Peclet numbers () in fluid states, g(r) can be expanded: . This result can be used to derive the relative distortion of the structure factor in the linear response regime (Strating 1999), 2.7 while g1(r) can be obtained via 2.8
3. Mode coupling theory in the integration through transients framework
The MCT-ITT approach generalizes the MCT of the glass transition to colloidal dispersions under strong continuous shear. It considers N spherical particles with arbitrary interaction potential that move by Brownian motion relative to a given linear shear profile. An equation of motion for a transient density correlator Φq(t) encodes structural rearrangements, and approximated generalized Green–Kubo laws relate stress relaxation to the decay of density fluctuations.
The transient density correlator is defined by Φq(t)=〈ρ*qρq(t)(t)〉/NSq, where the density fluctuation is as usual , and, in ITT, the average can be performed over the equilibrium Gibbs–Boltzmann ensemble. Thus, the normalization of the initial value Φq(0)=1 is given by the equilibrium structure factor Sq. The time-dependent or shear-advected wavevector appearing in the definition eliminates the affine particle motion with the flow field, and gives Φq(t)≡1 in the absence of Brownian motion. Shear flow coupled to random motion causes Φq(t) to decay, as given by an (exact) equation of motion containing a retarded friction kernel that arises from the competition of particle caging and the shear advection of fluctuations, 3.1 where the initial decay rate contains Taylor dispersion denoted by Γq(t)=q2(t)/Sq(t). The generalized friction kernel mq(t,t′), which is an autocorrelation function of fluctuating stresses, is approximated, following MCT precepts by an expression involving the structural rearrangements captured in the density correlators 3.2 with abbreviation p=q−k, n is the particle density and a vertex function is given by 3.3 where ck is the Ornstein–Zernicke direct correlation function ck=(1−1/Sk)/n. An additional memory kernel is neglected (Fuchs & Cates 2009). The equilibrium structure factor, Sk, encodes the particle interactions and introduces the experimental control parameters such as density and temperature. Similarly, the potential part of the stress in the non-equilibrium stationary state (neglecting the diagonal contribution that gives the pressure) is approximated assuming that stress relaxations can be computed from integrating the transient density correlations 3.4 Flow also leads to the build-up of shear-induced microstructural changes, which, again integrating up the transient density correlators, can be found from 3.5 A far smaller isotropic term in (see Fuchs & Cates 2009) is neglected here, as it is of importance for the plane perpendicular to the flow only.
4. Numerical aspects
The set of coupled equations (3.1)–(3.3) was solved self-consistently using modifications of standard algorithms and the Ng iteration scheme (Ng 1974). The functions were discretized on a two-dimensional Fourier grid consisting of 101 grid points in either coordinate direction using a cutoff of qD=30. More explicitly, the Fourier grid was discretized as qx,yD=−30,−29.4,…,0,…,29.4,30, with D being the particle diameter. In order to enhance the accuracy, the advected direct correlation function cq(t) in equation (3.3) was calculated from the input structure factor once it has been advected beyond the Fourier-grid cutoff up to qD=100. For the time integration of the convolution integral in equation (3.1), an initial time step of 10−9D2/D0 was used. The time integration in the calculation of derived quantities (3.4) and (3.5) was performed by dividing the integration interval into two sub-intervals, the first one containing times shorter and the second one containing times longer than . The integration was then carried out on the intrinsic quasi-logarithmic grid of the MCT equation solver for times 0≤t≤t* and on an equally spaced linear time grid with spacing Δt=t* for times t>t*.
5. Equilibrium and transient quantities
The MCT-ITT approach uses structural equilibrium correlations as input, computes transient structural density correlators to encode the competition between flow-induced and Brownian motion, and calculates all stationary properties from time integrals over the transient fluctuations. In this section, the input and intermediate quantities of the theoretical calculations are presented and discussed.
(a) Equilibrium structure factor
The equilibrium structure factor Sq is the only input quantity to the MCT-ITT equations. It varies smoothly with density or temperature, but leads to the transition from a shear-thinning fluid to a yielding glassy state at a glass transition density φc. Figure 1 shows modified hyper-netted chain structure factors of monodisperse hard discs in d=2 from Bayer et al. (2007) used in the MCT-ITT calculations. For comparison, the averaged structure factors obtained from the binary mixture simulation are also shown. In both cases, the value of the glass transition density φc is included in the curves in figure 1, and only smooth changes in Sq are noticeable. The short-range order at the average particle distance, as measured in the primary peak of Sq, increases with densification. Differences in the structure between the monodisperse and the bidisperse system become appreciable beyond the primary peak, and especially beyond the second peak in Sq. Also, the height of the primary peak in Sq at φc differs, indicating that the averaged structure factor in the bidisperse system does not well characterize the local structure and caging in the simpler system, and that a comparison with a bidisperse MCT-ITT calculation should be performed. As this is outside the present numerical reach, we proceed by comparing results for the characteristic wavevectors denoted by q1 to q4 in figure 1. The angular dependence of the anisotropic structure will be explored along the special directions indicated in the insets.
(b) Transient coherent density correlation functions
The central quantities in MCT-ITT, which encode the competition between shear-driven motion and random fluctuations, are the transient density correlators. Structural rearrangements manifest themselves as a (second) slow relaxation process in the Φq(t), whose relaxation time depends sensitively on the distance to the glass transition, measured by the (relative) separation parameter ε=(φ−φc)/φc, and on the magnitude of the shear rate. Figure 2 shows representative curves in a fluid state in figure 2a(i)–(iv) and in a shear-melted glassy state in figure 2b(i)–(iv) for the wavevectors and directions defined in figure 1, and for nine different shear rates. The denoted bare Peclet numbers measure the shear rate compared to the Brownian diffusion time estimated with the single-particle diffusion coefficient D0 at infinite dilution. Characteristically, for the present strongly viscoelastic system, shear affects the structural relaxation already at extremely small bare Peclet numbers Pe0≪1. The short time motion, which corresponds to the local diffusion of the particles within their neighbour cages, however, is not much affected as long as Pe0<1 holds.
In the fluid (ε<0), a linear response regime, where shear does not affect the decay of thermal fluctuations, is observed for small (dressed) Peclet or Weissenberg numbers Pe≤1, where measures the shear rate relative to the intrinsic relaxation time. The latter can be estimated from the relaxation of Φq2(t), viz. at the primary peak of Sq, where the structure relaxes most slowly in the quiescent case. The final (or α-) relaxation time τ increases strongly when approaching the glass transition with predicted by MCT (Bayer et al. 2007). For the fluid state in figure 2, only Pe0≤10−6 is small enough so that Pe<1 holds, and the final relaxation curves agree for different shear rates. In the glass (ε≥0), shear is the origin of the decay of the otherwise frozen-in density fluctuations, and all nine shear rates lead to different final decays. While for small and larger wavevectors the final decay is rather isotropic, around the main peak in Sq, some anisotropy is noticeable in the transient correlators. In the direction perpendicular to the flow, qx=0 (red), shear is not as efficient in decorrelating the density as in the other directions, which fall closer together.
One of the central predictions of MCT-ITT concerns the existence of a scaling law describing the yielding of glassy states, which manifests itself by an approach to a master function for decreasing shear rate, and ε≥0, where the rescaled time agrees with the accumulated strain. This scaling law, which is quite obvious in figure 2, is tested quantitatively for some wavevectors in figure 3, and holds extremely well (for q1), or with pre-asymptotic corrections (for q2). The shapes of the yielding master functions can be very well fitted with compressed exponentials, , but the parameters, including the exponent βq, depend on wavevector and orientation; see table 1 for representative values.
The quite good fit to compressed exponentials is at present only a numerical observation. Yet, close to the glass transition, asymptotic expansions are possible that analytically predict the initial part of the yield master functions (Fuchs & Cates 2002, 2003, 2009), As the critical glass form factor and the critical amplitude hq are isotropic, this result suggests a rather isotropic yielding process right at the glass transition density. Approximating by an exponential also provides an estimate for the final relaxation time under shear. Except for , all quantities in the above formula have been determined in quiescent MCT (Bayer et al. 2007), and our values (e.g. λ=0.698) differ only because of the coarser discretization in q space that is necessary under shear. Close to the glass transition, we estimate . As figures 2a and 3a show, this estimate for the final relaxation time is quite good close to the transition, but deeper in the glass, the correlators become somewhat more anisotropic. Specifically for qx=0 (red) with wavevector magnitudes around |q|≃q2, the correlator slows down relative to other orientations.
More detailed analytical predictions are possible around , where spatial and temporal dependences in the transient fluctuations decouple (Fuchs & Cates 2003), is a universal function, which only depends on λ, C=2.08 and is another MCT time scale that diverges at the glass transition (Bayer et al. 2007). This factorization property known from quiescent MCT generalizes to steady shear, and is an essential feature of the localization transition that underlies glass formation in MCT. In figure 4, the quantity Xq(t)=(Φq(t)−Φq(t′))/(Φq(t′)−Φq(t′′)) is plotted, which should become wavevector and orientation independent if factorization holds. This holds in the liquid and in the glass where, however, shear leads to strong (anisotropic) pre-asymptotic corrections already for .
6. Results and comparisons
Based on the approximated generalized Green–Kubo relations of equations (3.4) and (3.5), and the properties of the transient correlators discussed in §5, MCT-ITT makes a number of predictions for stationary stresses and structural correlations (Fuchs & Cates 2002, 2009). In the following, we will test these not only qualitatively, but also quantitatively, by comparing them to the two-dimensional simulation data.
(a) Stationary stresses
The quantity of most interest in nonlinear rheology is the shear stress . Flow curves giving the shear stress as a function of the shear rate can be obtained in the simulations and are shown in figure 5. The viscosity decreases below its Newtonian value upon increasing the shear rate already for small Pe0 values (figure 6). Shear thinning in which the stress increases less than linearly with , sets in at a Pe of the order of unity, and this crossover shifts to lower and lower Pe0 for increasing density. At the density around φc≈0.79, the crossover leaves the accessible shear-rate window. We use this as an estimate of an ideal glass transition density, where the final relaxation time and the quiescent Newtonian viscosity () diverge. Moreover, at this density, the flow curves change from a characteristic S shape in the fluid to exhibiting a rather -independent plateau at small shear rates. This change is the hallmark of the transition in MCT-ITT between a shear-thinning fluid and a yielding glass (Siebenbürger et al. 2009).
The numerical MCT-ITT solutions show the same transition scenario, but because of the approximations involved, exhibit a different critical density than that of the simulations. Even if theory and simulation considered the same (binary) system, a difference in critical packing fraction would be expected, as is well known also in three dimensions (Voigtmann et al. 2004). Thus, the following MCT-ITT calculations match the relative separation ε from φc in order to compare with the simulations. As MCT-ITT is aimed at describing the long-time structural motion, errors need to be anticipated in its description of short-time properties. This is obvious in real dispersions, where hydrodynamic interactions (neglected in MCT-ITT) affect the short-time diffusion coefficient Ds, but may also arise in the following comparisons. A rescaling of the effective Peclet number would correct for this change in Ds. In order not to introduce additional fit parameters we refrain from doing so, but anticipate that future comparisons may require Ds≠D0.
The approach to a yield scaling law, where the final decay of the transient correlators depends on the accumulated strain only (see figures 2 and 3), predicts the existence of a (dynamical) yield stress , which characterizes the shear-melted glass. In the bidisperse hard disc mixture, at the glass transition, it takes the (critical) value ; below the glass transition, the yield stress jumps to zero, . Its quantitative prediction is quite a challenge for theory, because equation (3.4) shows that it requires an accurate calculation of the shear-driven relaxation process. MCT-ITT overestimates the critical yield stress by roughly a factor of 10 because, presumably, the decay of the transient correlators is too slow. Yet, of course, the difference between the monodisperse system in the MCT-ITT calculation and the bidisperse simulated system contributes in an unknown way to the error. It appears reasonable to assume that mixing two species reduces the stresses under flow, which would explain part of the deviation.
At larger shear rates, the flow curves from simulation appear to approach a second Newtonian plateau that, presumably, strongly depends on the hard-core character of the excluded volume interactions and is outside the reach of the present MCT-ITT. The latter, by using its sole input Sq rather than pair potential, is not directly aware of hard-core constraints. We checked, however, that the states remain homogeneous and random up to the Pe0 values shown.
Reassuringly, the same rescaling factor of 0.1 as for the shear stress brings theoretical and simulational normal stress differences, σxx−σyy, into register also; see figure 7. The normal stress differences are positive (indicating that the dispersion would swell after flowing through a nozzle), and show similar behaviour to the stress, increasing as in the fluid for small shear rates, and levelling out onto a plateau in the yielding glass.
(b) Distorted microstructure
The macroscopic stresses in the flowing dispersion are experimentally most important, but provide only an averaged description of the local effects of shear. Spatially resolved information can be learnt from the distorted structure factor , which in MCT-ITT is connected to the stress via as follows from equations (3.4) and (3.5).
Figure 8 shows colour-coded structure factors as functions of the two-dimensional wavevector q, with qx being in the direction of flow and qy being along the gradient direction. In (a), panels with simulation data are compared to panels in (b) obtained in MCT-ITT. Scattering intensities are presented: (i) (figure 8a(i),b(i)) in the linear response regime in the fluid, effectively measuring the equilibrium structure factors already shown in figure 1, (ii) (figure 8a(ii),b(ii)) at high shear in the glass, where all densities are in the shear thinning region and (iii) (figure 8a(iii),b(iii)) in the glass at low shear rate, where the yielding glassy state is tested. While the fluid is isotropic for small Pe0 (case (i)), as required by linear response theory, increasing Pe0 to values around unity leads to an ellipsoidal scattering ring, which is elongated along the so-called ‘compressional axis’ qx=−qy and more narrow along the ‘extensional axis’ qx=qy. This indicates that shear pushes particles together along the compressional diagonal and pulls particles apart along the extensional diagonal (Vermant & Solomon 2005). Theory and simulation data in this representation qualitatively agree, except for the fact that MCT-ITT somewhat overestimates the anisotropic distortion of the glass structure at low Pe0.
A more careful look onto the distorted microstructure is provided by q-dependent cuts through along the directions colour coded in figure 1. Important for the direction qx=0 (red) is the need, present in both simulations and theory, to average over a small but finite angle, because exactly at qx=0, the structure oscillates noisily around zero.
Especially of interest are the intensities of case (iii), where the stationary structure of the shear-melted solid is studied. The density φ=0.79 is not yet high enough to lie in the glass, but close enough to the glass transition so that the correlations at a quite low rate, namely Pe0=2×10−4, closely resemble the ones of glassy states at very low shear rates. (While we obtain good statistics for stresses at Pe0=2×10−5 for all densities, structure factors cannot be sampled sufficiently there.) At this density of φ=0.79, the equilibrium structure factor Sq can also be obtained in long simulation runs, and the (relative) difference can thus be determined. This is shown in figure 9 for two states in order to investigate the distorted structure of the shear-melted glass in detail. We base the following discussion on the hypothesis that Pe0=2×10−4 at φ=0.79 captures a glass-like state in the limit of low bare Peclet number.
Figure 9 provides a sensitive test of the accuracy of the theoretical predictions. The lower left panel shows that the structure factor at vanishing shear rate jumps discontinuously at the glass transition; while holds in the fluid, holds in the glass. Relative deviations of 20 per cent remain. Simulation finds quite isotropic deviations that show a maximum on the low-q side of the primary peak in Sq. MCT-ITT predicts the absence of a linear response regime in as a function of the shear rate in the glass, and derives it from the existence of the yield scaling law in the transient correlators. Because sets the time scale for the final relaxation into the stationary state, the limit does not agree with the quiescent result .
Quantitatively, MCT-ITT overestimates the distortion again by a factor of up to 10, and finds a noticeable anisotropy, as discussed in context with figure 8. While the difference between the bidisperse and the monodisperse systems may influence the comparison, we believe that the major origin of the error is that MCT-ITT underestimates the speeding up of structural rearrangements caused by shear. The too slow transient correlators thus become anisotropic because the accumulated strain becomes too big before structural correlations have decayed.
Qualitatively, aspects of the anisotropy predicted by MCT-ITT can be seen in the simulations at only slightly larger shear rates, such as at Pe0=0.2 shown in figure 9a(ii). While along the two axis and the extensional diagonal direction, the low-q wing of the primary peak in becomes enhanced under shear, along the compressional axis (magenta), it gets suppressed, and the high-q wing is pushed up. Simulation also finds a suppression of the peak height along all directions, which MCT-ITT reproduces along the diagonal directions. Overall, the anisotropy and the magnitude of the distortions predicted by MCT-ITT remain too large, but the deviation decreases. The differences between the equilibrium structure factors Sq in the simulated and calculated systems should be taken into account in future work, and presently preclude comparisons at larger wavevectors.
7. Conclusions and outlook
The first fully quantitative solutions of the MCT-ITT equations for the distorted microstructure and the stresses in steadily sheared two-dimensional shear-thinning fluids and yielding glasses of Brownian hard discs exhibit all the universal features discussed within schematic MCT-ITT models (Fuchs & Cates 2003). They compare well qualitatively, but with appreciable errors quantitatively, with Brownian dynamics simulations of a bidisperse mixture without hydrodynamic interactions in a linear shear profile, which, for all states considered, remains in a homogeneous and disordered state. The non-analytic behaviour of the stationary properties, and the lack of a linear response regime throughout the (shear-melted) glass state, predicted by theory, can be found in the simulation.
This work was funded in part by DFG-IRTG 667 and EPSRC grants EP/045316 and EP/030173. M.E.C. holds a Royal Society Research Professorship.
One contribution of 12 to a Discussion Meeting Issue ‘Colloids, grains and dense suspensions: under flow and under arrest’.
- © 2009 The Royal Society