Emulsions in microfluidic channels with asymmetric boundary conditions and directional surface roughness: stress and rheology

Abstract

The flow of emulsions in confined microfluidic channels is affected by surface roughness. Directional roughness effects have recently been reported in channels with asymmetric boundary conditions featuring a flat wall, and a wall textured with directional roughness, the latter promoting a change in the velocity profiles when the flow direction of emulsions is inverted [D. Filippi et al., Adv. Mater. Technol., 2023, 8, 2201748]. An operative protocol is needed to reconstruct the stress profile inside the channel from velocity data to shed light on the trigger of the directional response. To this aim, we performed lattice Boltzmann numerical simulations of the flow of model emulsions with a minimalist model of directional roughness in two dimensions: a confined microfluidic channel with one flat wall and the other patterned by right-angle triangular-shaped posts. Simulations are essential to develop a protocol based on mechanical arguments to reconstruct stress profiles. Hence, one can analyze data to relate directional effects in velocity profiles to different rheological responses close to the rough walls associated with opposite flow directions. We finally show the universality of this protocol by applying it to other realizations of directional roughness by considering experimental data on emulsions in a microfluidic channel featuring a flat wall and a wall textured by herringbone-shaped roughness.

---

1 Introduction

In contrast to homogeneous fluids, complex structured fluids composed of soft domains, such as droplets in emulsions,¹ bubbles in foams² or blobs in gels,³ possess a non-linear rheological response to an applied stress. This non-linear response results in a fluidization scenario that profoundly differs from that of a Newtonian fluid, comprising heterogeneous flow structures,^4–6 shear-thinning behaviours^7,8 and also the emergence of yield stress rheology for the more concentrated systems.⁹ Furthermore, spatial cooperativity effects, arising from the presence of finite-sized soft domains, have been observed influencing the overall system fluidization.^10–13 This scenario has prompted several investigations in recent years, encompassing fundamental interests at the interface between physics, chemistry, and engineering.^9,14,15 Beyond this fundamental importance, these materials also show a wide variety of applications, from pharmaceutics¹⁶ to oil recovery,¹⁷ from cosmetics¹⁸ to food industry.¹⁹ The presence of interactions with the boundary conditions brings further complexities. The latter is the case – for example – of emulsions or foams flowing in microfluidic channels with rough walls. Indeed, surface roughness can promote droplet/bubble rearrangements close to the rough walls and alter local fluidization. These concepts were evidenced in experiments on concentrated emulsions^20–22 and later confirmed and further investigated in a variety of studies, both experimental and numerical.^23–26 In some of these studies, microfluidic channels with asymmetric boundary conditions have been employed to address the fluidization induced by surface roughness.^23,24 Roughness-induced fluidization hinges on a specific relationship between the wall stress and the local shear rate close to the wall, as pointed out in ref. 20 and 21. A net difference in the slip velocity – wall stress relation close to rough and smooth surfaces is also evidenced in ref. 27, where flow curves for rough and flat surfaces exhibit a difference for small stresses, with the flow curve for flat surfaces setting below the flow curve for rough surfaces. Similar behavior has been observed close to flat surfaces when the wettability of the wall changes, and the adhesion properties of the droplets/particles at the wall are altered.^26,28

The impact of the directional features of the wall texture geometry has not yet received enough attention. Recently, the role played by directional roughness has been reported for the pressure-driven flow of concentrated emulsions within confined microfluidic channels with asymmetric boundary conditions featuring a flat wall and another wall textured by herringbone-shaped roughness.²⁹ By directional roughness, we mean a realization of the surface texture so that the geometry is different with respect to the flow direction. Considerable effects of the directional roughness on the flow properties may be coupled with the rheology of the flowing materials due to the presence of a dispersed phase. To shed light on this interplay, detailed characterization of the rheological response of emulsions close to the rough wall is needed. In experiments, this analysis is hindered by the need for more information on the local stress field inside the microfluidic channel, which is not directly measurable in the presence of asymmetric boundary conditions since only velocity profiles are accessible.^23,30 This work aims to elaborate a validated protocol to reconstruct stress profiles in channels with asymmetric boundary conditions featuring directional roughness, thus opening the possibility to relate the directional effects on velocity profiles to different rheological responses close to the rough walls when the flow direction is inverted. Specifically, we focus on a minimalist model of directional roughness in two dimensions, with the roughness structured as a series of right-angle triangular-shaped posts equally spaced (see Fig. 1(a)) that make the flow direction-dependent (see Fig. 1(b) and (c)). The simplicity of the model allows us to take full control of all parameters of directional roughness. Numerical simulations capture different velocity profiles as well as rheological curves associated with the two flow directions. Importantly, they allow the establishment of a quantitatively validated protocol based on mechanical balance conditions to reconstruct stress profiles for pressure-driven flows in microfluidic channels with asymmetric boundary conditions. The protocol is quite universal and can also be applied to different roughness realizations. To corroborate this point, we show how the protocol can be applied to experimental data with directional roughness studied in ref. 29.

