# [

###### Abstract

In order to understand the flow profiles of complex fluids, a crucial issue concerns the emergence of spatial correlations among plastic rearrangements exhibiting cooperativity flow behaviour at the macroscopic level. In this paper, the rate of plastic events in a Poiseuille flow is experimentally measured on a confined foam in a Hele-Shaw geometry. The correlation with independently measured velocity profiles is quantified. To go beyond a limitation of the experiments, namely the presence of wall friction which complicates the relation between shear stress and shear rate, we compare the experiments with simulations of emulsion droplets based on the lattice-Boltzmann method, which are performed both with, and without, wall friction. Our results indicate a correlation between the localisation length of the velocity profiles and the localisation length of the number of plastic events. Finally, unprecedented results on the distribution of the orientation of plastic events show that there is a non-trivial correlation with the underlying local shear strain. These features, not previously reported for a confined foam, lend further support to the idea that cooperativity mechanisms, originally invoked for concentrated emulsions (Goyon et al. 2008), have parallels in the behaviour of other soft-glassy materials.

Plastic flow of foams and emulsions in a channel]Plastic flow of foams and emulsions in a channel: experiments and simulations

B. Dollet, A. Scagliarini and M. Sbragaglia]B. DOLLET, A.SCAGLIARINI and M.SBRAGAGLIA

## 1 Introduction

Foams and emulsions are dispersions of a fluid phase in a liquid phase, stabilised by surfactants. The dispersed phase is constituted of gas bubbles in foams, and liquid droplets in emulsions. These discrete objects are packed together and jammed, which makes foams and emulsions complex fluids: they exhibit a yield stress below which they do not flow, but deform elastically. Above yield stress, they flow like rheothinning fluids. Rheometric measurements in a Couette cell or in cone–plate geometry have shown that the shear stress and the shear rate obey an empirical Herschel–Bulkley law: , with the plastic viscosity and an exponent generally lower than 1, and often close to 0.5 (Princen & Kiss 1989; Marze, Langevin & Saint-Jalmes 2008; Denkov et al. 2009), with some dependence on the surfactants used (Denkov et al. 2009).

The aforementioned measurements did not give access to the microstructure under flow, and other techniques have been developed to visualise it. In emulsions, confocal microscopy on systems of matched optical index have recently enabled to measure the local structure (Jorjadze, Pontani & Brujić 2013) and the velocity field (Goyon et al. 2008; Goyon, Colin & Bocquet 2010; Mansard, Bocquet & Colin 2014). The latter could also be measured using magnetic resonance imaging (Ovarlez et al. 2013). In foams, index matching is not possible, and the route has been to devise bidimensional (2D) experiments, on either bubble rafts at the surface of a pool of soap solution, with or without a confining top plate (Lauridsen, Chanan & Dennin 2004; Dollet et al. 2005; Wang, Krishan & Dennin 2006; Katgert, Möbius & Van Hecke 2008; Katgert et al. 2010), or on bubble monolayers confined between two plates in a Hele-Shaw cell (Debrégeas, Tabuteau & di Meglio 2001).

Among many interesting features such as shear banding (see e.g. Schall & Van Hecke (2010) for a review), these studies have called the Herschel–Bulkley law found in rheometry into question. Among possible flow configurations, the Poiseuille flow in a straight channel is particularly interesting, since this geometry enforces a linear variation across the channel of the shear stress, which vanishes at the centre and reaches its maximum at the side walls. Together with an evaluation of the shear rate from the measured velocity profile, it gives access to the relation at the local scale. In particular, Goyon et al. (2008); Goyon, Colin & Bocquet (2010) have measured this relation in a series of experiments on emulsions, and they have shown that it did not collapse on a single Herschel–Bulkley law. This deviation from a single flow curve was ascribed to wall effects, more precisely to a nonlocal influence of plastic events happening in the vicinity of the boundaries. The velocity profiles were convincingly fitted by a fluidity model (Goyon et al. 2008; Goyon, Colin & Bocquet 2010). This model, based on a kinetic theory approach (Bocquet, Colin & Ajdari 2009), predicted that the fluidity, defined as , is proportional to the rate of plastic events and follows a nonlocal diffusion equation when it deviates from its bulk value. The range of influence appearing in this equation, called the spatial cooperativity, was shown to be of the order of a few times (typically, five) the size of the elementary microstructural constituent (the drop in the case of emulsions) (Goyon et al. 2008; Goyon, Colin & Bocquet 2010; Geraud, Bocquet & Barentin 2013). This picture was later applied to other complex fluids, such as Carbopol gels (Geraud, Bocquet & Barentin 2013), granular media (Amon et al. 2012; Kamrin & Koval 2012), and foams in a 2D cylindrical Couette geometry (Katgert et al. 2010). The fluidity model agrees with existing experiments, and provides a convenient framework to rationalise the flow of complex fluids. However, at least two points remain unclear and deserve further investigation. The first is the boundary condition at solid walls for fluidity. As a matter of fact, most experimentalists have set it as a free fit parameter, which certainly improves the agreement between the measurements and the predictions from the fluidity model, but does not provide any insight on the role of the walls. Only recently, Mansard, Bocquet & Colin (2014) explored the role surface boundary conditions for the flow of a dense emulsion. They show that both slippage and wall fluidisation depend non-monotonously on the roughness, a behaviour that has been interpreted with a simple model invoking the building of a stratified layer and the activation of plastic events by the surface roughness. These results are interesting and call for further verification in terms of numerical simulations (Benzi et al. 2013; Sbragaglia et al. 2012) and other complex fluids (Katgert et al. 2010). Second, the fluidity parameter has not been yet convincingly related to an independent measure of the local density of plastic events. In experiments, only indirect indications of such a relation have been proposed, based on the correlations of the fluctuations of the shear rate (Jop et al. 2012). Using numerical simulations based on the bubble model (Durian 1997), Mansard et al. (2013) were able to measure independently the fluidity and the density of plastic events, but they show that the two quantities are not proportional; more precisely, the rearrangement rate was found to be a sublinear power (with an exponent 0.4) of the fluidity.

