# Profile blunting and flow blockage in a yield stress fluid: A molecular dynamics study

###### Abstract

The flow of a simple glass forming system (a 80:20 binary Lennard-Jones mixture) through a planar channel is studied via molecular dynamics simulations. The flow is driven by an external body force similar to gravity. Previous studies show that the model exhibits both a static [Varnik et al. J. Chem. Phys. 120, 2788 (2004)] and a dynamic [F. Varnik and O. Henrich Phys. Rev. B 73, 174209 (2006)] yield stress in the glassy phase. These observations are corroborated by the present work, where we investigate how the presence of a yield stress may affect the system behavior in a Poiseuille-type flow geometry. In particular, we observe a blunted velocity profile across the channel: A relatively wide region in the channel center flows with a constant velocity (zero shear rate) followed by a non linear change of the shear rate as the walls are approached. The observed velocity gradients are compared to those obtained from the knowledge of the shear stress across the channel and the flow-curves (stress versus shear rate), the latter being determined in our previous simulations of homogeneous shear flow. Furthermore, using the value of the (dynamic) yield stress known from previous simulations, we estimate the threshold body force for a complete arrest of the flow. Indeed, a blockage is observed as the imposed force falls below this threshold value. Small but finite shear rates are observed at stresses above the dynamic but below the static yield stress. We discuss the possible role of the stick-slip like motion for this observation.

###### pacs:

64.70.Pf,05.70.Ln,83.60.Df,83.60.Fg## I Introduction

The so called soft glassy materials Sollich1997 (); Hebraud1998 (); Fielding2000 () exhibit a rich variety of interesting rheological phenomena. When compared to a Newtonian liquid (defined as a liquid whose viscosity, , does not depend on shear rate, ) the viscosity of a soft glassy system shows pronounced dependence on imposed shear rate. To be specific, let us consider disordered colloidal suspensions Laun1992 (); Petekidis2002 (); Petekidis2003 () as a typical example. It is well known that the shear viscosity of dense colloidal dispersions decreases with increasing (shear thinning) if one focuses the attention on low shear rates. In the limit of high shear rates, on the other hand, the shear viscosity starts to increase with (shear thickenning). While shear thinning is commonly attributed to the competition between the time scale imposed by the external flow and the time scale of inherent structural relaxation, shear thickenning phenomenon is rather understood as originating from hydrodynamic effects Cates () (whose contribution to the stress is negligible at low shear rates and high concentrations but increases considerably at high Phung1996 ()).

In this paper, we study a model system whose rheological properties can be rationalized without taking into account hydrodynamic effects Varnik2006 (). Previous studies of the model showed that it exhibits both a static and a dynamic yield stress. The main difference between the static and the dynamic yield stress lies upon the imposed quantity. While the static yield stress, , is measured in experiments upon imposed stress, the dynamic yield stress, , is measured in experiments upon imposed shear. This is so because a soft glassy material does not develop the same rheological response in the both mentioned cases.

To see this, let us consider a planar Couette cell. If a lateral force per unit area (=stress) is imposed to one of the walls of the Couette cell, the system may resist to the imposed stress if the latter is below some threshold value (which generally depends on temperature and density, say). For stresses higher than this threshold value, on the other hand, the system is ’liquidized’ and flows with a linear velocity profile across the channel Varnik2004 (). The static yield stress thus characterizes the response of an initially non-driven amorphous solid.

If instead of the stress an average shear rate is imposed (by, e.g., moving one of the walls with a constant velocity), the occurrence of a flow is unavoidable by construction. However, the flow profile need not be linear in this case and may exhibit strong spatial heterogeneity. In particular, a two-phase scenario may occur: A region of zero shear (’solid-like’ or ’jammed’) coexisting with a sheared (’liquid-like’) region Varnik2003 (). Interestingly, as the wall velocity approaches zero, the shear stress across the system does not converge to zero (as would be the case for a Newtonian fluid) but seems to saturate at a finite value Varnik2003 (); Fuchs2005 (); Varnik2006 (). This limit of the shear stress for vanishing shear rate is usually defined as dynamic yield stress. It follows from the above description that the dynamic yield stress characterizes the response of a shear molten amorphous solid.