Fig. 1 Simulation setups. Panel (a): a sketch of the directional roughness realized with right-angle triangular-shaped posts. The posts are structured with a height h comparable to the average droplet diameter d and periodically distributed at a fixed distance λ = g + w. The gap g between consecutive obstacles as well as the obstacle width w depends on the obstacle slope, tuned via angle α. All these parameters are kept fixed in simulations. Light domains refer to the dispersed phase inside the continuous one. Panels (b) and (c): sketches describing the numerical setups of a Couette (panel (b)) and a pressure-driven (panel (c)) setup. The microfluidic channel has height H and asymmetric boundary conditions in both setups. Indeed, the bottom wall of the microfluidic channel is structured with a directional surface roughness (detailed in panel (a)), while the top wall is flat. The definition of forward (FWD) and backward (BWD) flows is shown in both setups.

This paper is organized as follows: in section 2, we briefly describe the numerical methodology (section 2.1) and give details of experiments (section 2.2); numerical results obtained in the Couette and pressure-driven setups are discussed in section 3.1; in section 3.2, we report on the protocol to reconstruct stress profiles from velocity data; this protocol is applied to experimental data in section 3.3. Conclusions are finally drawn in section 4.

2 Materials and methods

2.1 Numerical simulations

Numerical simulations were performed via lattice Boltzmann models (LBMs)^31,32 using the open-source GPU-parallelized code TLBfind.³³ This code specializes in simulating concentrated emulsions in microfluidic channels structured with rough walls under shear flows, so we extended it to perform pressure-driven flow simulations. TLBfind simulates two-dimensional fluids made of immiscible non-coalescing droplets via a multi-component LBM, which was observed to be optimal for simulating problems involving emulsions^{13,25,34–37} (see Fig. 1). The desired number of droplets and the droplet diameters are input parameters in TLBfind, making it a flexible tool for investigating the dynamics of almost monodispersed emulsions with different droplet concentrations Φ. Note that the employed numerical method is a diffuse-interface method, with the interface separating the regions of the two components with a finite thickness. For this reason, we compute Φ as the fraction of domain size occupied by the dispersed phase, i.e., , where ρ_d is the dispersed phase density, ρ₀ is a reference density value, A = H × L is the domain size and Θ is the Heaviside step function (cf. ref. 36). Note that this protocol for the computation of droplet concentration does not ensure a one-to-one match with the experimental definition of Φ. In addition, by tuning the interaction between the multi-component fluid and the walls, we impose non-adhering boundary conditions for the droplets. Further details on the method can be found in ref. 33.

In this work, we simulate the flow of emulsions in a microfluidic channel of height H/d ∼ 11 and length L/d ∼ 40, where d is the average droplet diameter, which is kept fixed at varying droplet concentrations. The microfluidic channel presents periodic boundary conditions along the x-direction and asymmetric boundary conditions along the z-direction with a top flat and a bottom rough wall (see Fig. 1(b) and (c)). Then, the roughness on the bottom wall is structured as shown in Fig. 1(a): the posts are right-angle triangular-shaped, having a height h, being equally spaced by period λ = g + w, with a slope that can be varied via the inclination angle α. The latter establishes the free space length g between two consecutive obstacles as well as the obstacle width w. The choice of h, λ, and α is the result of systematic simulations at changing the post parameters probing the optimal roughness shape to highlight the role of directional roughness, corresponding to the case with h/d ∼ 1, λ/d ∼ 10, and α = 30°.

