## Abstract

We present numerical simulations based on a Boltzmann kinetic model with competing interactions, aimed at characterizing the rheological properties of soft-glassy materials. The lattice kinetic model is shown to reproduce typical signatures of driven soft-glassy flows in confined geometries, such as Herschel–Bulkley rheology, shear banding and hysteresis. This lends further credit to the present lattice kinetic model as a valuable tool for the theoretical/computational investigation of the rheology of driven soft-glassy materials under confinement.

## 1. Introduction

Many applications in modern science, engineering and biology have prompted the recent blossoming of research in the rheology of soft-flowing and non-ergodic materials, such as emulsions, foams and gels [1–8]. The theoretical understanding of such materials raises interesting questions on its own. Indeed, since soft materials share simultaneously many distinctive features of the three basic states of matter (solid, liquid and gas), their quantitative description does not fall within the traditional methods of equilibrium and/or non-equilibrium statistical mechanics. New concepts and theoretical paradigms are required to adjust many non-standard features, such as rheology, long-time relaxation, dynamic disorder, ageing and related phenomena. In particular, under simple shear conditions, some of these complex fluids may separate into bands of widely different viscosities. This phenomenon, known as ‘shear banding’ [9], involves inhomogeneous flows where macroscopic bands with different shear rates or shear stresses coexist in the sample. Although shear banding is attributed to a shear-induced transition from a microscopic organization of the fluid structure to another, it still raises lots of theoretical and experimental challenges [10,11].

Simultaneously, the need for new tools of analysis also emerges for computational studies, which typically involve many interacting space and time scales, the latter case being particularly acute in view of the aforementioned long-time relaxation. In the recent past, we presented a new conceptual/computational scheme for the numerical simulation of soft-flowing materials [12–15]. The scheme is based on a (lattice) Boltzmann (LB) formulation [16–18] for interacting binary fluids [19,20], in which, by a proper combination of short-range attraction and mid-range repulsion (competing self-interactions), an effective form of frustration was encoded within the physics of the binary mixture [12,13]. More specifically, by tuning the above interactions in such a way as to bring the surface tension down to nearly vanishing values, many typical signatures of soft-glassy behaviour, such as long-time relaxation, dynamical arrest, ageing and nonlinear Herschel–Bulkley rheology, have been clearly detected in numerical studies [12–15]. Loosely speaking, our model fluid can be paralleled to a micro-/nano-emulsion, say oil droplets in water with a surfactant additive, featuring very low surface tension. In this paper, we further elaborate on this model by performing numerical simulations of confined flows and looking at their rheological properties, i.e. the relation between the applied shear rate and the developed shear stress by the fluid. By performing numerical simulations over a wide range of shear rates, we detect the emergence of shear-banding effects, i.e. a shear stress plateau, which is shown to correspond to a wide range of stationary shear rates. These results lend further credibility to the Boltzmann kinetic formulation as a valuable tool for the systematic study of the emergence of nonlinear rheological properties from mesoscopic interactions.

## 2. The model: lattice Boltzmann equation with competing self-interactions

In this section, we review the main properties of the LB model used in the numerical simulations. Further details can be found in our recent work [12,13]. The starting point is the lattice transcription of the generalized Boltzmann equation [16–18] for a multi-component fluid with *S* species, as inspired by the work of Shan & Chen [19] and Shan & Doolen [21]:
2.1In the above, *f*_{is}(** r**,

*t*) is the probability density function of finding a particle of species

*s*=1,…,

*S*at site

**and time**

*r**t*, moving along the

*i*th lattice direction, defined by the discrete speeds

*c*_{i}with

*i*=0,…,

*b*(figure 1). For simplicity, the characteristic propagation time lapse Δ

*t*is taken equal to unity (Δ

*t*=1) in the following. The left-hand side of equation (2.1) stands for molecular free-streaming, whereas the right-hand side represents collisions as a simple time relaxation towards the local Maxwellian equilibrium on a time scale