Actually, fluidity models offer a potential explanation for the deviation from a unique relation between stress and strain rate, but they are not the only ones. Another approach has been to develop elasto-viscoplastic models (see e.g. Cheddadi, Saramito & Graner (2012) for a review of them) which, in essence, supplement the viscoplastic Herschel–Bulkley rheology by a description of elasticity. These models are local, but since they treat elastic deformation as an independent variable, they also predict deviations from a single Herschel–Bulkley relation. They have been compared with experiments in Couette flows (Cheddadi, Saramito & Graner 2012), but not for Poiseuille flows.

All these theoretical approaches rely crucially on the modelling of plastic events, and how they affect the elastic stress and the flow. However, although this connection between elasticity, plasticity and flow has been studied in foam flows in complex geometries (Dollet & Graner 2007; Dollet 2010; Cheddadi et al. 2011), there is no existing experimental measurement of the rate of plastic events in a Poiseuille flow. 2D foams are particularly well suited for such a study, because elementary plastic events (so-called T1 events) are well characterised by the neighbour swapping of four bubbles (Figs. 2 and 3) and are accessible by image analysis, more easily that in other soft glassy materials.

In this paper, we provide the first experimental measurements of the rate of plastic events in a Poiseuille flow, on a confined foam in a Hele-Shaw geometry. We show that it is closely related to the independently measured velocity profiles, and that there is still a non-vanishing plastic activity towards the centre of the channel. The study of the spatial distribution in the number of plastic events and the simultaneous analysis of the velocity profiles allows to bridge between the details of the irreversible plastic rearrangements and the corresponding cooperativity flow behaviour at the macroscopic level (Goyon et al. 2008; Goyon, Colin & Bocquet 2010; Geraud, Bocquet & Barentin 2013). We choose to explore this connection by looking at the relationship between the localisation length of the velocity profiles and the localisation length of the number of plastic events. In our experiments, because of wall friction, there is no simple relation between shear stress and shear rate. Therefore, we compare the experiments with simulations of emulsion droplets based on the lattice-Boltzmann method (Sbragaglia et al. 2012), which are performed both with, and without, wall friction. Numerical simulations also offer the possibility to test the robustness of some of the experimental findings versus a change in the viscous ratio between the dispersed phase and the continuous phase, this being set to in all the numerical simulations, whereas in foams; in that sense, the simulations look closer to emulsions. The numerical model possesses two advantages that are rarely present together. From one side, it gives a realistic structure of the emulsion droplets, like for example the Surface Evolver method (Cox & Janiaud 2008; Reinelt & Kraynik 2000; Kern et al. 2004); at the same time, due to the built-in properties, the model gives direct access to equilibrium and out-of-equilibrium stresses (Sbragaglia et al. 2012), including elastic and the viscous contributions. In contrast to other mesoscopic models, such as Durian’s bubble model (Durian 1997), our model naturally incorporates the dissipative mechanisms and the interfacial stresses.

The paper is organised as follows. In Sec. 2, we describe the experimental set-up along with the tools of image analysis for characterisation of the plastic events. In Sec 3, and supplementary material presented in Appendices A and B, we review our computational model based on the lattice Boltzmann models (LBM). The review of the computational model will be accompanied by further benchmark tests on the capability of the model to include crucial properties as disjoining pressure and friction. Results and discussions will be the subject of Sec. 4. The experimentally measured velocity profiles (Sec. 4.1) will be compared with local linear and nonlinear models (Sec. 4.2 and Appendix C). Results of numerical simulations and comparisons with the fluidity model (Goyon et al. 2008; Bocquet, Colin & Ajdari 2009) will be the subject of Sec. 4.3. In Sec. 4.4, we compare the localisations of the velocity profiles and of the rate of plastic events. In Sec. 4.5, we will finally report details on the orientation of the plastic rearrangements in the flowing material. Conclusions will follow in Sec. 5.

## 2 Experimental methods

### 2.1 Setup

We have adapted the setup described in Dollet (2010). The foam flows in a Hele-Shaw cell, made of two horizontal glass plates of length 170 cm and width 32 cm, separated by a gap mm thin enough that the foam is confined as a bubble monolayer (Fig. 1a). Two plastic plates of thickness 2 mm are inserted aside the Hele-Shaw cell, so that the width of the channel is reduced to 10.66 cm (Fig. 1b). These plates have a negligible roughness compared to the bubble size. The channel is connected upstream to a vertical chamber (Fig. 1a) in which a soap solution is fed at a prescribed flow rate thanks to a syringe pump (PHD2000, Harvard Apparatus). Nitrogen is continuously blown through injectors at the bottom of this chamber, producing rather monodisperse bubbles (Fig. 1b). The flow rate in each injector is independently controlled with an electronic flow-rate controller (Brooks). We identify the liquid fraction as the ratio of the liquid flow rate to the total flow rate: , with the gas flow rate. The resulting foam accumulates on top of the chamber, over a vertical distance where it drains, then is pushed through the channel. The transit time through the whole channel is less than 10 minutes in all experiments; we do not observe significant change of bubble size during this time, hence coarsening is negligible. The soap solution is a mixture of sodium lauryl-dioxyethylene sulfate (SLES), cocoamidopropyl betaine (CAPB) and myristic acid (MAc), following the protocol described in Golemanov et al. (2008): we prepare a concentrated solution of 6.6% wt of SLES and 3.4% of CAPB in ultra-pure water, we dissolve 0.4% wt of MAc by continuously stirring and heating at 60C for about one hour, and we dilute 20 times in ultra-pure water. The solution has a surface tension mN/m. The contraction region is lit by a circular neon tube, giving an isotropic and nearly homogeneous illumination over a diameter of about 20 cm. Movies of the foam flow are recorded with a CCD camera at a frame rate of 8 frames per second, with an exposure time of 8 ms. The movies are constituted of 1000 images of pixels. The pressure drop is measured across the observation zone, by a water–water differential manometer connected to two points of the channel (Fig. 1b) through tubes full of water. We have performed five different experiments, a summary of which parameters is provided in Tab. 1.

flow rate | bubble | symbol | ||||||||
---|---|---|---|---|---|---|---|---|---|---|

(ml/min) | (%) | area (mm) | (mm/s) | (mm) | (rad) | (rad) | (kPa/m) | (kPa/m) | ||