For the sake of simplicity, hereafter, we name “forward flow” (FWD) the one that follows the natural slope of the roughness posts given by angle α. In contrast, we name “backward flow” (BWD) when the emulsion moves towards the vertical side of the posts (see Fig. 1(b) and (c)). We perform a one-to-one comparison between a diluted and a concentrated emulsion, corresponding to a droplet concentration Φ = 0.384 and Φ = 0.629, respectively. The latter emulsion is sufficiently concentrated to show an incipient yield stress but enough fluid to flow in the microfluidic channel under a significant pressure gradient with no evidence of droplet coalescence. Furthermore, we have also been guided by the results of ref. 29 showing an optimum value of Φ to observe the flow gap between FWD and BWD directions. Velocity profiles have been estimated by applying a coarse-graining procedure at the droplet scale on averaged-in-time hydrodynamical velocity profiles. Concerning the rheological experiment in a Couette setup, we move the upper wall with velocity u = u_wx̂, while the lower one is immobile. The corresponding shear rate has been measured from the velocity profiles. For these simulations, we refer to the FWD flow when u_w > 0 and BWD when u_w < 0 (see Fig. 1(b)). To explore a more complex situation with a space-dependent fluid stress profile, we simulate the flow of emulsions driven by a pressure gradient ∇P = ΔP/L along the stream-flow direction x, where ΔP is the pressure difference across the channel (i.e., a pressure-driven setup). Specifically, the pressure gradient is applied to the system as a volume force with amplitude ∇P. In this setup, we identify the FWD flow when the emulsion moves in a left-to-right direction (i.e., ∇P < 0). In contrast, the BWD flow corresponds to the opposite situation (i.e., ∇P > 0) (see Fig. 1(c)). Notice that the physical quantities measured with the numerical simulations are hereafter shown in lattice Boltzmann simulation units (lbu).

Furthermore, the TLBfind code³³ allows for the local measurement of the stress tensor σ_ij. Specifically, for the numerical simulation data analyzed in this paper, we consider the off-diagonal component of the stress tensor σ_xz which is given by the sum of a viscous contribution (i.e., second order moment of fluid LBM populations) and some elastic contributions due to interactions between the two fluid components.^13,31 In more detail, we consider the interaction between the two components, which turns into lattice forces. Thus, each lattice node interacts only with the first nearest neighbors with a specific force. These forces contribute to the local stress, which we compute as a force for the unit area (further details can be found in the supplementary material of ref. 13). Then, after summing up all contributions, for each z coordinate, the stress tensor component σ_xz is averaged in time along the x coordinate. This average is denoted with σ(z). Mechanical balance implies that σ(z) is constant for the Couette setup while it is a linear function of the z coordinate in the pressure-driven setup, with a slope set by the pressure gradient.²⁰

2.2 Experiments

The microfluidic channel is realized by UV photolithography as described in ref. 23, 24 and 38. The channel is a rectangular capillary having width W = 4 mm, length L = 5 cm and height H = 120 μm. This provides a pressure-driven flow between extended parallel plates of area W × L. The bottom wall of the channel is textured with herringbone-shaped rectangular grooves that are textured with micrometer resolution. Referring to Fig. 2(a), the grooves are inclined by θ = 45° with respect to the longitudinal axis of the channel. As depicted in Fig. 2(b), the grooves have a rectangular shape with height h = 2.5 μm, width w = 17 μm, and they are separated from each other by a gap g = w. In parallel, we realize a non-textured microfluidic channel, with smooth walls, having the same size L, W, H as the textured microfluidic channel, where the stress profile can be experimentally obtained by measuring the pressure gradient in the flow direction. The microfluidic channel is filled with an emulsion consisting of silicone oil droplets dispersed in a refractive index matching phase made up of a mixture of glycerin/water 54% wt/wt and stabilized by tetradecyl-trimethyl-ammonium bromide (TTAB) at 1% wt/wt. The procedure for emulsion preparation is described in ref. 14, 20 and 23. The oil droplets have an average size of d ≃ 9 μm with a polydispersity index of 45%, and occupy a volume fraction Φ = 0.75, at which the emulsion behaves as a yield stress fluid.³⁹ The flow of the emulsion is pressure-driven by a microfluidic controller with a resolution of ±5 mbar in the range of 1 to 1000 mbar (MFCS series from Fluigent, France). The velocity profiles are measured using particle tracking velocimetry (PTV) by tracking fluorescent nanoparticles whose size is a ∼ 0.2 μm, dispersed in the glycerine/water mixture.^21–24 The presence of the herringbone roughness introduces an anisotropy in the flow directions. Depending on the pressure gradient, the flows within are labeled as forward (FWD) or backward (BWD) depending on the direction towards or against the herringbone tip. As described in ref. 29, the roughness affects the entire 3D velocity profiles. However, by considering a 35 μm × 7 μm region of interest (ROI) on the tip of the V-shaped grooves, which is thin compared to the channel width W, the emulsion flow can be considered to be probed over a 2D slice placed on the channel axis y = 0 (see Fig. 2(a)). By measuring the velocity profiles of the emulsion flow within the thin, quasi-2D ROI centered around y = 0, we provide a microfluidic device to assess the role of directional roughness and analyze the relation between stress and shear rate.