It is noteworthy that the static yield stress is found to be higher than its dynamic counterpart Varnik2004 (), a situation reminiscent of the difference between static and dynamic friction (without any claim for a strict analogy). For the present model and at a temperature and density of and our previous studies yield Varnik2006 () and Varnik2004 () (all quantities are given in Lennard-Jones (LJ) units; see below).

In a glass forming system, the static yield stress does in general depend on the system history and the way the stress is imposed. In particular, it depends on the waiting time (the time elapsed between the temperature/density quench and the beginning of stress ramp) as well as on the rate with which the stress is increased. In the present simulations, however, no stress ramp is applied. Rather, after a waiting time of , we instantaneously switch on a constant external force field exerting a body force of on each particle (see also section II). Furthermore, we focus on sufficiently large waiting times so that, within the time window accessible to our simulations, flow profiles are hardly affected by a dependence of the static yield stress upon (see section III for further discussion of this issue).

In contrast to the static yield stress, the above definition of a dynamic yield stress allows one to avoid these complications. This follows from the fact that reflects the system response to a small but finite steady state shear, where time translation invariance is recovered. Taking the limit of vanishing shear rate does not affect this property.

Obviously, the existence of a dynamic yield stress presupposes the existence of a plateau in the flow curve (shear stress versus shear rate) in the limit of low shear rates. However, although our previous studies support the existence of such a plateau in the case of the present binary LJ model, the issue of a dynamic yield stress in soft glassy materials still remains controversial (see e.g. Besseling2007 () for recent experimental studies of this topic in colloidal dispersions). Therefore, it is worth to check whether our model system also exhibits other features which follow from the existence of a yield stress. As will be shown in this report, the answer to this question is affirmative. In particular, we observe non trivial behavior such as profile blunting and flow blockage in a Poiseuille-type flow geometry, features which can consistently be described assuming the existence of a dynamic yield stress.

The paper is organized as follows. After an introduction of the model and the simulation method in the next section, the effect of yield stress on the behavior of the system in a planar channel flow driven by an external body force will be investigated. In particular, it will be shown that the velocity profile exhibits salient features of a two-phase system: A ’solid-like’ central part with zero velocity gradient and two lateral ’liquid-like’ sections between the channel center and the walls. A consequence of this property on the dependence of the mass flow rate upon the imposed force is worked out and compared to the case of the same system in the normal liquid state, where it behaves like a Newtonian liquid with a shear independent viscosity. A summary compiles our results.

## Ii A binary Lennard-Jones mixture

In order to address the above mentioned issue, we study via molecular dynamics simulations a generic glass forming system, consisting of a 80:20 binary mixture of Lennard-Jones particles (whose types we call A and B) at a total density of .

A and B particles interact via with , , , , and . In order to enhance computational efficiency, the potential was truncated at twice the minimum position of the LJ potential, .

The truncation of the LJ potential introduces a discontinuity in the force field, which could be corrected via a smoothing procedure Haile1992 (). However, the present model with the truncated version of the LJ potential has extensively been studied in the literature and has become a benchmark model for the studies of glassy systems. Therefore, and for the purpose of comparison with previous studies Kob1994 (); Kob1995 (); Kob1995a (); Nauroth1997 (); Berthier2002a (); Varnik2003 (); Varnik2004 (), we keep the model as it is. Note also that the use of a truncated LJ potential is not a priori a disadvantage, since we do not seek a comparison with analytic studies of this specific system. Rather, we are interested in generic features of a glass forming model system for which the present binary LJ mixture has indeed become a prototypical example.

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 . In order to test numerical accuracy, we also performed simulation runs using a smaller time step of . Since no deviations were found between the results obtained for and , we chose the larger time step for all the subsequent simulations.