27.5 | 4.8 | 2.3 | 3.3 | 0.62 | 1.08 | 0.50 | ||||

52.5 | 4.8 | 4.3 | 3.9 | 0.43 | 1.08 | 0.68 | ||||

102.5 | 4.8 | 8.5 | 8.7 | 0.34 | 2.23 | 0.51 | ||||

152.5 | 4.8 | 12.4 | 6.8 | 0.31 | 2.27 | 0.91 | ||||

160.2 | 16.9 | 12.9 | 4.6 | 0.90 | 0.91 | 0.087 |

### 2.2 Image analysis

To extract the relevant rheological information from the movies, we follow a home-made procedure very similar to that presented in Dollet & Graner (2007) and Dollet (2010); we report to these papers for full details. The velocity field is obtained after averaging of all the displacements of all individual bubbles between consecutive frames (about in total). Averaging is performed along 53 lanes aligned with the flow direction. The T1s are tracked as described in Dollet & Graner (2007). For the four bubbles concerned by a T1, we denote () the vector linking the centres of the two bubbles that lose (come into) contact, and the angle of these vectors with respect to the flow direction, that we can restrict to the interval because the orientations of and are irrelevant (Fig. 2 and Fig. 3), and () the position of the midpoint of the centres of the two bubbles that lose (come into) contact. In our program, the detection of appearing and disappearing contacts is first run independently. As a second step, to identify a T1 and minimise artefacts, we decide that a pair of an appearing and a disappearing contact constitute a single T1 if (i) they are on the same image or if the appearing contact happens in the image next to the disappearing one; the latter condition is necessary, because it happens that transient fourfold vertices are erroneously recognised as artificial small bubbles; (ii) the positions and are closer than a critical distance (that we choose to be of the order of the bubble size, to separate from T1s occurring in the neighborhood); (iii) is larger than a critical angle (we choose ), this condition being necessary because of the apparition of the aforementioned spurious bubbles. By visual inspection on 30 images, we estimate that this procedure leads to an uncertainty of no more than 5% on the number of T1s. We then define the quantity as the position of a T1, and we ascribe this information to the box where this position belongs. We thus compute the scalar field of the frequency of T1s per unit time and area:

where is the area of a box and the duration of a movie. Our T1 detection has two major advantages: (i) it is directly based on the topological rearrangements, contrary to indirect characterisations based on velocity correlations; (ii) it yields an unprecedented statistics, up to individual T1s, which enables to average over the same lanes as for the velocity and to perform quantitative analysis.

Finally, a specific advantage of 2D foams is the possibility to measure directly the elastic stress from image analysis. Neglecting the curvature of the bubble edges, the 2D elastic stress tensor writes (Batchelor 1970; Cantat et al. 2013):

(2.0) |

with the line tension, the areal density of bubble edges, and where the average is computed over all bubble edges . Assuming that the films are invariant between the top and bottom plates, the 3D elastic stress equals: .

## 3 Numerical method

For the numerical simulations, we adopt a dynamic rheological model based on the lattice Boltzmann method (LBM) (Benzi, Succi & Vergassola 1992; Chen & Doolen 1998; Aidun & Clausen 2010). Historically, the main successful applications of LBM in the context of computational fluid dynamics pertain to the weakly compressible Navier–Stokes equations (Benzi, Succi & Vergassola 1992; Chen & Doolen 1998) and models associated with more complex flows involving phase transition/separation (Shan & Chen 1993, 1994; Benzi et al. 2009). In particular, we will make use of a computational model for non-ideal binary fluids, which combines a positive surface tension, promoting the formation of diffuse interfaces, with a positive disjoining pressure, inhibiting droplet (or bubble) coalescence. The model has already been described in several previous works (Benzi et al. 2009, 2010, 2013). In this section we review the method and highlight its essential supramolecular features. The mesoscopic kinetic model considers two fluids and , each described by a discrete kinetic distribution function , measuring the probability of finding a particle of fluid at position and time , with discrete velocity . In other words, the mesoscale particle represents all molecules contained in a unit cell of the lattice. The distribution functions evolve in time under the effect of free-streaming and local two-body collisions, described, for both fluids (), by a relaxation towards a local equilibrium () with a characteristic time :

(3.0) |

Local equilibria are given by a low Mach number expansion of the Maxwellian distribution, namely:

(3.0) |

with a set of weights chosen in such a way to maximise the algebraic degree of precision in the computation of the hydrodynamic fields, while is a characteristic velocity (a constant in the model). Our lattice scheme features nine discrete velocities (Shan et al. 2006), whose details and associated weights are reported in Tab. 3 in Appendix A. Coarse grained hydrodynamical densities are defined for both species as well as a global momentum for the whole binary mixture , with . Non-ideal forces () and a body force term () are introduced with the source term in Eq. (3.0). The non-ideal forces include a variety of interparticle forces, . First, a repulsive () force with strength parameter between the two fluids

(3.0) |

is responsible for phase separation (Benzi et al. 2009). The parameter is a characteristic normalisation parameter, used as a free parameter in the model. The “short” range interaction in Eq. (3.0) is extended up to energy shells (lattice links have been normalised to a characteristic lattice velocity). Furthermore, both fluids are also subject to competing interactions whose role is to provide a mechanism for frustration () for phase separation. In particular, we model short range (nearest neighbor, NN) self-attraction, controlled by strength parameters , ), and “long-range” (next to nearest neighbor, NNN) self-repulsion, governed by strength parameters , )

(3.0) |

with a suitable pseudo-potential function. The pseudo-potential is taken in the form originally suggested by Shan & Chen (1993, 1994)

(3.0) |