Fig. 2 Experimental setup. Panel (a): top view sketch (xy-plane) of the directional roughness, textured on the bottom wall of the microfluidic channel with a herringbone design. Dotted lines highlight a representative V-shaped groove tilted at θ with respect to the flow directions. The emulsion flow is measured locally within a nearly 2D region of interest (ROI, red box) placed in the center of the channel cross section (y = 0, horizontal dashed line) and the middle of the channel length. Panel (b): side view (xz-plane) of the microfluidic channel close to the textured wall. The red rectangle is a vertical slice of the ROI. The scale bars are 2500 times shorter than the channel length L and 200 times shorter than the channel width W, allowing for local probing of the velocity profiles. Conversely, they are only six times shorter than the channel height H due to microfluidic confinement in the vertical direction.

3 Results and discussion

3.1 Numerical simulations

First, we explored the emulsion dynamics in a Couette flow induced by the directional roughness shown in Fig. 1(a). In this configuration, we performed rheological characterization of a diluted (Φ = 0.384) and a concentrated (Φ = 0.629) emulsion, comparing the FWD and BWD directional flow (see Fig. 1(b)). The results are reported in Fig. 3, showing that the flow in the two directions leads to a very similar mechanical response if Φ is low enough to enhance a Newtonian behavior (Φ = 0.384, Fig. 3(a)), while the FWD and BWD flow curves differ at large values of Φ and low values of the shear rate, creating a different rheological response (Φ = 0.629, Fig. 3(b)).

Fig. 3 Numerical results for Couette rheological measurements (see Fig. 1(b)): shear stress σ as a function of the shear rate [small gamma, Greek, dot above] for emulsions with different droplet concentrations Φ: a diluted (Φ = 0.384, panel (a)) and a concentrated emulsion (Φ = 0.629, panel (b)). Data for FWD (▶) and BWD (◁) directions are reported.

To delve deeper into the investigation, we perform a systematic analysis in a pressure-driven setup (see Fig. 1(c)), where the emulsions are driven by a pressure gradient ∇P and the shear stress σ(z) across the microfluidic channel is space-dependent. Velocity profiles u_x(z) as a function of the vertical direction z at fixed ∇P for the same two emulsions considered in the Couette setup are shown in Fig. 4(a) and (b). In this condition, we observed that only in the case of a concentrated emulsion (Fig. 4(b)) an evident gap between FWD and BWD velocity profiles appear at fixed ∇P, confirming the role of a directional roughness on the flow behavior of concentrated emulsions. We now want to show that the observed gap in the flow profiles for the pressure-driven setup is rooted in a different rheological response for more concentrated emulsion triggered by the directional roughness effects. With this aim, we investigated the relation between the shear stress profile σ(z) and the local shear rate [small gamma, Greek, dot above](z). In Fig. 4(c) and (d), we show the absolute value of the measured shear stress profiles |σ(z)| corresponding to velocity profiles in Fig. 4(a) and (b), respectively. The stress is a linear function of the z coordinate, as predicted by mechanical arguments. In the Newtonian case (Φ = 0.384), profiles for |σ(z)| overlap and are symmetric with respect to the microfluidic-channel center. However, this overlap and this symmetry are lost in the more concentrated case (Φ = 0.629). The latter evidence arises from the different droplet behavior close to the rough wall: in the FWD direction, droplets flow and slide on the obstacle slope and “jump” beyond the post, thus also triggering a vigorous plastic activity; contrariwise, in the BWD direction, droplets hit a vertical obstacle and remain blocked, undergoing a slowdown and a local deformation which turns into a more significant stress value. A more quantitative analysis is discussed in the ESI.† Furthermore, notice that, despite an evident shift in the stress profiles between FWD and BWD direction, the droplet concentration remains almost uniform across the channel section, with a slightly larger value close to the rough wall in the BWD direction, which is symptomatic of the droplet blockage when entering in contact with the vertical posts (see ESI†). To make progress, we needed to estimate the local shear rate [small gamma, Greek, dot above](z). In the most concentrated case, velocity profiles present some noise due to the simulation resolution, thus making a precise measurement of [small gamma, Greek, dot above](z) difficult. For this reason, a filtering procedure performed with a polynomial function is needed to compute the velocity gradient and extrapolate the corresponding rheological data. The resulting relation between |σ(z)| and the absolute value of the local shear rate |[small gamma, Greek, dot above](z)| is shown in Fig. 4(e) and (f) for the diluted and concentrated emulsion, respectively, confirming the evidence of a different surface rheological response induced by a directional roughness also in a pressure-driven setup. Notice that, since the asymmetric boundary conditions lead to asymmetric velocity and stress profiles with respect to the center of the microfluidic channel and we are interested in investigating the rheological properties of emulsions close to the rough wall, we consider only the values measured in the microfluidic channel region close to the rough wall where σ(z) > 0. Further, we do not consider the contribution of one droplet boundary layer close to z = −H/2 because data for BWD flows are affected by a concavity change in that region.