The present model has first been introduced by Kob and Andersen in the context of the dynamics of supercooled liquids Kob1994 (); Kob1995 (); Kob1995a (), who showed that it was suitable for an analysis of many aspects of the mode coupling theory of the glass transition (MCT) Bengtzelius1984 (); Gotze1989 (); Gotze1992 (). 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 Nauroth1997 (). A similar value is also obtained for a binary mixture of soft spheres Barrat1990 ().

The model has also been studied in the context of the rheology of disordered systems Berthier2002a (); Varnik2003 (); Varnik2004 (); Varnik2006 (); Varnik2006b (). In these studies, aspects such as flow heterogeneity Varnik2003 (), structural relaxation under external drive Berthier2002a () and the existence of a static Varnik2004 () and a dynamic Varnik2006 () yield stress were addressed.

An interesting consequence of the presence of a yield stress on the flow behavior of the system is presented here. 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. As will be shown below, this is a consequence of momentum balance equation Todd1995 (); Varnik2000 () and thus independent of the specific flow profile formed across the channel. If the fluid under consideration exhibits a finite yield stress, one may expect that the fluid portion in a certain region around the channel center (where the stress is below yield stress) should behave like a solid body while it should flow like a liquid further away from this region.

We first equilibrate a system of size (containing 92880 particles) at a temperature of . The temperature is then set to () corresponding to a glassy phase. Note that, the velocity distribution adapts itself very fast (within a time of one LJ unit) to the Maxwell distribution corresponding to the new temperature. Particle configurations, on the other hand, keep the memory of the old (high) temperature for far larger times (see e.g. Varnik2004 () for a short discussion of aging in the case of the present model system).

However, as time proceeds, particles gradually rearrange and the system moves towards that part of the configurational space which corresponds to the new (low) temperature. This process is fast in the beginning, where thermodynamic driving forces are the largest, but slows down with time. In fact, as a characteristic feature of a glassy phase, the system never reaches thermal equilibrium. Rather, it keeps evolving towards it endlessly (aging).

As a result of this aging process, time translation invariance is violated and those system properties which would be independent on the measurement time in an equilibrium state (such as structural relaxation time, diffusion coefficient, static structure function, etc.) show a dependence on the time elapsed between the temperature quench and the beginning of the measurement. Furthermore, quantities computed as time averages also show a dependence on the duration of the measurement, thus reflecting the fact that, in an aging system, ensemble and time averages are no longer equivalent.

In particular, in a glassy state, inherent system dynamics slows down upon aging and structural relaxation times grow (ideally) endlessly, eventually exceeding any other time scale in the problem (such as the time scale imposed by external shear). In this interesting limit, the system no longer behaves like a liquid but rather exhibits solid-like properties, as exemplified by the presence of a finite static yield stress. Since we wish to concentrate on this late time behavior, we first run simulations in the quiescent state for a sufficiently large waiting time of LJ units before imposing an external force.

A time of LJ units in our simulations is sufficiently large in the following sense. First, it is large so that the system exhibits a finite, measurable static yield stress Varnik2004 (). Second, it is large enough so that the increase of the static yield stress upon further aging is slow and does not lead to a qualitative change in the flow behavior (see section III for a more detailed discussion of aging effects on the flow behavior).

After this initial period of time, two solid walls are introduced parallel to the -plane by immobilizing all particles whose -coordinate satisfies (walls of three particle diameter thickness). A flow is then imposed along the -axis by applying on each particle a constant force, . This gives rise to a force density of . The force on a particle is the sum of the interaction forces arising from the Lennard-Jones potential and (recall that for ).

For our geometry with planar solid walls, it follows that , where is the streaming velocity (recall that is the flow direction, the neutral (vorticity) direction and the direction of the velocity gradient). In the steady state and in the absence of a pressure gradient (), the momentum continuity equation thus reduces to , which yields , where we used the symmetry of the shear stress with respect to the -plane ().