The parameter marks the density value above which non-ideal effects come into play. The prefactor in (3.0) is used to ensure that for small densities the pseudopotential is linear in the density . Despite their inherent microscopic simplicity, the above dynamic rules are able to promote a host of non-trivial collective effects (Benzi et al. 2009, 2010). The model gives direct access to the hydrodynamical variables, i.e. density and velocity fields, as well as the local (in time and space) stress tensor in the system, the latter characterised by both the viscous as well as the elastic contributions. When the strength parameter in the phase-separating interactions (3.0) is chosen above a critical value, the model achieves phase separation and promotes the emergence of diffuse interfaces. The use of competing interactions (3.0) is instrumental to achieve a positive disjoining pressure (Colosqui et al. 2013). To quantify the emergence of the surface tension and the disjoining pressure, one has to consider a 1D problem. For a planar 1D interface, developing along , the surface tension is a direct consequence of the pressure tensor developing at the non-ideal interface and is computed as the integral of the mismatch between the normal (N) and tangential (T) components of the pressure tensor. Such surface tension scales as as (Benzi et al. 2009)

(3.0) |

The quantity comes from a proper combination of the coefficients in the competing interactions. For repulsive interactions, () the second integral at the rhs is positive-definite, since . With a proper use of the competing interactions, one can choose , and the first term in the rhs of Eq. (3.0) is negative-definite; consequently, one can decrease the surface tension by simply increasing . The decrease of the surface tension goes together with an increase of the disjoining pressure at the thin film interface. The emergence of a positive disjoining pressure can be controlled in numerical simulations by considering a thin film with two non-ideal flat interfaces, separated by the distance . Following Bergeron (1999), we write the relation for the corresponding tensions

(3.0) |

where is the overall film tension. Similarly to what we have done for the surface tension , the expression for is known in terms of the mismatch between the normal and tangential components of the pressure tensor (Toshev 2008; Derjaguin 1989), , where, in our model, can be computed analytically (Shan 2008; Sbragaglia & Belardinelli 2013). All the detailed expressions for the interaction stress tensor are reported in Appendix A. From the relation it is possible to compute the disjoining pressure: a simple differentiation of permits to determine the first derivative of the disjoining pressure, . This information, supplemented with the boundary condition , allows to completely determine the disjoining pressure of the film (Sbragaglia et al. 2012). In Fig. 4 we analyse quantitatively some of these features. In particular we consider the interaction parameters , , with chosen in the interval . All numbers are reported in lbu (lattice Boltzmann units). As we can see, by increasing the value of , we enhance the energy barrier at the onset of the film rupture.

The body force in Eq. (3.0) contains the driving due to the imposed (constant) pressure gradient () and a drag force () mimicking the friction between bubbles and confining plastic plates, as in the experimental setup (Fig. 1). Such drag force is taken to be proportional to the velocity vector, as in Janiaud, Weaire & Hutzler (2006), i.e. . Once the droplets are stabilised with a positive disjoining pressure, different packing fractions and polydispersity of the dispersed phase can be achieved. In the numerical simulations presented in the following sections, the fraction of the continuous phase (i.e. the equivalent of the liquid fraction in the foam experiment) is kept approximately equal to . As already stressed in the introduction, the numerical model provides two basic advantages whose combination is not common. On one hand, it provides a realistic structure of the emulsion droplets, like for instance the Surface Evolver method (Cox & Janiaud 2008; Reinelt & Kraynik 2000; Kern et al. 2004); at the same time, due to its built-in properties, the model gives direct access to dissipative mechanisms in thin films. This latter point will be further discussed and detailed in Appendix B. The viscous ratio between the dispersed phase and the continuous phase is kept fixed to (the simulation parameters are summarised in Tab. 2). This choice is dictated by purely numerical reasons, as numerical instabilities emerge when one considers the case of a viscous ratio much smaller or much larger than unity. Nevertheless, we can use this as an advantage in our joint numerical and experimental study, as it offers the possibility to test the robustness of the experimental findings versus a change in the viscous ratio between the dispersed phase and the continuous phase. It is also comforting that the latest version of our GPU code (Bernaschi et al. 2009) allows for the simulation of emulsion droplets and their statistics in a reasonable amount of time. The current version runs on multiple-GPU and, by using a combination of CUDA streams and non-blocking MPI primitives, it is able to overlap completely the computation within the bulk of the domain with the exchange of the boundaries. Most simulations have been carried out on Kepler “Titan” GPUs, featuring 14 Streaming Multiprocessors, with a total of 2688 cores running at 0.88 Ghz and a memory bandwidth exceeding 200 GBytes/sec. Each run, spanning multi-million time steps for every single set of parameters, takes less than hours, to be compared with a running time of about hours on previous generation (Fermi) GPU cards. The speedup with respect to a highly tuned (multi-core) CPU version is above one order of magnitude. To develop a systematic analysis of plastic events, we perform a Voronoi tessellation (using the voro++ libraries (Rycroft et al. 2006)) constructed from the centres of mass of the droplets, a representation which is particularly well suited to capture and visualise plastic events in the form of droplets rearrangements and topological changes, occurring within the material.

RUN | |||||
---|---|---|---|---|---|

lbu | lbu | lbu | lbu | ||

P1 | |||||

P2 | |||||

P3 | |||||

P4 | |||||

P5 | |||||

P6 | |||||

C1 | |||||

C2 |

## 4 Results and discussions

### 4.1 Experimental velocity profiles

The velocity profiles measured for five experiments are shown in Fig. 5. They are quite flat at the centre () of the channel, although not completely flat as would be expected from a Herschel–Bulkley model, and decrease significantly close to the side walls (). They are well fitted by an exponential profile:

(4.0) |

with a set of three fitting parameters: either , and , or , and . We will retain the latter set of parameters, which has a clear physical meaning: is the centreline velocity, is the relative slip, i.e. the ratio of the slip velocity to the centreline velocity. The parameter , that we will henceforth call the velocity localisation length, describes the range of influence of the walls friction on the velocity profile. The values of the best fitting parameters are reported in Tab. 1. Among the four experiments run at constant control parameters except the driving flow rate, the relative slip tends to decrease, and the localisation length to increase, at increasing flow rate, except the experiment at flow rate 102.5 mL/min. The fifth experiment is run at larger liquid fraction that the four other: it shows a larger relative slip, and a smaller localisation length, than the experiment with comparable flow rate.

### 4.2 Comparison of experiments with a local model

#### 4.2.1 Linear model