Fig. 4 Numerical data for pressure-driven flows (see Fig. 1(c)) at fixed pressure gradient ∇P. In all panels, we show data for FWD (▶) and BWD (◁) flow directions. The left panels refer to the diluted emulsion with Φ = 0.384, while the right panels refer to the concentrated emulsion with Φ = 0.629. Panels (a) and (b): measured velocity profiles u_x(z) as a function of the vertical direction z, normalized by the microfluidic channel height H. Panels (c) and (d): absolute value of the corresponding measured shear stress profiles |σ(z)|. Panels (e) and (f): rheological curves obtained by computing the shear rate [small gamma, Greek, dot above] from panels (a) and (b), respectively.

The analysis of data shown in Fig. 4(e) and (f) is possible because we can directly measure the local stress in the simulations. We now want to devise a quantitatively validated protocol for the reconstruction of the stress profile indirectly from the velocity profiles, i.e., in those situations where the stress is not locally measurable but the velocity profile is. The latter may be the case of real experiments of emulsions driven by a pressure gradient in microfluidic channels with directional roughness in the presence of asymmetric boundary conditions at the channel walls.²⁹ With some validated protocol to reconstruct the stress profiles from experimental velocity profiles, we can then consider data from experiments and see if the picture portrayed by Fig. 4(d) and (f) still holds.

3.2 Stress profile in channels with asymmetric boundary conditions: a validated protocol from numerical simulations

The starting point is the mechanical-balance condition ∂_zσ(z) = ∇P, implying that the stress profile across the microfluidic channel is linear, with a slope that is given by the pressure gradient ∇P. The missing information to reconstruct the full profile σ(z) is the stress value in some given location: we choose the stress value on the flat wall σ_flat because if we know the relation between the slip velocity at the flat wall (u_s,flat) and the corresponding wall-stress σ_flat, the problem is closed. To estimate the relation between u_s,flat and σ_flat one can consider a pressure-driven flow in a flat–flat microfluidic channel (no roughness); in this case, the analytical solution for the stress profile predicts a linear profile that is symmetrical with respect to the microfluidic-channel center: the value of σ(z) on one of the two flat walls is σ(z = ±H/2) = σ_flat = ∓(∇P)H/2. Thus, by collecting measurements of u_s,flat for different values of ∇P, we can obtain a relation between u_s,flat and σ_flat. We show the corresponding data for the concentrated emulsion (Φ = 0.629) in a flat–flat channel as orange circles in Fig. 5. These data are in agreement with the experiments of Seth and coworkers²⁸ for concentrated emulsions in a Couette rheometer with flat and non-adhering walls, showing that the slip velocity u_s,flat scales linearly with the wall-stress σ_flat. Notice that no a priori prediction can be stated for estimating the constitutive laws relating wall stress and slip velocity: despite the linear behavior that would be expected for such boundary conditions, we cannot predict at this stage the scaling prefactor. If the stress–velocity relation σ_flat(u_s,flat) is specific to the coupling between the emulsion and the wall, then we expect to find the very same relation if we consider the emulsions in flat–rough microfluidic channels and we look at the stress–velocity relation close to the flat wall. This is a non-trivial point in that it validates the protocol that allows us to assign the stress at the flat wall in a flat–rough microfluidic channel once we know the relation between stress and slip velocity in a flat–flat channel. Numerical data reported in Fig. 5 confirm the expectations for both FWD and BWD flow directions.

Fig. 5 Slip velocity u_s,flat as a function of the shear stress σ_flat measured in the proximity of the flat wall. Data for FWD (▶) and BWD (◁) flow directions are compared with data for the case of a symmetric flat–flat microfluidic channel (). All data collapse on a linear scaling (dashed lines), in agreement with the non-adhering case in ref. 28. In the inset, experimental slip velocity as a function of the corresponding wall shear stress is reported () for the experimental emulsion at Φ = 0.75 driven in a symmetric flat–flat microfluidic channel having the same size W = 4 mm, L = 5 cm, H = 120 μm, at different pressure gradients from ∇P = L⁻¹ 100 mbar to ∇P = L⁻¹ 850 mbar (see section 3.3).