The external force does work on the system. This work is transformed into heat via viscous dissipation. In the absence of a thermostat, this would lead to a continuous increase of temperature with time. In order to keep the system temperature at a prescribed value, the viscous heat must be removed. For this purpose, we divide the system into parallel layers of thickness and rescale (once every 10 integration steps) the -component of the particle velocities within the layer, so as to impose the desired temperature . More precisely, we first compute the local kinetic energy per particle within a layer centered at . Here, is the mass of a particle, the number of particles in the layer and the sum runs over the particles in the layer only. A scale factor, , is then determined via the requirement that the new velocities satisfy . This gives . Finally, the new velocities are computed via multiplication of with .

Note that such a local treatment is necessary to keep a homogeneous temperature profile when the velocity profile is not linear. This is so because in this case the shear rate, , and hence the rate of heat production, may vary within the channel giving rise to a temperature gradient if only the average temperature in the channel is adjusted.

However, despite this local temperature control, the heat production close to the walls (where the shear rate is very high) is so fast that a temperature increase in this part of the channel can not be fully avoided (recall that viscous heat scales with ). The magnitude of the excess temperature strongly depends on the applied force per particle and is practically negligible for (see also the discussion of Fig. 6).

## Iii Results

Let us first examine possible effects of aging and flow time on the velocity profile. For this purpose, we prepare the system as described above for various choices of the waiting time . Recall that is the time interval between temperature quench and the time at which the external force is switched on. At this instant, we set the clock to zero, thus marking the onset of body force.

While the system would gradually ’solidify’ in the quiescent state, the external force tries to induce a flow and thus tends to ’fluidize’ (rejuvenate) the system. The extent and the rate of this rejuvenation does, however, strongly depends on the magnitude of the stress formed across the channel. As mentioned above, the stress in the present Poiseuille-type geometry is a linear function of the distance from the channel center: . For a given external force, the shear stress thus increases when approaching the walls.

As a consequence, the system flows first in the proximity of the walls, while the center of the channel behaves rather like a solid body. The effect of aging now reflects itself in the onset of the flow. The larger the slower the initial flow behavior. This feature is nicely born out in the left panel of Fig. 2. As a survey of the data corresponding to reveals, the system hardly flows for and (deforming rather like an elastic solid) while for it ’liquidizes’ in the proximity of the walls.

The data shown in the left panel of Fig. 2 also suggest that, for larger waiting times, the external force must be applied during a longer period of time (larger in the figure) in order to remove memory effects. This is seen by a closer look at the data corresponding to . Here, the curves belonging to and are practically identical, while that of shows significantly slower deformation behavior. However, applying the external force during a time of , aging effects disappear also in the case of .

For a given distance from the channel center, the shear stress increases (and hence the time scale imposed by the external force decreases) with increasing the magnitude of the applied force. This is an important aspect since memory effects (such as effect of aging) decay within a typical structural relaxation time, the latter being of the order of the externally imposed time scale in glassy systems (see e.g. Varnik2006b () for a detailed discussion of this issue). In fact, aging effects are expected to be observable as long as this externally imposed time scale is significantly larger than the waiting time. Since this time scale reduces upon increasing , aging effects are expected to also disappear faster. This property is nicely born out in the right panel of Fig. 2, where data similar to the left panel are shown for the choice of a larger external force per particle, . As seen from this panel, already at a time of , the flow behavior of the system is practically independent of , even for a waiting time of .

Unless otherwise stated, the waiting time is for all the results presented below. After the onset of external force, we wait another period of time (of duration ) before starting the data collection. This helps to reduce the above discussed transient effects related to a competition of aging and shear. Note, however, that in a Poiseuille-type flow of a glassy system, transient effects will always be present to some extent. This is closely related to the fact that the shear rate approaches zero as one moves from the walls toward the center of the channel, thus giving rise to progressively large relaxation times. As a consequence, the time necessary to establish steady state ideally diverges in the center of the channel. This behavior is enhanced in a yield stress fluid, where the zero-shear zone has a finite extension comprising a region around the channel center where the local shear stress is below the fluid’s yield stress. The ’jammed’ region corresponds to this part of the system.