To provide analytical reference equations for the velocity profiles and place our work in the context of the existing literature, we start by comparing our velocity profiles to local models. We start by a comparison to the model of Janiaud, Weaire & Hutzler (2006). It is appealing, owing to its simplicity, and it has been shown to reproduce well experimental velocity profiles for foam flows in plane Couette geometry (Katgert, Möbius & Van Hecke 2008). The model considers a steady unidimensional flow, where inertia vanishes identically. We also neglect end effects, hence assume that flow is streamwise invariant. Hence, the flow profile writes: , with the streamwise direction and the spanwise one. The streamwise invariance implies a constant pressure drop: with constant. From momentum conservation, , with the foam/wall friction force per unit area. Taking the streamwise component of the equation, we get:

(4.0) |

where is the component of the stress tensor. The model assumes that for the shear stress: with the shear strain and the shear rate, and the yield stress and the yield strain, respectively, a plastic viscosity of the foam, and a function quantifying the variation of the elastic stress with the shear strain. Inserting this model in (4.0) yields:

where is assumed to be proportional to the velocity, and defined as:

(4.0) |

If we neglect the elastic term for simplicity, we get the following ODE for the velocity:

where:

(4.0) |

A first boundary condition comes from the fact that is a symmetry axis, hence is an even function of , and we recover the exponential profile (4.0): , first proposed in Sec. 4.1 as an empirical fit, with the characteristic velocity proportional to the pressure gradient:

(4.0) |

As shown in Sec. 4.1, it turns out that this functional form, with , and as free fitting parameters, reproduces very well the experimental profiles (Fig. 5). However, there is a second boundary condition, coming from a force balance of the foam at the wall:

(4.0) |

A macroscopic, visible signature of the balance (4.0) is the angle between the bubble edges and the side walls in Fig. 1b, see also Dollet & Cantat (2010). Inserting (4.0), (4.0) and (4.0) in (4.0), , hence:

which we insert in (4.0) to get:

(4.0) |

This new functional form, with and as free fitting parameters, does not fit the experiments. Actually, there is a major discrepancy with the relative slip; setting in (4.0), we get:

(4.0) |

since is very close to 1 for all our experiments. This prediction is much higher than the experimental value, except for the wet foam (Fig. 6).

#### 4.2.2 Nonlinear model

A possible reason for the discrepancy lies in the fact that the friction force is nonlinear in velocity, and the viscous stress is nonlinear in shear rate. Following Denkov et al. (2009), the internal viscous stress for a 3D foam equals (film and interface contribution respectively) with the bubble radius, and:

(4.0) |

with the capillary number ( Pa s: bulk viscosity), and with S.I. an empirical constant for SLES/CAPB/MAc foams. Using as orders of magnitude from the experiments mm and s, we get , hence the film term is dominant, although the surface term is not negligible. Keeping only the film term for simplicity, Eq. (4.0) shows that the viscous stress scales sublinearly with the shear rate:

(4.0) |

where the prefactor (primed to distinguish it from the plastic viscosity in the linear law used in Sec. 4.2.1) is:

(4.0) |

For solutions giving rigid interfaces, like SLES/CAPB/MAc (Golemanov et al. 2008), foam/wall friction is quantified by the force per unit area (or equivalently the wall stress) (Denkov et al. 2009):

(4.0) |

with two empirical constants and , another capillary number, and:

For cm/s, the ratio of the second term to the first term in (4.0) is 5, hence we neglect the second term. Eq. (4.0) then shows that the friction force scales sublinearly with the velocity:

(4.0) |

with the following value of the friction constant (primed to distinguish it from its counterpart in the linear law used in Sec. 4.2.1):

(4.0) |

Like in Sec. 4.2.1, see Eq. (4.0), we can construct a characteristic length from and . To do so, it is convenient to replace the exponent 0.47 by in (4.0), recasting the factor in the definition (4.0) of ; this factor is almost constant, and equal to 1, for all our experiments. The characteristic length is then . We compute with all the experimental values of the parameters appearing in (4.0) and (4.0): mm for , and 2.7 mm for . These orders of magnitude are compatible with the experimental values of the localisation lengths (Tab. 1).

The effect of nonlinear friction and viscous stress on the velocity profile has been theoretically considered for Couette flows (Weaire et al. 2008; Weaire, Clancy & Hutzler 2009), but not for Poiseuille flows. Therefore, in Appendix C, we compute analytically the velocity profile using the nonlinear laws (4.0) and (4.0), and we show that the role of these nonlinearities on the velocity is negligible.

Hence, the model of Janiaud, Weaire & Hutzler (2006) is too simple to model wall slip in our experiments, which suggests that the role of elastic stresses is crucial. This is qualitatively supported by the fact that the only experiment for which the local model is quite accurate in predicting the amount of slip is for a wet foam, which stores less elastic energy (Cantat et al. 2013). To further support this idea, we plot the shear component of the elastic stress and the normal elastic stress difference in Fig. 7. The shear elastic stress is indeed about four times weaker for the wet foam than for the four other experiments. For these experiments at given bubble area and liquid fraction, its variation across the channel is as follows: towards the centre of the channel, although with a significant asymmetry for some experiments, there is a zone of quasilinear increase around . The width of this region decreases slightly at increasing flow rate. Outside this region, the elastic shear stress plateaus to a value which does not depend much on the flow rate. Interestingly, there is still some velocity variation, and a significant plastic activity, outside those regions where the shear elastic stress plateaus. Except the experiment for the wet foam, the normal elastic stress difference is always positive, i.e. the bubbles are elongated streamwise, an effect which is clearly visible on Fig. 1b. It tends to increase towards the wall.

Eq. (4.0) expresses the balance between the driving pressure gradient, the foam/wall friction, and the gradient of elastic and viscous stresses. Close to the middle of the channel, the velocity gradient is very weak, hence the viscous stress is negligible, and the gradient of the shear elastic stress is roughly constant (Fig. 7). It is interesting to compare the value of this gradient and the pressure gradient . Their experimental values are reported in Tab. 1; the pressure gradient is always larger than the gradient of shear elastic stress, the missing part being the friction. This is a major difference with Poiseuille experiments in 3D channels (Goyon et al. 2008; Goyon, Colin & Bocquet 2010; Geraud, Bocquet & Barentin 2013), where the friction force is absent. This prevents us from measuring directly the spanwise stress from the pressure gradient, contrary to the aforementioned studies.

