Structural relaxation and rheological response of a driven amorphous system
The interplay between the structural relaxation and the rheological response of a simple amorphous system (an 80:20 binary Lennard-Jones mixture [W. Kob and H.C. Andersen, PRL 73, 1376 (1994)] is studied via molecular dynamics simulations. In the quiescent state, the model is well known for its sluggish dynamics and a two step relaxation of correlation functions at low temperatures. An ideal glass transition temperature of has been identified in the previous studies via the analysis of the system’s dynamics in the frame work of the mode coupling theory of the glass transition [W. Kob and H.C. Andersen, PRE 51, 4626 (1995)]. In the present work, we focus on the question whether a signature of this ideal glass transition can also be found in the case where the system’s dynamics is driven by a shear motion. Indeed, the following distinction in the structural relaxation is found: In the supercooled state, the structural relaxation is dominated by the shear at relatively high shear rates, , whereas at sufficiently low the (shear-independent) equilibrium relaxation is recovered. In contrast to this, the structural relaxation of a glass is always driven by shear. This distinct behavior of the correlation functions is also reflected in the rheological response. In the supercooled state, the shear viscosity, , decreases with increasing shear rate (shear thinning) at high shear rates, but then converges toward a constant as the is decreased below a (temperature-dependent) threshold value. Below , on the other hand, the shear viscosity grows as suggesting a divergence at . Thus, within the accessible observation time window, a transition toward a non-ergodic state seems to occur in the driven glass as the driving force approaches zero. As to the flow curves (stress versus shear rate), a plateau forms at low shear rates in the glassy phase. A consequence of this stress plateau for Poiseuille-type flows is demonstrated.
Although simple at first sight, suspensions of spherical colloidal particles under shear exhibit a rich phenomenology. In the dilute regime, at temperatures corresponding to the liquid state, forced Rayleigh scattering experiments Qiu-PRL61-1988 () show an increase of diffusion constant upon shearing (shear thinning), distinct from Taylor dispersion Taylor-ProcRSocLondonA-219::1953 (). Indeed, Taylor dispersion (displacement of particles along the flow direction, , as they move in the direction of shear gradient) would give rise to an enhancement of mobility via , whereas shear thinning leads to , with a diffusion coefficient, , which increases with increasing in a non-trivial fashion.
When shearing starts in the crystalline phase, however, shear induced melting of the crystalline structure is observed both in light scattering experiments Ackerson-Clark::PRL46-1981 (); Ackerson::JRheol34-1990 () as well as in Brownian dynamics simulations Xue-Grest::PRL64::1990 (). Furthermore, in the shear melted regime, experiments show evidence for shear thinning due to the presence of freely slipping two dimensional crystalline layers Ackerson-Clark::PRL46-1981 (); Ackerson::JRheol34-1990 (). Simulations also show a shear thinning regime below a ”critical” shear rate, , followed by a transition to a string-like order for Xue-Grest::PRL64::1990 ().
On the other hand, studies of disordered suspensions of hard spheres show that shear thinning and shear melting phenomena may also occur in the absence of a crystalline structure Laun-et-al::JRheol36-1992 (); Petekidis::Physica-A306-2002 (); Petekidis::PRE66-2002 (); Petekidis-Vlassopoulos-Pusey::FaradayDiscuss::2003 (). Similar observations have also been made in light scattering echo studies of (disordered) dense emulsions Hebraud-Lequeux-Munch::PRL78-1997 (). Brownian dynamics simulations show that shear thinning in concentrated colloidal suspensions is related to the fact that, in the limit of low shear rates, the main contribution to the shear stress originates from the Brownian motion of colloidal particles and that this contribution decreases with Phung-Brady-Bossis-JFluidMech::1996 ().
Many theoretical approaches have been proposed in recent years aiming at a description of the rheology of disordered materials. Some examples are studies of driven -spin glasses Berthier-Barrat-Kurchan::PRE (), the so called “soft glassy rheology” model (SGR) Sollich::PRE-1998 () and non-equilibrium extensions of the mode coupling theory (MCT) of the glass transition Fuchs-Cates::PRL-2002 (); Fuchs-Cates::FaradayDiscuss-2003 (); Miyazaki-Reichman-PRE66 (); Miyazaki-Reichman-Yamamoto-PRE70 ().
These approaches have the common feature of predicting (or, more precisely, “reproducing”, as the phenomenon was known before the theories have been developed) the shear thinning phenomenon in an ergodic disordered system, i.e. in a system where temporal correlations decay to zero in a finite time. In these systems, shear thinning is commonly attributed to the competition between a typical structural relaxation time (of the quiescent system) and the time scale imposed by the external drive. A shear thinning behavior is expected as soon as the latter becomes significantly shorter than the former, in which case the structural relaxation is said to be driven by the external force Berthier-Barrat-Kurchan::PRE (). This idea is also supported by the observation that even a simple liquid may exhibit shear thinning provided that the shear rate is high enough in order for shear disturbance to dominate the system dynamics Evans-Morriss::StatMechNonEqLiq1990 ().
However, the statements of different theories diversify as it comes to the rheological response of a “jammed system” such as an amorphous solid. The discrepancy is best seen when considering the behavior of the yield stress, i.e. the stress response in the limit of vanishing shear rate (driving force), . While calculations based on the -spin glasses in the thermodynamic limit predict that no yield stress exists in the system Berthier-Barrat-Kurchan::PRE () ( vanishes with a power law, with ), the SGR model Sollich::PRE-1998 () predicts the onset of a (dynamic) yield stress at the “jamming transition”, . Here, is a noise temperature which controls the ”degree of jamming” or the distance from the glass transition. corresponds to the glass transition (or “jamming”) temperature, and characterizes the glassy or “jammed” phase (for noise temperatures above the transition, , a power law decrease of the shear viscosity with the applied shear rate is found Sollich::PRE-1998 ()).
Interestingly, the presence of a finite yield stress in an amorphous solid has also been predicted within the simplest (or idealized) version of the non-equilibrium MCT approach Fuchs-Cates::FaradayDiscuss-2003 (); Fuchs-Cates::PRL-2002 (). A related MCT approach to the fluctuations around the steady state has recently been proposed in reference Miyazaki-Reichman-PRE66 (); Miyazaki-Reichman-Yamamoto-PRE70 (). The issue of yield stress, however, could not be addressed in that approach.
Despite the above mentioned qualitatively different predictions on the stress response of a glass, both the -spin approach and the non-equilibrium MCT commonly predict that the presence of an external drive leads to a melting of the amorphous solid.
In reference Berthier-Barrat-Kurchan::PRE () it is argued that the main reason for the absence of a yield stress in the -spin model lies in its mean-field nature so that free energy barriers are impenetrable. If the system is quenched from a high temperature to a temperature in the glassy phase while at the same time being driven by a nonconservative force, it remains sliding above the free energy threshold below which the free energy surface is split into exponentially many disconnected regions Berthier-Barrat-Kurchan::PRE (). This threshold does not exist at high temperatures but emerges at a critical temperature, , separating the liquid like region from the glassy phase. On the other hand, if the system is prepared in an energy state above the free energy threshold, it starts evolving towards the threshold free energy as the temperature is reduced below . This process slows down continuously due to the decreasing connectivity of the visited free energy landscape so that the free energy threshold is never reached. Therefore, even a vanishingly small driving force is sufficient in order to stop aging and keep the system sliding over this threshold free energy Berthier-Barrat-Kurchan::PRE ().
This picture motivated a later work, where the existence of a yield stress could be shown for a -spin system with a finite number of spins Berthier2003 (). The work follows the idea that free energy barriers are finite at finite system size, thus allowing the thermal activations to play a role which is not possible in the case of an infinite system. In a system with finite free energy barriers, thermal activation becomes important since it allows the system to leave local free energy minima where it is trapped so that deeper valleys can be reached. A finite force must then be applied to the system in order to stop this aging process and to impose a steady sliding motion. For , the finite-sized version of the model has been investigated by Monte Carlo simulations supporting the existence of a critical driving force below which the system is trapped (’solid’) and above which it “flows” (’liquid’) Berthier2003 ().
A manifest of shear melting of an amorphous solid can be seen in the behavior of particle displacements. While in the quiescent state a particle is practically eternally trapped in its local environment (the ’cage” formed by its nearest neighbors) it can explore far larger distances as the system is exposed to an external drive such as a shear motion. Thus, the particles in a driven glass are not localized but behave as in a fluid where the whole available space can be explored via diffusive motion.
This “freedom” in particle motion is also reflected in the behavior of the structural correlations. While in the quiescent system the structural relaxation times continuously grow as the system is quenched into a glassy phase (aging), eventually exceeding available experimental time window, all correlations decay to zero within a finite time in a driven glass.
Based on the above picture of shear-driven melting of an amorphous solid, one could possibly expect that, regardless of the shear rate, a finite stress must always be applied in order to enable the system to overcome free energy barriers so that structural relaxation can take place. If this conclusion is true, a stress plateau is expected as soon as the time scale of the inherent structural relaxation is far beyond the accessible simulation time, since, in this case, the finite relaxation time in the driven system can be considered as uniquely imposed by the external drive, a situation very similar to shear-induced melting.
In order to test this idea, we performed large scale molecular dynamics simulations of a simple molecular glass (see below for a description of the model) focusing on the interplay between the structural relaxation and the rheological response of the system both in the supercooled state and in the glassy phase. As published in a recent paper Varnik2006 (), our simulation results clearly support the existence of a stress plateau at low shear rates in an amorphous solid. The focus of that paper was on the stress response and on a comparison of simulated results with MCT-based theoretical calculations. Here, we go a step further and study the dynamics of structural relaxation in the same range of temperatures and shear rates for which the flow-curves have been studied in reference Varnik2006 (). Furthermore, we will investigate an interesting consequence of the presence of a stress plateau (which in this case plays the role of a yield stress) on a Poiseuille-type flow in a planar channel. It will be shown that in a region around the channel center where the stress is below the plateau stress the system behaves as an amorphous solid whereas it exhibits liquid-like behavior beyond this central region.
It is noted that the rheological response of exactly the same model has already been studied in reference Berthier-Barrat::JChemPhys2002 (). However, due to a rather restricted range of studied shear rates, no conclusion could be made on the existence or not of a dynamic yield stress in the glassy phase. The novelty of our simulations is the extended range of shear rates allowing a conclusive statement on the existence of a stress plateau for the present model in the glassy state. When restricted to the range of shear rates accessible to our simulations, this stress plateau may be regarded as the yield stress of the system.
The paper is organized as follows. After an introduction of the model and the simulation method in the next section, the effect of shear on the structural relaxation of an amorphous system will be discussed in section III. In particular, it will be shown that the structural relaxation of an amorphous system exposed to a shear exhibits signature of the quiescent state so that a study of the structural relaxation alone is sufficient in order to find out whether the quiescent state of the (driven) system belongs to a liquid-like (supercooled) state or to a glassy phase. In section IV the rheological response of the system is studied. The focus of this section is to show that a similar information on the quiescent state of the system may also be gained via a study of the stress response. This is possible since the stress-shear rate curves of a glass are qualitatively different compared to those of a supercooled liquid. The information on the glass transition is, therefore, also encoded in the rheological response. Section V compiles our results.
Ii model and simulation method
A generic glass forming system, consisting of an 80:20 binary mixture of Lennard-Jones particles (whose types we call A and B) at a total density of and in a cubic box of length ( particles) is studied Kob-Andersen::PRL73::1994 ().
A and B particles interact via , with , , , , and . The potential was truncated at twice the minimum position of the LJ potential, . The parameters , and define the units of energy, length and mass. All other quantities reported in this paper are expressed as a combination of these units. The unit of time, for example, is given by and that of stress by . Equations of motion are integrated using a discrete time step of .
The system density is kept constant at the value of for all simulations whose results are reported here. This density is high enough so that no voids occur at low temperatures and low enough so that system dynamics remains sensitive to a variation of temperature (see references Starr-PRE60-1999 (); Sastry-PRL85-2000 () for effects of high density/pressure on the liquid-glass transition).
The present model was found suitable for an analysis of many aspects of the mode coupling theory of the glass transition Kob-Andersen::PRE51 (); Kob-Andersen::PRE52 (). In particular, at a total density of , equilibrium studies of the model showed that the growth of the structural relaxation times at low temperatures could be approximately described by a power law as predicted by the ideal MCT, . Here, is the mode coupling critical temperature of the model and is the critical exponent. For the present binary Lennard-Jones system, numerical solution of ideal MCT equations yields a value of Nauroth-Kob::PRE55 (). A similar value is also obtained for a binary mixture of soft spheres Barrat-Latz:JPCM2-1990 ().
Simulation results are averaged over 10 independent runs. For this purpose, ten independent samples are equilibrated at a temperature of (above ) and serve as starting configurations for all simulated temperatures and shear rates. The temperature is controlled via Nosé-Hoover thermostat Nose::JCP81 (); Hoover::PhysRevA31 (). It is set from to the desired value at the beginning of shear, whereby only the -component of particle velocities is coupled to the heat bath ( being the streaming and the shear gradient directions, see also below).
The temperature quench is done only in one step, i.e. without a continuous variation from to . However, as the numerical value of is changed, it takes a time of the order of the velocity autocorrelation time for the new temperature to be established. During this period of time the Maxwell distribution of velocities undergoes changes in order to adapt itself to a distribution determined by the new temperature. This time is of order unity (in reduced units) and quite short compared to all other relevant timescales in the problem.
Previous studies of the stress-strain relation of the same model showed that the initial transient behavior is limited to strains below 50% Varnik-Bocquet-Barrat::JChemPhys2004 (). Indeed, by shifting the time origin in measurements of various correlation functions, we verified that the time translation invariance was well satisfied in sheared systems for strains larger than 50%. We neglected strains before starting the measurements. Unless otherwise stated, all simulations reported below had a length corresponding to a strain of 7.8 (780%). In the steady state, correlation functions were averaged both over independent runs and over time origins distributed equidistantly along each simulation run. The shear stress is calculated using the virial expression Evans-Morriss::StatMechNonEqLiq1990 ()
where stands for statistical averaging, and are the velocity and the mass of -th particle, and the -component of the force of particle on .
Recent studies of the present model in the glassy state showed that the system may exhibit shear-localization if the shear rate is imposed by using a conventional Couette cell with moving atomistic walls Varnik-Bocquet-Barrat-Berthier::PRL2003 (). A shear banding is, however, undesired in the context of present analysis since we are interested in the effects of a homogeneous shear. On the other hand, simulations of the present model using the so-called Lees-Edwards boundary conditions along with the SLLOD equations of motion first proposed by Evans and coworkers Evans-Morriss::StatMechNonEqLiq1990 () show that a linear velocity profile forms across the system thus leading to a spatially constant velocity gradient Berthier-Barrat::JChemPhys2002 (). We, therefore, also used this approach for our simulations (see Fig. 1 for an illustration of the Lees-Edwards boundary condition). Within this simulation method, we do indeed observe a spatially constant velocity gradient in all studied cases (Fig. 2).
The external shear does work on the system. Therefore, in order to keep the system temperature at a prescribed value, this extra heat must be removed. This is done by coupling the -component of particle velocities to the so-called Nosé-Hoover thermostat Nose::JCP81 (); Hoover::PhysRevA31 (). We impose a flow in the -direction with a shear gradient in the -direction. The -component of particle velocities are therefore not affected by the local flow velocity, thus simplifying the computation of instantaneous kinetic energy which enters the thermostating part of equations of motion. This choice also ensures that possible artefact due to a profile biased thermostat Evans-Morriss::PRL56-1986 () are absent. Further details about the simulation and the model can be found in reference Varnik-Bocquet-Barrat::JChemPhys2004 ().
The above choice of the simulation method leads to the following set of equations of motion,
Here , and denote coordinates, momenta and the mass of -th particle respectively. is the total number of particles, the Boltzmann constant and the total force on a particle. It is noted that the momenta occurring in the above equations are the peculiar ones, i.e., they correspond to the particles’ momenta in a (local) frame of reference moving with the flow () Evans-Morriss::StatMechNonEqLiq1990 (). The variable, , plays the role of a friction (acceleration) coefficient, since a positive (negative) tends to decrease (increase) . The variation of , on the other hand, is controlled by the deviation of the actual kinetic energy of the system from that prescribed by the temperature, , of the heat bath. The parameter controls the strength of the coupling of the particles’ momenta to the heat bath: the smaller Q the stronger the coupling (see reference Varnik::Dissertation::Mainz2000 () for more details on thermostating methods).
Iii Structural relaxation in a driven system
It is one of the aims of the present section to show that the structural relaxation of an amorphous system exposed to a shear already exhibits signature of the quiescent state of the system. In other words, a study of the structural relaxation in the driven case allows one to find out whether the quiescent state of the system belongs to a liquid-like (supercooled) state or to a glassy phase.
As will be shown below, this possibility of distinguishing between a sheared glass and a driven supercooled liquid is closely related to the fact that, in the supercooled state, the structural relaxation is affected by the external drive at high shear rates only. As the shear rate tends to zero, the (shear-independent) equilibrium relaxation is recovered. In contrast to this, the structural relaxation of a sheared glass is always driven by the imposed shear, since, in an amorphous solid, there is no equilibrium relaxation in the available experimental (or simulation) time window.
obtained from molecular dynamics simulations at various shear rates and temperatures. The use of -component naturally eliminates undesired contributions to particle rearrangements arising, for example, from the motion of particles with the flow and the so-called Taylor dispersion (affecting the displacement of particles along the flow direction as they move in the direction of shear gradient) Taylor-ProcRSocLondonA-219::1953 ().
We measure for a wave vector of , corresponding to the principal peak of the static structure factor. thus reflects collective particle rearrangements on the scale of average nearest neighbor distance, the relevant length scale for a study of the cage effect.
In Fig. 3, is depicted for two temperatures from the supercooled regime. In each panel, the incoherent scattering function is shown for various shear rates ranging from a regime, where the system dynamics is accelerated by the imposed shear, to lower where the structural relaxation is independent of shear rate. First note that, at short times, all curves overlap regardless of the imposed shear rate indicating that the short time dynamics is unaffected by the shear. The situation is, however, different when long time behavior of is considered. At high , the final relaxation of strongly depends on the magnitude of the imposed shear rate: The higher the faster the decay of . However, as is progressively decreased, the dependence of on becomes weaker, eventually disappearing for the lowest shear rates shown in the figure.
A comparison of for a shear rate of in both panels of Fig. 3 shows that a shear rate which does not affect the equilibrium relaxation at a high temperature (here, ) may induce an acceleration of the system dynamics at a lower temperature (). Interestingly, as will be shown in the next section (Fig. 8), the stress response exhibits a linear (non linear) behavior for roughly the same shear rates where is unaffected (affected) by shear.
The above comparison follows a similar plot for in reference Berthier-Barrat::JChemPhys2002 () (see the upper left panel in Fig. 8). The new aspect in our plots is the study of considerably wider range of shear rates allowing to decrease sufficiently in order to recover equilibrium behavior at a temperature as low as , a feature suggested but not shown in reference Berthier-Barrat::JChemPhys2002 ().
The discussion of Fig. 3 underlines the idea that the linear or non linear nature of the system response to an imposed shear rate is the outcome of the competition between the time scale imposed by the flow and the inherent (structural) relaxation time of the system. The faster the inherent system dynamics, the larger the range of shear rates in the linear regime Berthier-Barrat-Kurchan::PRE (); Fuchs-Cates::FaradayDiscuss-2003 (); Miyazaki-Reichman-PRE66 ().
Figure 4 illustrates the effect of shear on the relaxation behavior of a glass. For this purpose, is depicted for a temperature far below the mode coupling critical temperature of the model. Again, at short times, all curves are identical irrespective of the imposed shear rate. As for the long time behavior, exhibits a dependence on at all simulated shear rates. The plot in Fig. 4 is also motivated by the upper right panel in Fig. 8 of reference Berthier-Barrat::JChemPhys2002 (). In addition to a wider range of shear rates shown in our plot, Fig. 4 also compares the sheared to its quiescent counterpart, , measured after a waiting time of (and averaged over independent initial samples). This allows to unambiguously demonstrate the fact that, even at the lowest studied shear rate (), the final decay of the relaxation function at is indeed a result of imposed shear.
In the linear response regime, transport coefficients can often be expressed as time integrals of equilibrium correlation functions. The shear viscosity, for example, is given by the well-known Green-Kubo relation , where ( denotes the system volume and a component of the stress tensor with zero mean, ). Thus, in a range where becomes independent of the applied shear rate, the shear viscosity must also be a constant (Newtonian behavior).
As can be seen from Eq. (1), the shear stress can be expressed as a function of particle positions and (much faster variables) momenta and forces, whereas depends on the particle positions only. Thus, at least for wave vectors corresponding to the average interparticle distance, one expects a faster decay of stress autocorrelation function compared to a relaxation of . This expectation is born out in Fig. 5. However, despite the fact that the time scale for the final decay of is roughly by one order of magnitude shorter than the decay time of (for being the maximum of the structure factor), the cross over from the shear thinning behavior to equilibrium relaxation occurs at practically the same for the both types of correlation functions (compare the data for in Figs. 3 and 5).
This observation underlines the relevance of as an appropriate quantity for at least qualitative studies of non-linear rheology. In a -range where becomes independent of the shear rate, the shear viscosity is also expected to be independent of . The study of the flow curves presented below (Fig. 8) confirms this expectation.
We examine the validity of the so-called ”time-shear superposition principle”, predicted both within the spin glass theory Berthier-Barrat-Kurchan::PRE () and the non-equilibrium MCT Fuchs-Cates::FaradayDiscuss-2003 (). This property indicates that the shape of at sufficiently large times is independent of shear rate. More precisely, , where is a dimensionless time. The shear rate dependent time, , characterizes the time scale of the final decay of the correlator (the so-called -relaxation).
This property is checked for in Fig. 6. In this figure, the same data as shown in Fig. 3 are depicted versus rescaled time ( is estimated from ). As seen from the both panels of Fig. 6, there is a group of curves following a master curve. These curves belong to low , whereas at higher the time-shear superposition principle is violated. This is reminiscent of similar deviations observed in the non-driven supercooled state where the role of is played by temperature Kob-Andersen::PRE52 (); Bennemann-Paul-Binder-Duenweg::PRE57 ().
For the both temperatures shown in Fig. 6, there is an intermediate time window (“intermediate” in the sense that it is long compared to the initial decay time but short compared to the final decay time of ), known as MCT-relaxation regime, for which the master curve can be well fitted to the so-called von Schweidler law,
Here, the non ergodicity parameter, , determines the height of the plateau (to which the short time decay leads) and is the so-called von Schweidler exponent. We find practically the same numerical value for and very close values for at both temperatures. This is consistent with ideal MCT, which yields a theoretical justification of this empirical law Goetze::LesHouches::1989 (). Indeed, here we use the results of equilibrium MCT, since at low shear rates (where the theory is supposed to apply) recovers its equilibrium behavior.
For the glass, a similar plot as in Fig. 6 is shown in Fig. 7. Obviously, also in the glass, the low shear rate data confirm the validity of the time-shear superposition principle. It is worth mentioning that, Eq. (3) remains valid also in the case of a driven amorphous solid. However, in contrast to the case of the quiescent system, where the von Schweidler exponent, , may depend on the system properties at equilibrium, the non-equilibrium MCT predicts for the case of a driven glass Fuchs-Cates::PRL-2002 (); Fuchs-Cates::FaradayDiscuss-2003 (). Using , a fit to Eq. (3) is also shown in Fig. 7. As also expected from the enhanced solid-like character of the system at compared to and , the value obtained for the non ergodicity parameter, , is significantly larger than , obtained in the supercooled state (Fig. 6).
Iv Rheological response
In this section, the rheological response of the system to a spatially constant (homogeneous) shear rate is studied. It will be shown that the stress response in the supercooled state is qualitatively different from that in the glass. In the supercooled regime, the non linear response at high shear rates is always followed by a linear relation between the shear stress and the shear rate (resulting in a constant viscosity) as the shear rate is sufficiently lowered. However, as temperature is decreased, the linear regime is shifted towards progressively lower shear rates and eventually disappears in the glassy phase. In other words, in the glass, any shear rate, as low as it might be, leads to a non linear response.
In order to elucidate this aspect, we investigate the behavior of the flow curves (shear stress versus shear rate) across an ideal glass transition.
Simulation results on the stress response to a time independent shear (flow-curves) are shown in the top panel of Fig. 8 for temperatures ranging from the supercooled state to the glassy phase. In the supercooled state, the shear thinning behavior at high shear rates crosses over to the linear response as . However, as the system is cooled toward the glassy phase, the linear response is reached at progressively lower shear rates. In the glass, on the other hand, the shear stress becomes independent of at low shear rates. For many practical purposes (see e.g. Fig. 9) this stress plateau plays the role of a yield stress: If the system is subject to stresses below this plateau stress, it behaves like a solid body (zero velocity gradient) while for higher stresses it exhibits a liquid-like character. Let us mention that a stress plateau is also observed in recent experiments on the rheology of dense colloidal dispersions Petekidis2004 (); Fuchs2005 ().
The bottom panel of Fig. 8 depicts the same data as in the top panel in a different way. Here, the apparent shear viscosity, , is plotted against . The presence of a linear response regime in the supercooled state is now illustrated as a constant viscosity line, whereas the apparent -divergence of the viscosity in the glass (for ) is also nicely born out. From this panel, one can easily extract, for temperatures in the supercooled state, the shear rate at which the crossover from linear to non linear response takes place. This gives for and for . Interestingly, these values also mark the beginning of the shear rate dependence of for the corresponding temperatures (Fig. 3).
Let us work out an interesting consequence of the presence of a yield stress on the flow behavior of the system. For this purpose, we simulate a Poiseuille type flow in a planar channel. The study of such a situation is interesting since the stress in a Poiseuille-type flow is zero in the channel center and increases linearly with the distance from it. If the fluid under consideration exhibits (at least within the time window accessible for the simulation) a finite yield stress, one may expect that the fluid portion in a certain region around the channel center (defined by the condition that the stress in this part of the channel be lower than the yield stress) should behave like a solid body while it should flow like a liquid further away from this region.
The above discussed Poiseuille-type simulation is performed in the following way. We first equilibrate a system of size (containing 92880 particles) at a temperature of and then quench it to a temperature of . The system is then aged for a waiting time of in order to suppress the fast aging process which occurs directly after the -quench and allow the system to solidify. After this period of time, two solid walls are introduced by immobilizing all particles whose -coordinate satisfies (note that in general; this gives rise to walls of three particle diameter thickness). A flow is then imposed by applying on each particle a constant force of (in LJ units). In the case of a liquid at equilibrium, this would give rise to a parabolic velocity profile. However, as Fig. 9 clearly demonstrates, the velocity profile rather exhibits a behavior expected for a yield stress fluid: In the central region defined by the condition of (see the stress plateau at in the upper panel of Fig. 8) the velocity profile is flat with a zero gradient while it gradually departs from this constant behavior (shear rate becomming nonzero) beyond this central part of the channel. An interesting consequence of this behavior is that, if for a given driving force (pressure gradient) the channel width is too small to allow the formation of stresses above the system’s yield stress, no flow will occur and the whole system will behave like a solid body.
We report results of large scale molecular dynamics simulations on the structural relaxation and the rheological response of a well established glass forming model, an 80:20 binary Lennard-Jones mixture first introduced by Kob and Andersen in the context of the dynamics of supercooled liquids Kob-Andersen::PRL73::1994 ().
Previous studies of the present model Kob-Andersen::PRL73::1994 (); Kob-Andersen::PRE51 (); Kob-Andersen::PRE52 () showed that the model was suitable for an analysis of many aspects of the so-called mode coupling theory of the glass transition (MCT) Goetze::LesHouches::1989 (); Goetze-Sjoegren::RepProgPhys55 (). In particular, for a total density of , an ideal computer glass transition (in the sense that the relaxation times approximately obey a power law divergence predicted by ideal MCT) was observed at a mode coupling critical temperature of .
The present work is strongly motivated by recent theoretical progress on the field of the rheology of dense amorphous systems. While various theoretical approaches Sollich::PRE-1998 (); Berthier-Barrat-Kurchan::PRE (); Fuchs-Cates::FaradayDiscuss-2003 () provide a coherent description of the rheological response of a supercooled liquid (prediction of shear thinning at high shear rates (or driving force), , followed by linear response as approaches zero), they make different predictions with regard to the stress response of an amorphous solid at low shear rates.
In particular, the issue of dynamic yield stress, , (shear stress in the limit of ) remains controversial. While approaches based on the numerical studies of disordered -spin systems in the thermodynamic limit Berthier-Barrat-Kurchan::PRE () suggest that identically vanishes in the glassy phase, both the semi-phenomenological “soft glassy rheology” model Sollich::PRE-1998 () (SGR, an extension of the trap model Bouchaud::JPhysI-1992 () taking into account the presence of an external drive) as well as the idealized version of the non-equilibrium MCT proposed in reference Fuchs-Cates::PRL-2002 () predict the existence of a finite dynamic yield stress in the glassy phase.
Nevertheless, qualitatively similar statements on the relaxation behavior of the correlation functions in a driven disordered solid are made both within the -spin model and in the non-equilibrium MCT. In particular, shear melting of a glass and the recovery of the time translation invariance (shear stops aging Kurchan::condmat1998-2000 ()) is commonly predicted by both these approaches.
An explanation for the absence of dynamic yield stress in the -spin model is given in reference Berthier-Barrat-Kurchan::PRE () relating this property to the presence of infinite free energy barriers occurring in the glassy phase in the thermodynamic limit (). This picture has been checked by studying via Monte Carlo simulations a -spin model with a finite number of spins, whereby considering a case with finite free energy barriers. Indeed, a yield stress in the glass is observed within these studies Berthier2003 ().
The effect of shear on the relaxation behavior of the intermediate scattering function, , as well as the stress autocorrelation function, , is studied both in the supercooled state and in the glass. In the supercooled state, high shear rates lead to an accelerated decay of the correlation function (shear thinning) compared to the equilibrium relaxation which is recovered at sufficiently low shear rates.
Interestingly, despite the fact that decays much faster than (for corresponding to the inverse of the nearest neighbor distance), the cross over shear rate from shear thinning to equilibrium relaxation is practically the same for the both types of the correlation functions and (compare the data for in Figs. 3 and 5). Furthermore, this cross over shear rate also marks the change from linear to non-linear behavior in the stress response (Fig. 8). These observations underline the role of the intermediate scattering function as an appropriate observable for a study of at least qualitative features of non-linear rheology.
In the glass, on the other hand, our simulations clearly indicate that the dynamic yielding (shear melting) of an amorphous solid (Fig. 7) is accompanied by the presence of a finite stress plateau (Fig. 8), which we identify as the yield stress of the model given the restriction to the accessible observation time window. It is worth mentioning that a stress plateau is also found in recent experiments on the rheology of dense colloidal dispersions Petekidis2004 (); Fuchs2005 ().
An interesting consequence of the presence of a yield stress is discussed for the case of a flow driven by an external body force (equivalent of a pressure driven Poiseuille-type channel flow; see Fig. 9). This choice is motivated by the fact that in such a situation the stress in the channel center is zero and linearly increases with distance from the middle of the channel. One thus expects a solid-like behavior in a central region (defined by ; denoting the transverse coordinate measured from the channel center) and a fluid-like behavior in the region between this central part and the channel walls (). As shown in Fig. 9, this expectation is nicely born out. An interesting consequence of this behavior is that, if for a given driving force (pressure gradient) the channel width is too small to allow the formation of stresses above , the whole system will behave like a solid body sticking to the walls.
We thank M. Fuchs, M.E. Cates and J.-L. Barrat for useful discussions. Many thanks to J. Horbach for careful reading of this manuscript and his instructive comments. During this work, F.V. was supported by the Deutsche Forschungsgemeinschaft (DFG), Grant No VA 205/3-2. Generous grants of simulation time by the ZDV-Mainz, PSMN-Lyon and IDRIS (project No 031668-CP: 9) are also acknowledged.
- X. Qiu, H. Ou-Yang, D. Pine, and P. Chaikin, Phys. Rev. Lett. 61, 2554 (1988).
- G. Taylor, Proc. R. Soc. London A 219, 186 (1953).
- B. Ackerson and N. Clark, Phys. Rev. Lett. 46, 123 (1981).
- B. Ackerson, J. Rheol 34, 553 (1990).
- W. Xue and G. S. Grest, Phys. Rev. Lett. 64, 419 (1990).
- H. M. Laun et al., J. Rheology 36, 743 (1992).
- G. Petekidis et al., Physica A 306, 334 (2002).
- A. M. G. Petekidis and P. Pusey, Phys. Rev. E 66, 051402 (2002).
- G. Petekidis, D. Vlassopoulos, and P. Pusey, Faraday Discuss. 123, 287 (1999).
- P. Hébraud, F. Lequeux, J. Munch, and D. Pine, Phys. Rev. Lett. 78, 4657 (1997).
- T. Phung, J. Brady, and G. Bossis, J. Fluid Mech. 313, 181 (1996).
- L. Berthier, J.-L. Barrat, and J. Kurchan, Phys. Rev. E 61, 5464 (2000).
- P. Sollich, Phys. Rev. E 58, 738 (1998).
- M. Fuchs and M. E. Cates, Phys. Rev. Lett. 89, 248304 (2002).
- M. Fuchs and M. E. Cates, Faraday Discuss. 123, 267 (2003).
- K. Miyazaki and D. Reichman, Phys. Rev. E 66, 050501 (2002).
- K. Miyazaki, D. Reichman, and R. Yamamoto, Phys. Rev. E 70, 011501 (2004).
- D. J. Evans and G. P. Morriss, Statistical Mechanics of Non Equilibrium Liquids (Academic Press, London, 1990).
- L. Berthier, J. Phys.: Condens. Matter 15, S933 (2003).
- F. Varnik and O. Henrich, Phy. Rev. B 73, 174209 (2006).
- L. Berthier and J.-L. Barrat, J. Chem. Phys. 116, 6228 (2002).
- W. Kob and H. Andersen, Phys. Rev. Lett. 73, 1376 (1994).
- F. Starr, M.-C. Bellissent-Funel, and H. Stanley, Phys. Rev. E 60, 1084 (1999).
- S. Sastry, Phys. Rev. Lett. 85, 590 (2000).
- W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (1995).
- W. Kob and H. C. Andersen, Phys. Rev. E 52, 4134 (1995).
- M. Nauroth and W. Kob, Phys. Rev. E 55, 657 (1997).
- J.-L. Barrat and A. Latz, J. Phys. Condens. Matter 2, 4289 (1990).
- S. Nosé, J. Chem. Phys. 81, 511 (1984).
- W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
- F. Varnik, L. Bocquet, and J.-L. Barrat, J. Chem. Phys. 120, 2788 (2004).
- F. Varnik, L. Bocquet, J.-L. Barrat, and L. Berthier, Phys. Rev. Lett. 90, 095702 (2003).
- D. J. Evans and G. P. Morriss, Phys. Rev. Lett. 2172 (1986).
- F. Varnik, Ph. D. thesis, University of Mainz, 2000, available from http://ArchiMeD.uni-mainz.de/pub/2001/0007/.
- C. Bennemann, W. Paul, K. Binder, and B. Dünweg, Phys. Rev. E 57, 843 (1998).
- W. Götze, in Les Houches 1989, Session LI, edited by J. P. Hansen, D. Levesque, and J. Zinn-Justin (North-Holland, Amsterdam, 1989), pp. 287–503.
- G. Petekidis, D. Vlassopoulos, and P. Pusey, J.Phycs: Condens. Matter 16, S3955 (2004).
- M. Fuchs and M. Ballauff, J. Chem. Phys. 122, 094707 (2005).
- F. Varnik and K. Binder, J. Chem. Phys. 117, 6336 (2002).
- W. Götze and L. Sjögren, Rep. Prog. Phys. 55, 241 (1992).
- J. Bouchaud, J. Phys. I 2, 1705 (1992).
- J. Kurchan, Proceedings of ”Jamming and Rheology: constrained dynamics of microscopic and macroscopic scales”, ITP, Santa Barbara, 1997, Report-no LPENSL-th-1498 (see also cond-mat/9812347, cond-mat/0011110).