In the case of a Newtonian liquid, the above mentioned approach to drive the flow would give rise to a parabolic velocity profile of the form . However, as Figs. 2 and 3 clearly demonstrate, the velocity profile in the case of the present model in the glassy phase exhibits a quite different behavior: In the central region, the velocity profile is flat with a zero gradient while it gradually departs from this constant behavior (shear rate becoming non zero) beyond this central part.

It is interesting to note that similar (blunted) shapes of the velocity profile are also observed in pressure driven flows of both neutrally buoyant suspensions of spheres (with a size of the order of 1mm) Hampton1997 (); Han1999 () as well as red blood cells (biconcave disks of m thickness and m diameter) Tangelder1986 (); Bishop2001a (). In these cases, however, the profile blunting is usually accompanied by a migration of particles from the wall region towards the center of the channel (’wall migration’) a phenomenon, which is absent in the case of present studies (see the right panel of Fig. 3).

Other examples of blunted velocity profiles occur in systems with a nematic order parameter (see e.g. Callaghan2001 () and references therein). In these systems, the non-trivial rheological response is closely related to a variation of the system structure accompanied by a change in the nematic order parameter. However, such an order parameter-related structural change is absent in the case of our glass forming model.

Nevertheless, a qualitative similarity to the case of present simulations may be found when red blood cells are concerned. This similarity rests upon the fact that profile bunting in red blood cells occurs only if a certain amount of aggregation among red blood cells is present (see e.g. Fig. 7 in reference Bishop2001a ()). The aggregation enhances the solid character of the suspension and leads to a higher yield stress, similar to a reduction of temperature in our model.

The width of the “jammed” region can be estimated from a knowledge of the yield stress in the system via which gives . The question arises whether the static or the dynamic yield stress should be used for an estimate of . Results of our simulations suggest that the dynamic yield stress yields a better estimate of the width of the solid-like region in the channel center. This is exemplified in Fig. 3. In this figure, two vertical dashed lines mark the bounds for the solid-like region estimated via dynamic yield stress whereas the bounds denoted by short vertical solid lines are obtained using the static yield stress. As a survey of the velocity gradient reveals, the use of overestimates the width of the solid-like region.

It is interesting to check the origin of the observed finite shear rates at stresses above the dynamic but below the static yield stress. For this purpose, we performed a series of long simulation runs of a smaller system () at constant imposed stress for a temperature of (far below ). All simulations started after an initial aging of the system during a time of . Interestingly, our data reveal the presence of a stick-slip like plastic deformation. This occurs not only at stresses between and but also at stresses below .

This behavior is illustrated in the left panel of Fig. 4. The panel shows the center of mass position of the whole fluid versus time for some selected values of the imposed stress ranging from (below the dynamic yield stress) up to the static yield stress (). The corresponding center of mass velocity is depicted in the right panel of the same figure. As seen from this panel, for stresses below , the contribution of the stick-slip like motion to the flow velocity drops by roughly an order of magnitude. It is noteworthy that a qualitatively similar stick-slip like dynamics has also been observed in our previous simulations under imposed shear Varnik2003 () as contrasted to the present case of imposed stress.

Using the data shown in Fig. 4, we can estimate the contribution of this intermittent motion to the overall shear rate via the relation . Using , this gives and for and respectively. The finite shear rate observed in the case of stresses between and is therefore closely related to the onset of significant stick-slip motion at these stresses.