In summary, to measure the stress profile from velocity data in a pressure-driven microfluidic channel with asymmetric (i.e., flat–rough) boundary conditions, we need to follow these steps: (i) we can first estimate the relationship between the slip velocity and the stress on the flat wall σ_flat(u_s,flat) from dedicated experiments in microfluidic channels with symmetric boundary conditions (both flat walls); (ii) the mechanical-balance condition implies a linear stress profile σ(z) as a function of the vertical position z with a slope given by the pressure gradient ∇P; (iii) from the measure of the slip velocity at the flat wall (u_s,flat) in experiments in microfluidic channels with asymmetric boundary conditions, we compute σ_flat(u_s,flat) from the relation obtained in (i), thus determining the stress profile as:

In the next section, we apply the operative protocol outlined above by points (i)–(iii) to experimental data collected in a microfluidic channel textured by a directional roughness consisting of herringbone-shaped grooves.²⁹

3.3 Protocol application to experiments with a different directional roughness

To address the universality of the protocol described in section 3.2 to reconstruct the stress profiles in microfluidic channels with different directional roughness realizations, we considered a microfluidic channel with flat–flat walls separated by the same height H of the flat–rough microfluidic channel shown in Fig. 2. Notice that the geometry of the directional roughness addressed by the experiment differs from that of the numerical simulation. Nevertheless, the protocol that we have numerically validated is based on the mechanical balance hence it is expected to hold true also on different roughness realizations.

In the inset of Fig. 5, we report the slip velocity measured in a flat–flat channel (u_s,flat) as a function of wall stress σ_flat = L⁻¹ΔP(±H/2). Data refer to an emulsion with droplet concentration Φ = 0.75 for different pressure gradients. The values of slip velocities are the average over values measured on both flat walls to avoid noise fluctuations. The resulting wall slip velocity–stress relation is as linear as that obtained from numerical data for a flat–flat channel. The latter data show that the wall slip is proportional to the wall stress on the flat wall. We then apply the protocol described in section 3.2 to our experimental data: we first measured u_s,flat from velocity profiles reported in Fig. 6(a) in the flat–rough microfluidic channel (described in section 2.2) at fixed pressure gradient ∇P = L⁻¹ 700 mbar; then, we estimated the corresponding value of σ_flat from the inset of Fig. 5; finally, we applied the key protocol steps discussed in section 3.2 to obtain the stress profiles σ(z). Fig. 6(b) reports the resulting absolute value of the shear stress profiles |σ(z)| for FWD and BWD flows shown in Fig. 6(a), and it confirms that the existence of different slip velocities on the flat wall results in a shift between the FWD and the BWD shear stress profiles. The experimental flow curves shown in Fig. 6(c) are achieved by plotting |σ(z)| as a function of the absolute value of the shear rate |[small gamma, Greek, dot above](z)| obtained as the derivative of the velocity profile [small gamma, Greek, dot above](z) = ∂_zu(z). Also in experiments, we focus only on the half of the channel close to the rough wall, i.e., the spatial region where σ(z) > 0, since we are interested in investigating the rheological properties close to the rough wall. Flow curves in Fig. 6(c) verify the emergence of different rheological responses in the pressure-driven flow with a different directional texturing. This result further suggests that the directional roughness geometry triggers a different rheological response between FWD and BWD flow directions due to the different directional stress on the droplets. Note that it has been essential to perform a fitting procedure on velocity profiles to obtain experimental flow curves; otherwise, the noise introduced by the numerical derivatives of the raw velocity data is so high that it masks the stress difference of the FWD and BWD flow curves.

Fig. 6 Experimental data. FWD (▶) and BWD (◁) flow directions for the emulsion at Φ = 0.75 in the flat–rough microfluidic channel, driven by ∇P = L⁻¹ 700 mbar. Panel (a): experimental velocity profiles u_x(z) measured by PTV as a function of the vertical direction z, normalized by the height H of the microfluidic channel. Panel (b): absolute value of the stress profile |σ(z)| within the microfluidic channel obtained from eqn (1) by running the wall slip–stress protocol described in section 3.2. Panel (c): flow curves extracted from the stress data reported in (b) as a function of shear rate [small gamma, Greek, dot above] extracted from panel (a).

4 Conclusions