### 4.3 Numerical simulations and comparison with the fluidity model

Having shown the inaccuracies of a local model without elasticity to capture our experimental data thanks to the inspection of the boundary condition at the wall, we now want to test the effect of elasticity. Local visco-elastoplastic models could be used (Cheddadi, Saramito & Graner 2012), but it is not straightforward to deduce from them testable predictions. We thus test a kinetic elastoplastic model (Bocquet, Colin & Ajdari 2009), which encompasses the effect of elasticity through the nonlocal relaxation of elastic stress induced by plastic events. It predicts that the rate of plastic events should be proportional to the fluidity , a prediction that we can readily test. However, there is a difficulty in testing this nonlocal model against our experiments. The role of friction is crucial in experiments, whereas the nonlocal model has been set up and tested in its absence, although recent studies have considered the coupled role of friction and nonlocality (Barry, Weaire & Hutzler 2011). Indeed, friction complicates the stress profile across the channel, as discussed in Sec. 4.2, and it is thus not straightforward to extract relevant flow curves from our experiments. Hence, it is interesting to run numerical simulations, where the friction can be set off and tuned at will.

Various sets of numerical simulations have been performed in the parameter space (see Tab. 2 for the numerical values used), where is meant to be the value of made dimensionless with the channel width and viscosity , i.e.

(4.0) |

where the last equality is based on the definition of the localisation length given in (4.0). A flat velocity profile in the bulk is shown by all curves (Fig. 8), including the case with , witnessing the presence of a non trivial bulk rheology (Goyon et al. 2008; Goyon, Colin & Bocquet 2010; Geraud, Bocquet & Barentin 2013). We also report the experimental velocity profile with flow-rate 152.5 ml/min (see Tab. 1), just to show that we are able to tune the friction parameters in the numerics to achieve the same localisation observed in the experiments for which an equivalent friction parameter can be estimated based on the friction constant (4.0) and the plastic viscosity (4.0). A direct comparison of the velocity profiles of experiments and numerical simulations is somehow irrelevant, because of the boundary conditions: as observed in the experimental data of Fig. 5, slippage is found to occur at the surfaces of the experiments, while the numerical simulations are performed by imposing no-slip at the walls. Boundary conditions have an impact on the microscopic dynamics at the wall (Mansard, Bocquet & Colin 2014). The numerical simulations will be, on the other hand, quite useful in validating the picture of the plastic flow (Sec. 4.4) at changing the friction constant, a feature that is freely tunable in the numerics.

The top panel of Fig. 9 indeed provides some indications that the extra friction (i.e. the confining plates) does not seem to dramatically affect the distribution of plastic events. There we plot the rate of plastic rearrangements, normalised by the total number of events, from experiments and numerics (for the three ’s). Data show a moderately good collapse onto each other. At a given driving pressure drop, increasing wall friction results in a decrease of the total number of plastic events . We could estimate the number in the numerical simulations and it is reported in the bottom panel of Fig. 9: for the same simulation time (see caption of Tab. 2) plastic events diminish from to for increasing from to . A similar trend is observed for the centreline velocity which is reported in the inset of the bottom panel of Fig. 9. To make this statement more quantitative, we notice that the overall decrease in the number of plastic events can be well captured by the function

(4.0) |

a scaling behaviour that can be obtained from the expression of the centreline velocity, in (4.0), with (no-slip boundary condition for the numerics). The parameter in equation (4.0) sets the number of plastic events in the limit whereas the argument of the hyperbolic function in Eq. (4.0) is inversely proportional to the localisation length. The choice of equation (4.0) as a fitting function is suggested by the consideration that the total number of plastic events is dominated by events occurring in boundary regions where the shear stress is approximately constant and, hence, the fluidity is basically proportional to the shear rate . Consequently, being the number of events, by definition, equal to the integral of the corresponding rate and since , we can assume that , which equals the centreline velocity. Interestingly, the estimate of that we get from a best fit procedure () is greater than the estimate of based on the friction localisation in (4.0), which would yield . This is an indication that the localisation lengh in the numerical simulations is larger than the localisation length induced by the extra friction . This is not a surprise, because our numerical simulations have already confirmed the presence of a cooperativity length scale (Sbragaglia et al. 2012) without wall friction. This supports the idea that an effective localisation length results from the sum of the friction localisation plus cooperativity length (Barry, Weaire & Hutzler 2011), an issue that we will further explore in Sec. 4.4.

The study of the spatial distribution in the number of plastic events and the simultaneous analysis of the localisation in the velocity profiles, allows to bridge between the “microscopic” details of local irreversible plastic rearrangements and the macroscopic flow. It has been suggested that such a bridge can be established by introducing a cooperativity scale which determines correlations (nonlocal effects) in the flow rheology (Goyon et al. 2008; Goyon, Colin & Bocquet 2010; Geraud, Bocquet & Barentin 2013). The underlying idea is that correlations among plastic events exhibit a complex spatio-temporal scenario: they are correlated at the microscopic level with a corresponding cooperativity flow behaviour at the macroscopic level. The cooperativity scale (Goyon et al. 2008), directly affects the fluidity, which has been claimed to follow a nonlocal diffusion equation where the diffusivity is directly proportional to (Goyon et al. 2008; Bocquet, Colin & Ajdari 2009)

(4.0) |

The quantity is the bulk fluidity, i.e. the value of the fluidity in the absence of spatial cooperativity (). The nonlocal equation (4.0) has been justified (Bocquet, Colin & Ajdari 2009) based on a kinetic model for the elastoplastic dynamics of a jammed material, which takes the form of a nonlocal kinetic equation for the stress distribution function. In the steady state, under the hypothesis of low cooperativity, the model predicts rheological equations which take the form of equation (4.0), plus an equation predicting a proportionality between the fluidity and the rate of plastic events, a prediction that is interesting to test. A connection between the rate of plastic events and the fluidity field is indeed visible in the top panel of Fig. 9. The dashed line indicates , which is the “synthetic” fluidity profile (up to an unessential numerical scaling factor) based on the hyperbolic cosine fit of the velocity profile (see Sec. 4.1) and a linear variation of the shear stress across the channel. Interestingly, a significant plastic activity remains towards the centre of the channel, and it is well correlated to the fluidity field, which remains finite in such regions, whereas the strain-rate goes to zero. Moreover, a closer inspection reveals that the decrease in the number of plastic events is affected by the friction constant , with a steeper decrease associated with the larger . Fig. 9 calls therefore for a deeper understanding with regard to the link between the rate of plastic events and the local flow properties.