One interesting consequence of the presence of a yield stress is that, for a given driving force (pressure gradient) a flow blockage may occur if the channel width is too small to allow the formation of stresses above the system’s yield stress. This ’critical’ channel width is simply estimated via . Obviously, a similar situation would also occur at constant channel width via decreasing the applied body force below . Let us estimate this ’threshold’ . At the temperature and density studied here ( and ) the system exhibits a dynamic yield stress of . Using this value as the yield stress along with we obtain for the minimum force per particle required to induce a flow in the system.

This aspect is illustrated in the left panel of Fig. 5, where velocity profiles are shown for various choices of . As expected, the width of the ’jammed’ region increases with decreasing . For , for example, there remains a narrow liquid-like region close to the walls while the rest of the system is in a ’jammed’ state. This observation is made more quantitative in the right panel of the same figure. For this purpose we determine for each the width of a region with a shear rate larger than a small but finite value. We find that the choice is a good choice for an accurate determination of the position of the ’interface’ between the liquid-like and the solid-like regions. As can be seen from the right panel of Fig. 5, results of this analysis obey well the expected relation if is used. The use of static yield stress () would overemphasize the size of the solid-like part of the channel.

The left panel of Fig. 6 illustrates the velocity gradients, , for exactly the same values of the external force per particle as shown in the left panel of Fig. 5. Here, velocity gradients are compared with profiles of shear rates, , estimated from a knowledge of the local stress and the flow curve upon homogeneous shear: For each , the corresponding stress is first evaluated via . In order to estimate the shear rate corresponding to this shear stress, we use the homogeneous flow curves () determined in our previous simulations Varnik2006 (). In these simulations, an algorithm capable of ensuring a constant shear rate across the system was used (see e.g. Varnik2006b () for details of the simulation method).

As the inset of the left panel of Fig. 6 demonstrates, and agree well within a certain region in the channel, whose extension increases upon decreasing the externally imposed force, . This trend is probably related to the variation of the temperature across the system.

Indeed, the right panel of Fig. 6 shows that the system temperature is not constant overall in the channel but slightly deviates from the prescribed value in a region between the channel center and the walls. The extent of the central region with a constant temperature increases at lower while at the same time the magnitude of the excess temperature with respect to the prescribed value decreases. Noting that significant temperature excesses occur only as the shear stress reaches values comparable to twice the dynamic yield stress of the model, it is not surprising that, at such high stresses, even a local thermostat is not able to ensure a constant temperature profile across the system.

In the above, the knowledge of the homogeneous flow curve is used in order to obtain an independent estimate of the local shear rate across the channel. Similarly, one can use the fact that the shear stress is well known in the channel, , along with the knowledge of the velocity gradient in order to obtain the flow curves, i.e. shear stress as a function of shear rate (velocity gradient). For this purpose, we plot in Fig. 7 for each value of the corresponding values of versus . Due to the symmetry around the mid-plane of the channel, the above procedure would yield identical flow curves using the data from the left and right halves of the channel provided that no statistical uncertainty is present. In reality, however, there is a finite statistical scatter. In order to illustrate this fact, we do not average the results but depict the individual flow curves obtained from the analysis of each half of the channel. The so obtained curves are then compared to the result of simulations at homogeneous shear Varnik2006 (). The comparison is restricted to shear stresses, where the relative deviations between the local temperature and the prescribed one are below one percent. With this restriction, the flow curve obtained from the present simulations agree well with that of homogeneous shear.

Finally, Fig. 8 illustrates how the mass flow rate, , varies with applied force per particle. As a survey of the velocity profiles (Fig. 5) already suggests, the mass flow rate rapidly decreases as the ’critical’ force per particle is approached (note the logarithmic scale for the -axis in Fig. 8). In order to highlight the effect of a finite yield stress, the data corresponding to the same model in the normal liquid state (where the yield stress is identically zero) are also shown. For the normal liquid state, it is straightforward to show that . This relation is nicely born out by by the simulated data. We are, however, not aware of an analytic expression for the mass flow rate in a yield stress fluid. Therefore, we fit the data to the simplest polynomial in which describes the simulated data best.