*τ*

_{s}. The local Maxwellian is truncated at second order, an approximation that suffices to recover the correct athermal hydrodynamic balance

with *c*^{2}_{S} the square of the sound speed in the model and *δ*_{ab} the Kronecker delta with *a*,*b* indicating the Cartesian components (repeated indices are summed upon). The ’s are equilibrium weights used to enforce isotropy of the hydrodynamic equations [17,18]. The equilibrium for the *s* species is a function of the local species density
and the common velocity, defined as
The common velocity receives a shift from the force *F*_{s} acting on the *s* species [19,21]. This force may be an external one or the result of intermolecular interactions. The pseudo-potential forces embed the essential features to achieve phase separation (non-ideal equation of state and non-zero surface tension), as well as a mechanism of frustration through competing self-interactions. More specifically, within each species, the forces consist of an attractive (denoted with (att)) component, acting only on the first Brillouin region (belt, for simplicity), and a repulsive (denoted with (rep)) one acting on both belts, whereas the force between species (*X*) is short-ranged and repulsive. In equations:
where
2.2
2.3
and
2.4
In the above, the groups ‘belt 1’ and ‘belt 2’ refer to the first and second Brillouin zones in the lattice and *c*_{i}, *p*_{i},*w*_{i} are the corresponding discrete speeds and associated weights (figure 1 and table 1). The interaction parameter *G*_{ss′}=*G*_{s′s}, *s*′≠*s*, is the cross-coupling between species, *ρ*_{0} is a reference density to be defined shortly and, finally, *r*_{i}=** r**+

*c*_{i}are the displacements along the

*i*th direction. These interactions are sketched in figure 1 for the case of a two component fluid (say species

*A*and

*B*). This model bears similarities to the next-to-nearest-neighbour frustrated lattice spin models [22,23]. However, in our case, a high-lattice connectivity [24,25] is required to ensure compliance with macroscopic non-ideal hydrodynamics, which is at the core of the complex rheology to be discussed in this work. To this purpose, the first belt is discretized with eight speeds, whereas the second with 16, for a total of

*b*=25 connections (including rest particles). All the weights take the values illustrated in table 1.

The present model has already proven capable of reproducing several signatures of soft-glassy behaviour, such as long-time relaxation, ageing and non-Newtonian rheology. As a general remark, the onset of non-trivial behaviour has been found to be associated with the regime of very low surface tension, a physical quantity that the present scheme allows to tune through the coupling parameters *G*’s and *ρ*_{0}. In particular, once the couplings *G*’s are fixed, surface tension can be brought down to nearly vanishing values by increasing the value of *ρ*_{0}, because the repulsive intra-species force, , contributing a positive surface tension, is scaled by a factor (see equation (2.4)). Full details can be found in our previous publications. To date, the aforementioned effects have been explored only in homogeneous, boundary-free geometries. Given the importance of relating to experimental results, it is of great interest to assess whether the above phenomenology can also be reproduced for the case of soft-glassy flows under geometrical confinement. This is precisely the task undertaken in the following section.

## 3. Numerical results