### 4.4 Localisation lengths: comparison of plasticity and shear rate

To go further, we choose to explore the connection between the rate of plastic events and the local flow properties, by looking at the relationship between the localisation length of the velocity profiles, , and the localisation length of the number of plastic events, . This connection enables to compare experiments and simulations, despite their different boundary conditions. The localisation length of the velocity profiles is estimated by a hyperbolic cosine function , from which the decay length is extracted (see Sec. 4.1). Simultaneously, the plastic localisation length is computed out of an exponential fit of the symmetrised rate of plastic events across the channel (Top Panel of Fig. 10). Since our numerical simulations have already confirmed the presence of a cooperativity length scale (Sbragaglia et al. 2012) without wall friction, they are good candidates to complement the experimental findings, showing how the spatial distribution of plastic events is affected by a change in the friction . Hence, in the Bottom Panel of Fig. 10 we also look at the localisation in the numerics, by fixing the pressure gradient and changing , something that cannot be easily done in experiments with the data at hand. At fixed pressure gradient, we show the Log-lin plot of the rate of plastic events from simulations with different . The extracted is found to be a decreasing function of .

In Fig. 11 we report a scatter plot of the shear localisation length versus the plastic localisation length for three sets of data: experiments (symbols as in Tab. 1), simulations with fixed pressure drop and various ’s (filled squares) and simulations with fixed and various pressure drops (filled circles). Fig. 11 shows that the two localisation lengths are indeed close to each other. The fact that the values and agree, confirms the picture of the “plastic flow”; it is also compatible with the fact that the rate of plastic events and the fluidity seem to be proportional (Sec. 4.3). More precisely, since the localisation length is much smaller than the channel width in our experiments and simulations, the stress does not vary much across the localisation zone, hence in most of the channel there is a good correlation between the shear rate and the rate of plastic events. Barry, Weaire & Hutzler (2011) have combined the local model presented in Sec. 4.2 with a nonlocal constitutive equation for the fluidity field, in the case of a Couette flow with linear laws for the viscous stress and the friction (as in Sec. 4.2.1). They predicted that the localisation length of the velocity profile is an increasing function of both the cooperativity length , and of the friction length defined by (4.0). Here, this theoretical prediction can be tested for the first time versus our experiments and simulations. Some care is required in doing so, because of the nonlinear nature of friction and viscous stress in experiments and of the difference between Couette and Poiseuille flows. However, we have shown in Appendix C that the effect of nonlinear friction and viscous stress is very weak. Since the localisation length is much smaller than the channel width in our experiments and simulations, the comparison with the Couette predictions is relevant. For this reason, we also repeated some numerical simulations in a Couette flow geometry (see Tab. 2). The associated data nicely collapse on the same master curve, stressing even more the robustness of our findings at changing the load conditions. The simulations with a Poiseuille flow (inset of the Bottom panel of Fig. 9) show that the localisation length is a decreasing function of at fixed pressure gradient. This is compatible with the local model, more precisely with the expression (4.0) of the friction length. However, the latter tends to diverge at vanishing friction, whereas the localisation length remains finite in this limit; this suggests that the model of Barry, Weaire & Hutzler (2011) breaks down at vanishing friction.

In experiments, for the foam of given liquid fraction and bubble area, the localisation length increases with increasing flow rate. The friction length being fixed, our result may suggest that the cooperativity length could increase with the flow rate. However, looking at the signature of individual plastic events in elastic stress (or displacement) redistribution in their surroundings, we could not find a conclusive signature of an increasing range. This question thus remains open and should be further addressed in the future. On the other hand, Tab. 1 shows that at given flow rate (up to 5%), the localisation length is lower for a wet foam than for a dry one. This is in qualitative agreement with the experiments of Goyon, Colin & Bocquet (2010) for emulsions, who showed that the cooperativity length is an increasing function of the packing fraction above the jamming transition. We remark that testing the effect of the cooperativity length is highly nontrivial; whereas it has been shown to increase towards the jamming transition (Bocquet, Colin & Ajdari 2009; Jop et al. 2012; Nicolas & Barrat 2013), the influence of shear rate remains rather unclear. In other experiments on concentrated emulsions (Goyon et al. 2008; Bocquet, Colin & Ajdari 2009), a unique length was found to account for all experimental data for the flow profiles and local flow curves for a given volume fraction of the emulsion, independently of the pressure drop, confinement and surface nature. However, the flow behaviour in the limit of low shear rates is difficult to access experimentally. These conclusions apply equally well in the case of our numerical simulations where, for a given friction constant, we could not find any signature of an increase of the cooperativity length at changing the shear rate. However, also for the numerical simulations, the situations with low shear are very difficult to analyse, as they require a very large statistic. In some sense, the work presented here bypasses the problem of an accurate measurement of the cooperativity length close to the yield stress and directly explores the link between localisation phenomena in the velocity profiles and the micro-dynamics, characterised by the rate of plastic events and their localisation.

### 4.5 Orientation of the plastic events

The importance of plastic rearrangements has been stressed in that the occurrence of these events induces long range correlations within the soft-glassy material. It is also acknowledged (Picard et al. 2004; Schall et al. 2007) that T1s possess a non trivial angular structure with a quadrupolar topology. It seems, then, reasonable to argue that for a full understanding of the way they determine nonlocal effects inside the system, not only the distribution of their locations in space, but also their orientational properties need to be addressed. Therefore, we go further with the description of plastic events, and study their angular statistics from experiments and simulations. More precisely, focusing on the four bubbles involved in a T1, we define as a disappearing link the segment connecting the centres of the two bubbles which were in contact before the event (and which are then far apart), and as an appearing link the connector between the other two bubble centres (see also Sec. 2 and Figs. 2 and 3); we then measure for each event the angle between such links and the flow direction. We have observed that the angles are reversed between both sides of the channel, consistently with the fact that is an axis of symmetry. Therefore, we choose to analyse the statistics of the quantities and (see Fig. 1 for the sign convention of ). We did not observe a significant variation of the distribution of these angles across the channel, hence we analyse the distributions of these angles for all T1s, whatever their location across the channel. Fig. 12 shows the histogram of and for one experiment and one simulation, while the average and standard deviation of these quantities are summarised for all experiments in Tab. 1. This analysis shows that T1s have preferential orientations: is peaked around 0.5 rad, with a small dispersion, and around rad, with a larger dispersion. The average values do not depend significantly on the flow rate. For the wet foam, is larger, and slightly smaller.