## Iv Summary

In this paper, we report on the flow behavior of a simple glass forming system. The model consists of a 80:20 binary Lennard-Jones mixture first introduced by Kob and Andersen in the context of the dynamics of supercooled systems Kob1994 (); Kob1995 (); Kob1995a (). It is well known for its capability to form a disordered solid at low temperatures/high densities Kob1997 ().

Previous studies of the rheological response of the model suggest the existence of a finite yield stress, , in the glassy phase. In particular, at a temperature of (deep in the glassy phase), a finite static yield stress of was found via simulations upon imposed stress Varnik2004 (). More recent studies of the same model under imposed shear, on the other hand, showed the existence of a stress plateau in the low shear rate limit of the flow curve (stress versus shear rate) in the glassy phase Varnik2006 (). As the present work also underlines, this stress plateau paly the role of a (dynamic) yield stress.

Here, we study a consequence of the presence of a finite yield stress as a flow is induced in a planar channel via the application of an external force. The stress in such a Poiseuille-type flow is a linear function of the distance from the channel center. A two phase behavior may, therefore, occur in the glassy phase provided that the channel width is sufficiently large in order to ensure that stress close to the walls (where it reaches its maximum) is higher than the yield stress of the system. While the system response is expected to be solid-like (zero shear rate) in a central part of the channel (defined as a region where the stress is below the system’s yield stress), it should flow in the ’wings’ delimited by this central solid-like region and the walls. This expectation is born out by our simulations (Fig. 3).

Furthermore, using the velocity gradient across the channel (Fig. 6), we define the width of the ’jammed’ phase as the size of the region with a shear rate of . The accuracy of this estimate is demonstrated in the right panel of Fig. 5, where it is shown that the relation is well satisfied by the data obtained from the above analysis. As to the numerical value of used in the above formula, our simulation results are consistent with a value of (=stress plateau upon imposed shear Varnik2006 ()) while the use of Varnik2004 () seems to overemphasize the size of the solid-like region. Our simulation results also clearly show that, as the stress increases above the dynamic yield stress, the contribution of a stick-slip like motion to the overall shear rate increases significantly (Fig. 4).

As a consistency check, flow curves obtained from the present Poiseuille-type simulations are compared to the result of previous simulations under homogeneous shear Varnik2006 (). With the exception of relatively high driving forces where uncontrollable viscous heat significantly biases the present simulation results in the the vicinity of the walls, good agreement is found between both approaches (Figs. 6 and 7).

Finally, the dependence of the mass flow rate on the externally imposed force is studied. It is shown that the flow fully stops as . Again, simulated data are well described by this formula provided that is used. In particular, the use of leads to for a complete arrest of the flow, in contrast with our simulations where a finite flow is observed for this value of the external force per particle (Fig. 8).

It is interesting to note that similar (blunted) shapes of the velocity profile are also observed in pressure driven flows of both neutrally buoyant suspensions of spheres (with a size of the order of 1mm) Hampton1997 (); Han1999 () as well as red blood cells (biconcave disks of m thickness and m diameter) Tangelder1986 (); Bishop2001a (). In these cases, however, the profile blunting is usually accompanied by a migration of particles from the wall region towards the center of the channel (’wall migration’) a phenomenon, which is absent in the case of present studies (see the right panel of Fig. 3).

Nevertheless, a qualitative similarity to the case of present simulations may be found when red blood cells are concerned. This similarity rests upon the fact that profile bunting in red blood cells occurs only if a certain amount of aggregation among red blood cells is present (see e.g. Fig. 7 in reference Bishop2001a ()). The aggregation gives rise to a finite yield stress, an important feature whose effects on the flow profile are the focus of the present work.

## Acknowledgments