We have simulated a two-dimensional channel flow driven by the upper wall, moving at speed *U*_{w} along the mainstream direction, *x*. The channel measures *L*=256 lattice sites in length and *H*=512 lattice sites in width. The initial conditions correspond to a modulated density, , with *ρ*_{m}=0.612 and zero flow. Periodic boundaries are imposed at the inlet and outlet sections. The coupling parameters are as follows: , , , and *G*_{AB}=0.405. The time relaxation is *τ*=1 for both species, corresponding to a kinematic viscosity *ν*=1/6 in lattice units. Under ordinary (Newtonian) flow conditions, the wall drive generates a linear flow profile *u*_{x}(*y*)=*U*_{w}*y*/*H*, with an associated homogeneous shear stress *σ*_{xy}=*μU*_{w}/*H*, *μ*=*ρν* being the dynamic viscosity of the fluid. This is precisely the situation our model is expected to reproduce in the regime of sufficiently high surface tension.

### (a) High surface tension: Newtonian rheology

In figure 2, we report the average shear stress *σ*_{xy}, as collected over the set of slices at *x*=1,*L*, as a function of the applied shear rate *S*_{w}=*U*_{w}/*H*. The simulations have been performed with *ρ*_{0}=0.70, corresponding to a sizeable value of surface tension (*γ*∼0.05), for which Newtonian behaviour is expected. As one can appreciate, the shear stress–shear rate relation, on a logarithmic scale, is well fitted by an exponential curve, indicating a linear relation.

### (b) Low surface tension: non-Newtonian rheology and shear banding

In order to investigate the emergence of non-Newtonian behaviour, we have repeated the simulations with a higher *ρ*_{0}=0.83, corresponding to a much lower surface tension *γ*∼5×10^{−3}. The corresponding shear stress–shear rate relation is shown in figure 3. Several comments are in order. At low shear rates, *S*_{w}<2×10^{−6}, the curve still follows a linear, Newtonian relation, as in the high surface tension case described previously. Upon increasing the driving shear rate, however, no further increase of the shear stress is observed, up to *S*_{w}=2×10^{−5}. This is the so-called *shear-banding* effect, signalling qualitative rearrangements in the structural response of the flow pattern, to be discussed shortly. By further increasing the driving shear rate, the stress curve regains an exponential (on a logarithmic scale) behaviour, corresponding, however, to a smaller value of the effective viscosity. The emerging picture is that of a yield-stress, shear-thinning fluid, whose effective viscosity can be mapped into a Herschel–Bulkley relation of the form , with *A*=1.48×10^{−6}, *B*=0.000441 and *α*∼0.25. The appearance of a non-zero yield stress, *σ*_{Y}≡*A*, prepares the ground for the emergence of shear banding. Indeed, the fluid can only flow in regions where the stress exceeds *σ*_{Y} (the shear bands), while in the rest of the domain the fluid stands basically still, like a solid.

### (c) Low surface tension: hysteresis

We have also investigated the occurrence of *memory-dependent phenomena*, such as hysteresis. To this purpose, we have repeated a series of simulations taking as an initial condition the steady state at *S*_{w}=9.76×10^{−5}, and then scanning *S*_{w} backwards to decreasing values. As again apparent from figure 3, the corresponding values of the stress do *not* lie on the same curve as before, but trace in fact a lower lying curve. A similar behaviour is observed by repeating the same backward scan, starting from a higher value *S*_{w}=1.5×10^{−4}, although no crossover is observed in this latter case. By and large, all of the non-Newtonian features portrayed in this figure reveal a striking similarity with recent experiments on soft-glassy flows in microchannels [10].

### (d) Spatial velocity profiles

Finally, we inspect the spatial distribution of the average flow velocity, *u*_{x}(*y*), inside the channel, see figure 4. From this figure, it is apparent that in the low surface tension, non-Newtonian regime, the velocity profiles tend to form a steep layer next to the two walls (the shear bands), with a substantial flattening in the bulk region of the flow, where little dissipation takes place. This flattening is responsible for the plateau in the shear rate–shear stress relation. Much like a turbulent flow, under increasing shear rate, the non-Newtonian fluid reacts in such a way as to confine most of the shear, and stress, to an increasingly thinner near-wall layer. This is again in line with recent experimental observations [10].

## 4. Conclusions

Summarizing, we have shown that the two-component LB model with competing interactions is capable of reproducing distinctive features of driven soft-glassy flows in confined geometries, such as non-Newtonian rheology, shear-banding and hysteresis. Simulation data also show striking similarities with recent experimental results on driven soft-glassy flows in Couette geometries. Future work shall be directed to the quantitative comparison with experimental data and more quantitative analysis in terms of plastic and elastic events [26] from the order parameter snapshots underlying the non-trivial rheology.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society