Emulsions flowing in microfluidic channels with rough walls can exhibit directional roughness effects.²⁹ This scenario manifests when the realization of the wall roughness produces different flow profiles depending on the flow direction. When the textures are arranged asymmetrically, a protocol is needed to reconstruct stress profiles inside the microfluidic channel for pressure-driven flows. This protocol is key to studying how the emulsion rheology close to the rough walls depends on the directionality of the flow. In this study, we established this protocol with a minimalist model, featuring lattice Boltzmann numerical simulations of emulsions moving within a 2D microfluidic channel textured at one wall with right-angle triangular shaped posts, using the open-access TLBfind code.³³ We show that directional effects on the velocity field are associated with the onset of a different rheological response in the two flow directions at large droplet concentrations. Numerical simulations are essential to develop an operative protocol to obtain the stress profile inside the microfluidic channel starting from the velocity profile of pressure-driven flows. The robustness of this protocol hinges on the underlying mechanical balance conditions. To show the universal application of the presented protocol to other directional roughness realizations, we applied the protocol to experimental data taken on a microfluidic channel textured by herringbone roughness, considering the emulsion flow at a volume fraction above the yielding on the tip of the V-shaped grooves.²⁹

Notably, as reported in Fig. 3(b), 4(f) and 6(c), at a fixed shear rate, the shear stress is more significant in the backward direction than in the forward direction. This rheological difference echoes previous findings reporting larger shear stress values in a concentrated emulsion when flowing through a rough microfluidic channel compared to a flat one.²¹ It also echoes the observations of a different rheological response when boundary hydrophobicity is changed and an increase in shear stress is observed with a hydrophobic wall compared to a hydrophilic one.^26,28 However, in the present study, the rheological response is not due to a change of the wettability pattern on the wall but rather by a change of flow direction with respect to a fixed suitable physical texture of the wall roughness. Notice that the developed protocol for stress estimation proves an invaluable tool for experimentally estimating the stress within a microchannel with asymmetric boundary conditions on opposite walls.

Resuming, we confirm that boundary conditions consisting of textures asymmetrically arranged with respect to the flow direction can induce different rheological responses resulting in different macroscopic flow properties. As expected, either the bulk rheology or the presence of a dispersed phase with a characteristic size comparable to the roughness plays a crucial role in setting the local stress close to the asymmetric rough elements. This property can be fruitfully exploited to design a convenient geometry to induce a local change in the rheological response of the flows.

Data on directional roughness are scarce in the literature; therefore, there are many perspectives for future advancements. For example, it would be interesting to understand the optimal design of a directional roughness that allows for passive control of the response of the emulsions in microfluidic channels. This obviously depends on the micro-mechanics of the emulsion droplets close to the wall, and in this perspective a synergistic work between experiments and numerical simulations could help to shed light on the matter. Furthermore, detailed characterization of the role played by the pressure gradient in the pressure-driven setup may better clarify the parameter combination that maximizes the directional roughness effect on emulsion flow. Finally, it would be interesting also to investigate the stress–slip relationships as a function of the shape of the roughness in a corresponding symmetric scenario, i.e., with a microfluidic channel patterned with two identical directional roughness. In this case, the general protocol explained in section 3.2 is trivial because the channel symmetry establishes the stress profile to cross zero in the center of the channel, but the velocity profiles and slip velocities may differ because the emulsion interacts with a different roughness realization in the two directions.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work has been carried out within the TEXTAROSSA project (G. A. H2020-JTI-EuroHPC-2019-1 No. 956831). Partial funding from BIRD 2021 “BiodivSeq” of Padua University is kindly acknowledged. D. F. acknowledges funding from Programma Operativo Nazionale (PON) “Ricerca e Innovazione” 2014–2020 (ReactEU-RUAPON-Viganov, No. 19-G-12542-4) of the Italian Ministry of University and Research. MS gratefully acknowledges the support from the National Center for HPC, Big Data and Quantum Computing, Project CN_00000013 – CUP E83C22003230001, Mission 4 Component 2 Investment 1.4, funded by the European Union – NextGenerationEU.

Notes and references