We now derive some reference values for these angles from a microstructural analysis. Since our foams are rather monodisperse, it is interesting to use the simple geometrical model of a sheared 2D hexagonal foam (Princen 1983) (see also Khan & Armstrong (1986)). In this model, the unit cell drawn in dashed lines in Fig. 13 is sheared, and the location of the vertices is computed to comply with the equilibrium rule that the three edges meet at equal angles. To account for the finite liquid fraction, the vertices are decorated with Plateau borders which radius is an increasing function of the liquid fraction (Fig. 13, left): , with the area of one hexagon. This structure can be sheared up to the point where two neighbouring Plateau borders meet, which defines the onset of the T1 (Fig. 13, right).

The two angles and can be computed from simple geometry, when the two Plateau borders come into contact (Fig. 13, right). The length of the edge between the two bubbles about to detach is equal to . Now at a given strain , this length equals (Khan & Armstrong 1986): , where is the side length of the undeformed hexagon. Setting in the latter equation yields the strain at which the T1 occurs; is a decreasing function of . Moreover, (Princen 1983), and we compute from geometrical considerations in the right panel of Fig. 13: , and . Qualitatively, these two expressions show that both and are decreasing functions of , hence increasing functions of , hence of the liquid fraction. The extreme values are found for a dry foam at , for which (Princen 1983), and for the jamming transition for which : varies between rad (dry foam) and rad (jamming transition), and between and . Our measured values are indeed in these ranges. The values of the disappearing angles for the four experiments with are compatible (within experimental dispersion) with the dry foam prediction, the latter indicated with a vertical dashed line in Fig. 12. The predicted increase of the angles with liquid fraction is compatible with the experiments for , but not for .

Although the model by Princen (1983) gives useful reference values, it is difficult to make a more quantitative comparison based on liquid fraction, because the distribution of liquid is specific to each system. In simulations, the films between droplets are thick, and contain a significant proportion of the liquid. In experiments, the distribution of water is complex because of the 3D structure of the bubbles; there is relatively more water close to the confining plates than in the midplane in between (Cox & Janiaud 2008). The hexagonal foam model of Princen (1983) is a good approximation of the structure of our experimental foams across the midplane between the top and bottom confining plates, but the liquid fraction across this plane, relevant in the hexagonal model, is significantly lower than the experimental liquid fraction. Moreover, the measurement of the appearing angle is less precise than that of the disappearing one, because the relaxation of the four bubbles after a T1 is fast; hence, the measurements made on the image after the topological rearrangement may not be representative of the configuration at the instant of a T1. This also explains why the dispersion is larger for than for .

## 5 Conclusions

We have reported on the first experimental study measuring the rate of plastic events in Poiseuille flows of foams. Experiments have been supplemented by numerical simulations, capable to capture the realistic foam structure and to incorporate naturally the expected mesoscopic dynamics. We have addressed the relation between T1 distribution and macroscopic rheology and revealed a link between the localisation lengthscale of the velocity profiles and that of plastic events across the channel, confirming the relevance of cooperativity for foams (Katgert et al. 2010). The use of numerical simulations allowed to study in a controlled way (something not easily feasible in the experiments) the effect of wall friction, helping to confirm its role in the emergence of an extra localisation for the velocity profiles, as predicted theoretically (Barry, Weaire & Hutzler 2011). Our study highlighted that the elasticity gives rise to a complex near-to-wall dynamics which calls for focused studies both experimentally (in the spirit of the recent work by Mansard, Bocquet & Colin (2014)) and numerically, and for a more refined theoretical modelling of the boundary conditions. Finally, unprecedented results on the distribution of the orientation of plastic events show — with good agreement between experiments and numerics — that there is a non-trivial correlation with the underlying local shear strain; this suggests that more complex forms for the propagators invoked in theoretical models of soft-glassy materials (Bocquet, Colin & Ajdari 2009) may be needed, with an explicit angular structure, especially in situation of non-homogeneous stress (as it is for Poiseuille flows).

The authors kindly acknowledge funding from the European Research Council under the EU Seventh Framework Programme (FP7/2007-2013) / ERC Grant Agreement no[279004]. We acknowledge computational support from CINECA (IT). MS and AS gratefully acknowledge M. Bernaschi for computational support and F. Bonaccorso for helpful visualisations of the flowing emulsions from the numerical simulations.

## A Pressure tensor in LBM simulations

In this Appendix we provide the technical details for the lattice Boltzmann pressure tensor used in equations (3.0) and (3.0) to compute both the surface tension and the disjoining pressure at the non-ideal interface. Given the mechanical model for the lattice interactions described in (3.0)-(3.0), an exact lattice theory is available (Shan 2008; Sbragaglia & Belardinelli 2013) which allows to connect the interaction forces to the lattice Pressure Tensor. The exact pressure tensor is given by

The term represents an internal contribution to the pressure tensor, while is a contribution coming from the interactions. We can separately write the contributions coming from the repulsive (r) phase separating interactions (see equation (3.0)), and those coming from competing interactions providing a mechanism of frustration (F) (see equation (3.0))

(A 0) |

The contribution coming from the phase separating interactions is (Sbragaglia & Belardinelli 2013)

(A 0) |

while the contributions coming from the frustrating interactions are given by various terms, , , , , , labeled with the number of the “energy shell” (see Tab. 3)

(A 0) |

(A 0) |

(A 0) |

(A 0) |

Phase Separating Interactions eq. (3.0) | Shell | lattice Links |
---|---|---|