We thank Jean-Louis Barrat, Lyderic Bocquet for fruitful discussions who strongly motivated the present studies. Parts of the MD simulations related to this work have been performed during the stay of F.V. in Laboratoire de Physique de la Matière Condensée et Nanostructure (LPMCN), Université Claude Bernard, Lyon 1. The financial support by the Max-Planck-Initiative “Multiscale Materials Modelling of Condensed Matter” (MMM) and by the DFG Priority Programme Nano-& Microfluidics (SPP1164, project Va 205/3-2) during the preparation of this manuscript is also acknowledged.

## References

- (1) P. Sollich, F. Lequeux, P. Hébraud, and M. Cates, Phys. Rev. Lett. 78, 2020 (1997).
- (2) P. Hébraud and F. Lequeux, Phys. Rev. Lett. 81, 2934 (1998).
- (3) S. M. Fielding, P. Sollich, and M. E. Cates, J. Rheol. 44, 323 (2000).
- (4) H. Laun et al., J. Rheology 36, 743 (1992).
- (5) G. Petekidis, A. Moussaid, and P. Pusey, Phys. Rev. E 66, 051402 (2002).
- (6) G. Petekidis, D. Vlassopoulos, and P. Pusey, Faraday Discuss. 123, 287 (2003).
- (7) M. Cates, C. Holmes, M. Fuchs, and O. Henrich, , [cond-mat/0310579].
- (8) T. Phung, J. Brady, and G. Bossis, J. Fluid Mech. 313, 181 (1996).
- (9) F. Varnik and O. Henrich, Phy. Rev. B 73, 174209 (2006).
- (10) F. Varnik, L. Bocquet, and J.-L. Barrat, J. Chem. Phys. 120, 2788 (2004).
- (11) F. Varnik, L. Bocquet, J.-L. Barrat, and L. Berthier, Phys. Rev. Lett. 90, 095702 (2003).
- (12) M. Fuchs and M. Ballauff, J. Chem. Phys. 122, 094707 (2005).
- (13) R. Besseling, R. E. Weeks, A. B. Schofield, W. C. K. Poon, Phys. Rev. Lett. 99, 028301 (2007).
- (14) J. M. Haile, Molecular Dynamics Simulation (Wiley, New York, 1992).
- (15) W. Kob and H. C. Andersen, Phys. Rev. Lett. 73, 1376 (1994).
- (16) W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (1995).
- (17) W. Kob and H. C. Andersen, Phys. Rev. E 52, 4134 (1995).
- (18) M. Nauroth and W. Kob, Phys. Rev. E 55, 657 (1997).
- (19) L. Berthier and J.-L. Barrat, J. Chem. Phys. 116, 6228 (2002).
- (20) U. Bengtzelius, W. Götze, and A. Sjölander, J. Phys. C 17, 59115 (1984).
- (21) 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.
- (22) W. Götze and L. Sjögren, Rep. Prog. Phys. 55, 241 (1992).
- (23) J.-L. Barrat and A. Latz, J. Phys. Condens. Matter 2, 4289 (1990).
- (24) F. Varnik, J. Chem. Phys. 125, 164514 (2006).
- (25) B. D. Todd, D. J. Evans, and P. J. Daivis, Phys. Rev. E 52, 1627 (1995).
- (26) F. Varnik, J. Baschnagel, and K. Binder, J. Chem. Phys. 113, 4444 (2000).
- (27) R. Hampton et al., J. Rheol. 41, 621 (1997).
- (28) M. Han, C. Kim, M. Kim, and S. Lee, J. Rheol. 43, 1157 (1999).
- (29) G. Tangelder et al., Circ. Res. 59, 505 (1986).
- (30) J. J. Bishop, A. S. Popel, M. Intaglietta, and P. Johnson, Biorheology 38, 263 (2001).
- (31) E. Fischer and P. T. Callaghan, Phys. Rev. E, 64, 011501 (2001).
- (32) F. Varnik and K. Binder, J. Chem. Phys. 117, 6336 (2002).
- (33) W. Kob et al., Phys. Rev. Lett. 79, 2827 (1997).