1. S. R. Derkach, Adv. Colloid Interface Sci., 2009, 151, 1–23.
2. X. Ji, X. Wang, Y. Zhang, W. R. Peltier and P. P.-C. Wu, J. Phys.: Condens.Matter, 2005, 17, R1041.
3. D. R. Picout and S. B. Ross-Murphy, Sci. World J., 2003, 3, 105–121.
4. P. Jop, V. Mansard, P. Chaudhuri, L. Bocquet and A. Colin, Phys. Rev. Lett., 2012, 108, 148301.
5. B. Dollet and C. Raufaste, C. R. Phys., 2014, 15, 731–747.
6. B. Dollet and C. Bocher, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, 38, 1–11.
7. R. Pal, AIChE J., 1996, 42, 3181–3190.
8. R. Pal, J. Colloid Interface Sci., 2000, 225, 359–366.
9. D. Bonn, M. M. Denn, L. Berthier, T. Divoux and S. Manneville, Rev. Mod. Phys., 2017, 89, 035005.
10. A. Nicolas and J. L. Barrat, Phys. Rev. Lett., 2013, 110, 138304.
11. V. Mansard, A. Colin, P. Chaudhuri and L. Bocquet, Soft Matter, 2013, 9, 7489–7500.
12. A. Scagliarini, B. Dollet and M. Sbragaglia, Colloids Surf., A, 2015, 473, 133–140.
13. B. Dollet, A. Scagliarini and M. Sbragaglia, J. Fluid Mech., 2015, 766, 556–589.
14. V. Mansard and A. Colin, Soft Matter, 2012, 8, 4025–4043.
15. N. J. Balmforth, I. A. Frigaard and G. Ovarlez, Annu. Rev. Fluid Mech., 2014, 46, 121–146.
16. E. Bouyer, G. Mekhloufi, V. Rosilio, J. L. Grossiord and F. Agnely, Int. J. Pharm., 2012, 436, 359–378.
17. Y. Zhou, D. Yin, W. Chen, B. Liu and X. Zhang, Energy Sci. Eng., 2019, 7, 1046–1058.
18. D. Venkataramani, A. Tsulaia and S. Amin, Adv. Colloid Interface Sci., 2020, 283, 102234.
19. B. Ozturk and D. J. McClements, Curr. Opin. Food Sci., 2016, 7, 1–6.
20. J. Goyon, A. Colin, G. Ovarlez, A. Ajdari and L. Bocquet, Nature, 2008, 454(7200), 84–87.
21. J. Goyon, A. Colin and L. Bocquet, Soft Matter, 2010, 6, 2668–2678.
22. V. Mansard, L. Bocquet and A. Colin, Soft Matter, 2014, 10, 6984–6989.
23. L. Derzsi, D. Filippi, G. Mistura, M. Pierno, M. Lulli, M. Sbragaglia, M. Bernaschi and P. Garstecki, Phys. Rev. E, 2017, 95, 052602.
24. L. Derzsi, D. Filippi, M. Lulli, G. Mistura, M. Bernaschi, P. Garstecki, M. Sbragaglia and M. Pierno, Soft Matter, 2018, 14, 1088–1093.
25. F. Pelusi, M. Sbragaglia, A. Scagliarini, M. Lulli, M. Bernaschi and S. Succi, Europhys. Lett., 2019, 127, 34005.
26. J. Paredes, N. Shahidzadeh and D. Bonn, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 042313.
27. S. P. Meeker, R. T. Bonnecaze and M. Cloitre, Phys. Rev. Lett., 2004, 92, 198302.
28. J. R. Seth, C. Locatelli-Champagne, F. Monti, R. T. Bonnecaze and M. Cloitre, Soft Matter, 2011, 8, 140–148.
29. D. Filippi, L. Derzsi, F. Nalin, A. Vezzani, D. Ferraro, A. Zaltron, G. Mistura and M. Pierno, Adv. Mater. Technol., 2023, 2201748.
30. A.-L. Vayssade, C. Lee, E. Terriac, F. Monti, M. Cloitre and P. Tabeling, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 052309.
31. T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer International Publishing, 2017, pp. 61–104.
32. S. Succi, The Lattice Boltzmann Equation: For Complex States of Flowing Matter, 2018.
33. F. Pelusi, M. Lulli, M. Sbragaglia and M. Bernaschi, Comput. Phys. Commun., 2022, 273, 108259.
34. X. Shan and H. Chen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 1815.
35. R. Benzi, M. Sbragaglia, S. Succi, M. Bernaschi and S. Chibbaro, J. Chem. Phys., 2009, 131, 104903.
36. F. Pelusi, M. Sbragaglia, R. Benzi, A. Scagliarini, M. Bernaschi and S. Succi, Soft Matter, 2021, 17, 3709–3721.
37. F. Pelusi, S. Ascione, M. Sbragaglia and M. Bernaschi, Soft Matter, 2023, 19, 7192–7201.
38. L. Brigo, M. Natali, M. Pierno, F. Mammano, C. Sada, G. Fois, A. Pozzato, S. Dal Zilio, M. Tormen and G. Mistura, J. Phys.: Condens. Matter, 2008, 20, 354016.
39. J. R. Seth, L. Mohan, C. Locatelli-Champagne, M. Cloitre and R. T. Bonnecaze, Nat. Mater., 2011, 10, 838–843.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm00041b

---