Shear dynamics of confined bijels

Bijel represents a remarkable example of the composite soft porous material made of bicontinuous domains of two immiscible fluids (such as water and oil) frozen in a rigid state by a monolayer of colloidal particles irreversibly trapped at the fluid interface.

Such a material, first predicted by pioneering the lattice Boltzmann (LB) simulations¹ and then realized in the lab,^2–8 has been proposed for a variety of applications, ranging from energy storage and molecular encapsulation^4,6,9 to catalysis and tissue engineering.^2,3,8,10 A distinctive property of bijels is the arrested domain coarsening, a process in which the demix of a binary fluid is inhibited whenever colloidal particles, wetting both liquids, are adsorbed onto the fluid interface. Since the energy cost necessary to remove these particles is generally several orders of magnitude higher than the thermal energy,^1,11 the colloidal assembly jams at the interface and locks the system into an amorphous glassy-like configuration.¹² 

Manufacturing and design of a bijel remains a rather difficult task since its realization depends on careful control of several key physical parameters, such as the size and volume fraction of the particles, their affinity toward the two liquids, and liquid–liquid surface tension. Numerical simulations are essential to this scope since they provide direct access to experimental parameters, which can be selectively tuned in order to explore the kinetic pathway leading to the formation of the bijel. In addition, the physics of this material can be almost exclusively investigated by means of suitable computational models due to the complex structure of the governing equations involved.

Theoretical studies performed so far have investigated the routes to achieve a stable bijel in bulk samples^1,11,13,14 and only more recently in a confined region such as the thin fluid film;^15–17 whereas its rheological response under shear has been experimentally examined in Refs. 18 and 19. In this context, it has been shown that the breakup of the structure is achieved through an intermediate step in which the interface is fluidized and the bijel is allowed to move.¹⁹ 

In this paper, we present new large-scale simulations addressing the case where a bijel, confined within a three-dimensional rectangular channel, is subject to a symmetric shear flow sufficiently strong to cause its breakup. The physics is numerically studied by using LBsoft,²⁰ a recent multiscale lattice Boltzmann (LB) code built for large-scale simulations of colloidal fluids. Our results show that the shear is generally capable of removing a considerable number of colloids from the fluid interface toward the bulk fluid, thus yielding to further domain coarsening and loss of interfacial area. Fluid domains are found to preferentially align along the direction imposed by the shear flow and to attain a size comparable to the length of the channel. Such dynamics flattens the initial corrugated interface, and at high shear, it considerably diminishes the global curvature of the resulting colloidal fluid. Finally, particles favor the formation of a shear banding arrangement of the flow profile far from the walls. Those adsorbed at the interface maintain their orientation set by the wetting angle, while those accumulating in the bulk constantly spin and move around due to coupling with the surrounding fluid.

This paper is organized as follows. In Sec. II, we outline the LB approach used to simulate the physics of a bijel by also providing further details about computational performances. In Sec. III, we discuss numerical results, considering a bijel confined in a three-dimensional channel. Its response under shear is quantified in terms of the structure of the fluid flow, colloidal position, and orientation in the fluid and curvature dynamics of the fluid interface. Some final remarks conclude the paper.

The physics of the bijel is simulated using LBsoft,²⁰ an open-source software developed for solving the hydrodynamics of colloidal systems.^{13,14,21–23} LBsoft essentially couples a Shan–Chen lattice Boltzmann approach,^24–26 used to integrate the continuity and the Navier–Stokes equations for multiple fluid components with surface tension, with a discrete particle dynamics scheme to solve the Newton’s equations of momentum and angular momentum, tracking the dynamics of an arbitrarily shaped rigid body. Here, we shortly outline the method and refer to Ref. 20 for more details.

In the Shan–Chen model, the two components of a binary fluid are described by two discrete sets of distribution functions $f i k r , t$ (k = 1, 2) representing the probability of finding a “fluid particle” on a lattice site r at time t moving at velocity c_i along a predefined set of discrete speeds. LB soft employs a three-dimensional 19-speed cubic lattice (D3Q19) in which 19 discrete velocities c_i (i = 0, …, 18) link mesh points located at distance Δx (first neighbors) and $2 Δx$ (second neighbors). The local density ρ^k(r, t) of the kth component and the total momentum of the mixture ρu = ∑_kρ^ku^k are given by the zeroth and the first order moment of the distribution functions, i.e., $ρ k ( r , t ) = ∑ i f i k ( r , t )$ and $ρu= ∑ i ∑ k f i k ( r , t ) c i$⁠.

The dynamics of $f i k$ is governed by the following set of discrete Boltzmann equations:

in which $f i e q$ are local equilibrium distribution functions computed as a second order expansion in the fluid velocity u,

where $c s =1/ 3$ is the speed of sound, Δt is the lattice time step, $ω=2 c s 2 / ( 2 ν − c s 2 )$ is related to the viscosity ν of the mixture, I is the unit matrix, and ω_i is a normalized set of weights. In a D3Q19 lattice, one has ω₀ = 1/3, ω_1–6 = 1/18, and ω_7–18 = 1/36.²⁷ 

Finally, the term $S i k$ is a forcing contribution defined as

in which F^k is an extra forcing accounting for the interaction between the two fluid components. It is given by

where G_c is a coupling constant controlling the interaction strength between fluid particles. It is set to a positive value to promote phase separation. Finally, ψ^k is an effective density term, which coincides with the physical density of the fluid, ψ^k(r, t) = ρ^k(r, t).

Solid particles are modeled by following the approach initially introduced by Ladd^21,22 and further extended by Harting,^13,14 in which a closed curve mimicking the surface of a colloid is implemented by means of stick boundary conditions, bouncing back the density of the moving fluid and leaving the internal region of the colloid fluid-free.

Their wettability at a fluid interface is modeled by filling each solid node of a particle frontier with a virtual fluid density, taken as the average value at the fluid nodes in the nearest lattice positions (i_NP). The virtual density of the kth component is given by

where the factor ζ(r) tunes the local wettability of the solid node located at r. If ζ is positive, the particle surface prefers the kth fluid component, whereas if negative, it prefers the other one.

The equations of motion of a colloid of mass m_p are then given by

where u_p is the velocity of the colloid, ω_p is its angular momentum, and I_p is the moment of inertia. Finally, F_p and T_p are the force and the torque acting on the particle, respectively.

The former one is, in turn, the sum of five different contributions accounting for (i) the force F_l,p acting on the particle due to momentum transfer of the surrounding fluid, (ii) the forces F_d,p and F_c,p due to the change of a lattice node from fluid to solid (destruction of fluid) and from solid to fluid (creation of fluid), (iii) a repulsive force F_rep, which inhibits interpenetration of pairs of colloids as well as between colloids and walls, and, finally, (iv) a lubrication force F_lub necessary to model the fluid flow effects between two particles when their distance falls below the lattice resolution. These two terms are, respectively, given by

and

where K is an elastic constant, r_ij = r_j − r_i is the distance between i-th and j-th particles of radius R, h = |r_ij| − 2R, and h_n is the cutoff distance.

Similarly, the torque is the sum of three different terms, T_l,p, T_d,p, and T_c,p, accounting for the angular momentum transfer of the surrounding fluid and for the destruction and creation of fluid due to particle motion over the lattice. Further details are given in Ref. 20.

Simulations are run over a three-dimensional domain of size L_x = 128, L_y = 128, and L_z = 1024 under periodic boundary conditions along the mainstream axis z [see Fig. 2(a)]. Flat walls are set along the yz planes at x = 0 and x = L_x, and along the xz planes at y = 0 and y = L_y. Here, no-slip conditions hold for the velocity field, and neutral wetting is set for the phase field $ϕ= ρ y − ρ b ρ y + ρ b$ (ranging between −1 and 1), where ρ_y and ρ_b are the densities of the two components. In addition, we have set G_c = 0.65, K = 20, and h_n = 1. The wall-particle repulsion has the same functional form of Eq. (8) and is switched on when the distance between the center of mass of particle and the wall is <7.5 lattice units.

We initialize the system in a mixed state with a small uniformly distributed noise with −0.1 < ϕ < 0.1 and with walls at rest, and include N = 3600 solid spheres, of radius r = 5.5 lattice units, randomly dispersed in the channel. Their wetting angle is θ = 100° (i.e., the angle where the Shan–Chen force is maximum), which is assumed to be equal for all particles. Such a configuration corresponds to a volume fraction,

a value sufficiently high to form a stable bijel at the steady state.

Thereafter, the walls at y = 0 and y = L_y are moved along the opposite directions, the former with speed −u leftward and the latter with speed u rightward, while walls at x = 0 and x = L_x are kept at rest. This sets a shear rate $γ ̇ =2u/ L x$⁠. Once a quasi-stationary state is achieved, the shear force is turned off and the mixture is let to relax toward its new equilibrium state.

The dynamic scheduling essentially consists of three separate stages. In stage (1), after ∼1.5 × 10⁵ time steps, a bijel is formed, and this configuration is used as a starting point for the stage (2), in which a symmetric shear, lasting for t = 5 × 10⁵ time steps, is applied. Finally, the shear is switched off and the simulation is run for further t = 5 × 10⁵ times steps.

The behavior of the bijel is then inspected by varying the speed of the walls in the range 0.005 < u_w < 0.02, corresponding to a shear rate $γ ̇$ ranging from ∼4 × 10⁻⁵ to ∼1.5 × 10⁻³. If not otherwise stated, parameters of the model are set equal to the ones reported in Appendix D of Ref. 20.

All simulations have been performed on the UK Archer supercomputing facility composed of a cluster of 4920 computing nodes each with two 2.7 GHz, 12-core Intel E5-2697 v2 (Ivy Bridge) CPU, and 64 GB of RAM. The nodes are connected by using a Cray Aries network, with a DragonFly topology with a peak bisection bandwidth of over 19 013 GB/s over the whole system. We have used the Intel fortran compiler v18 with an optimization level of −O3.

In Table I, we report the wall-clock time measured as a function of the number of cores for two lattices of size 128 × 128 × 1024 (top, used in our simulations) and 256 × 256 × 2048 (bottom), to show the weak scaling behavior of the code. Simulations were run using 128 cores for 1 × 10⁶ time steps, resulting in 3.5 days of running time for each case.

The kinetic route as well as the mechanical properties to realize a well-formed and stable bijel is rather well established and is known to generally depend on surface tension of the mixture as well as on the size, volume fraction, and wetting angle of the colloids. The structural integrity of the material is ensured, for example, with colloids of size larger than the width of the interface, at a volume fraction higher than 15%, and with a wetting angle around 90°.^11,28–31

Numerical simulations performed so far^1,13,14 have studied bijels within cubic lattices under periodic boundary conditions, but only, more recently, in a confined region, such as a thin fluid film.¹⁵ In Fig. 1 (left), we show a typical example of such structured fluid in a three-dimensional cubic box of linear lattice size L = 128, obtained by dispersing N = 450 particles of radius r = 5.5 lattice sites (which sets a volume fraction of ∼15%) adsorbed onto the interface and wetting both components of a bicontinuous fluid. The characteristic bulged soft porous matrix results from the arrest of the domain coarsening due to irreversible particle adsorption and occurs after approximately t = 2 × 10⁵ time steps, in good agreement with previous works.^1,14,20

Recent experiments have investigated the response of this material to an external solicitation, such as an oscillatory shear stress,¹⁹ and have found that the bijel essentially undergoes a two-step yielding process: in the first one, interfaces are fluidized by the shear flow and destabilize the structure, while in the second one, the material is permanently broken.

In Sec. III A, we focus precisely on the dynamics observed after the breakup and provide a selected number of statistical indicators to quantitatively assess the mechanical properties.

In Fig. 2, we show, for example, the evolution of a confined bijel under a symmetric shear (⁠$γ ̇ ≃3×1 0 − 4$⁠, v_w = 0.02) in a dynamic cycle. The shear is turned on once the bijel is close to its equilibrium (a) (a weak domain coarsening still occurs, see Fig. 6), and it is turned off when the material resulting after rupture attains a quasi-stationary state (c). Note, in particular, that the morphology of the equilibrated bijel in the channel (a) is characterized by bicontinuous fluid domains mixed with discrete ones stabilized by spherical colloids, and closely resembles those reported in Ref. 15, in which a detailed study of shapes obtained in the thin fluid film has been reported.

Due to the high shear rate, two concurrent effects take place. While highly anisotropic fluid domains, of size approximately half of the channel length, form and align along the flow, a fraction of colloids, initially sequestered at the fluid interface, detaches and migrates away dragged by the fluid flow (b). This promotes contacts among interfaces and allows for domain coalescence, favored by the detachment of particles accumulating in the bulk of the fluid. At the steady state, a single highly stretched fluid domain lies along the channel, while colloids ceaselessly spin and travel around transported by the fluid (c). Once the shear is turned off, the relaxation only mildly affects the shape of the remaining fluid domain, while particles progressively slow down, until the whole mixture equilibrate (d). The tenuous differences between sheared (c) and un-sheared (d) state suggest that the mixture has been driven toward a long-lived configuration stable over weak perturbations, likely, a further minimum of the complex free energy landscape.

In the following, we quantitatively characterize the dynamic properties of such complex processes in terms of (i) fluid structure and velocity profile within the channel, (ii) particle position and orientation, and (iii) curvature dynamics of the fluid interface.

In Fig. 3, we show a three-dimensional cross section of the velocity field within the microchannel for $γ ̇ ≃3×1 0 − 4$⁠. The bidirectional flow is higher near the walls and gradually weaker toward the center of the device. Green rounded regions identify the instantaneous position of the colloids (whose interior is fluid-free), located both in the bulk and close to the fluid interface. A typical velocity profile, computed at different positions along the channel and mediated along the y-direction, is provided in Fig. 4. It is “smoothed” over the particle location, or, in other words, it is reconstructed by using the truncated Fourier method, which cuts high frequencies of the Fourier transform of the phase field emerging in the surrounding of each colloid. Such a technique proves useful to unveil a clear shear banding-like signature in the middle section of the channel, where the speed decreases, while it correctly captures the rapid increase in the velocity near the walls. This effect is very much likely due to the colloids, whose presence prevents the formation of a fully developed linear profile, such as the one of a single fluid under Couette flow. Indeed, the shear banding effect persists as long as particles obstruct the flow (far from the walls) and exhibits local deformations facing the peaks of P_N(x), capturing the position of the colloids.

This picture suggests that colloids play a crucial role in controlling the mechanical response of a bijel (and of the resulting sheared colloidal fluid), subject to an externally applied forcing. Hence, there is a significant scope for monitoring their dynamics, in particular, position and orientation at the fluid interface and in the bulk of the mixture.

In Fig. 5, we show a typical arrangement of colloids under shear, both at the interface (in green) and in one of the fluid components (in dark).

As discussed in Sec. III A, an intense shear flow favors the detachment of a relevant fraction of particles initially stuck at the fluid interface and their migration toward the bulk of the fluid. This dynamics can be quantitatively monitored by the cumulative distribution of colloids C_N in the mixture, defined as C_N(ϕ) = ∫N(ϕ, r)dr (see Fig. 6). If $γ ̇ =0$⁠, such a process is very mild. Indeed, the initial sigmoidal profile undergoes mild modifications over time, since a weak domain coarsening still occurs and a few particles move toward the bulk (where ϕ ≃ ±1). In the presence of shear, the detachment is enhanced as $γ ̇$ increases, and the profile of C_N significantly sharpens, becoming essentially “digitized” for high values of the shear rate. Once the shear is turned off, particle migration proceeds, which is driven only by domain growth, except where the high shear rate has led the mixture to a stable state, such as the one in Fig. 2(d).

Such complex dynamics also significantly affects the orientation of the particles in the mixture. In Fig. 5, we show, for instance, a typical configuration of colloids both at the fluid interface and in the bulk, each one with a local frame of reference. While, for the particles segregated at the fluid interface, it is fixed (i.e., the local x-axis is essentially parallel to the normal n at the interface), for those in the bulk it looks randomly oriented, since it generally depends on the complex coupling between the local transfer of momentum and torque of the surrounding fluid to the particle. This behavior is nicely captured by the “eyeball” scatterplots shown in Fig. 7, where the cosines of the angles formed between each axis r_i (i = 1, 2, 3) of the reference frame of a particle and the normal to the fluid interface n are reported. The angles are defined as $cos ( α i ) = r i ⋅ n | r i ‖ n |$⁠, where r₁ = (1, 0, 0), r₂ = (0, 1, 0), and r₃ = (0, 0, 1). If the alignment is perfect, α = 0, hence, cos(α) = 1.

The red spot in the lower right part of the eyeball indicates that the x-axis of the particles located at the fluid interface lies approximately parallel to n, since cos(α₁) ≃ 1, cos(α₂) ≃ 0, and cos(α₃) ≃ 0 (see also Fig. 5), while those located in the bulk distribute more uniformly on the sphere. The progressive shrinking of the spot for increasing values of the shear rate suggests, once again, that a significant fraction of colloids move from the interface toward the bulk, although the orientation of the former ones, set by the wetting angle θ = 100°, remains essentially unvaried under shear.

The results discussed so far outline a picture in which the complex interplay between solid particles and fluid interfaces dramatically influences the topological properties of the material. The latter ones are also critically controlled by the fluid interface, which exhibits a rather complex dynamics under shear. Further insight about its behavior can be gained by tracking the evolution of its interface curvature over time and for different values of shear rates. Section III D is dedicated to elucidate this point.

Previous studies on colloidal mixtures³² have shown that the probability distribution function (pdf) of the local curvature p(k) provides a further indicator capturing features such as fluid domain coarsening, inhomogeneities in the local particle volume fraction, and the underlying dynamics of the fluid interface.

In Fig. 8, we investigate the dynamics of the magnitude of the unsigned curvature, defined as $k=|∇⋅ ∇ ϕ ‖ ∇ ϕ ‖ |$⁠. In order to provide an approximate global evaluation of k in the system, we reconstruct the interface in regions occupied by solid particles through a truncated Fourier method. This allows us to remove inevitable artificial effects due to the corrugated structure of the particles, which would overestimate the contribution of higher curvatures and to smooth the profile, thus providing an easier assessment of its trend.

In the absence of shear, the curvature undergoes negligible modifications over time, indicating that, although a weak domain coarsening still occurs, the bijel is overall stable within the temporal window investigated. Increasing the shear rates favors a transition from the initial broad-shaped structure toward a more localized distribution at lower values of curvature, since the shear stretches and flattens the interface along the flow direction (z in our simulations). A perfectly flat interface would result in a Dirac-like distribution centered in k ∼ 0. Once the shear is turned off, the pdf only mildly shifts toward lower values of k, except for the case shown in Fig. 8(d) in which a long-lived stationary state has been already attained under shear and no further relevant shape modifications occur.

It is finally worth noting that the late time unsheared pdfs can be best fitted by a stretched exponential function whose exponent is generally higher than 1, regardless of the presence of a shear force. A similar trend has been also reported in unbound bicontinuous colloidal mixtures in which dumbbells (modeled by the immersed boundary method) are confined at the interface. Here, however, unlike the bijel, the exponent of the steady-state interfacial curvature was found sensibly lower than 1. Although we currently lack a fully convincing explanation, a potentially relevant reason giving rise to such a difference may lie in the fact that the colloidal–dumbbell mixture was not found to exhibit any arrested phase separation, thus, probably allowing for an alternative dynamic evolution of the fluid interface. In addition, the presence of confining walls, as well as different wetting angles and shapes of colloids, may also partially influence the relaxation of the interface.

In summary, we have numerically investigated, by means of large-scale lattice Boltzmann simulations, the dynamic response of a bijel confined within a three-dimensional rectangular channel, subject to a symmetric shear sufficiently intense to induce the rupture of the material.

Our results show that the shear flow promotes the “melting” of the material by suppressing the arrested coarsening (the mechanism leading to the formation of the bijel), essentially because a sizable fraction of the colloidal particles fails to stick to the interface. While particles detach from the fluid interface, coarsening proceeds and domains stretch and elongate along the shear direction, attaining a steady-state of size comparable with the longest dimension of the channel. Such an effect is increasingly prominent at higher rates.

The complex dynamic behavior of the material is studied in terms of a series of quantitative indicators, such as the fluid structure in the channel, the spatial distribution of colloids, and the curvature of the fluid interface. The fluid flow exhibits a characteristic shear-banding-like signature, basically due to colloids, which hinder the formation of a linear velocity profile. Particles also display distinctive features. While those sequestered at the interface, perhaps counterintuitively, maintain their orientation (essentially set by the wetting angle) even under shear, the ones in the bulk, whose number increases at higher shear rates, spin around due to the momentum and torque transfer from the surrounding fluid. These processes produce significant modifications of the fluid interface, which turns from a highly corrugated entangled arrangement with large curvature to a smooth and rounded one with low curvature.

Our work sheds light on the mechanical response of a confined bijel subject to shear flow capable of jeopardizing the structural integrity of the material. A route for preventing such a dramatic event, already partially explored in stabilized monojels,³³ may consider the use of attractive colloids (or a combination of attractive and repulsive ones), which could potentially suppress the melting of the material by providing further resistance to shear deformations.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors acknowledge funding from the European Research Council under the European Union’s Horizon 2020 Framework Programme (Grant No. FP/2014–2020), ERC Grant Agreement No. 739964 (COPMAT), and the PRACE 16DECI0017 RADOBI project. Supercomputing time on ARCHER is provided by the “UK Consortium on Mesoscale Engineering Science (UKCOMES)” under the UK Engineering and Physical Sciences Research Council, Grant No. EP/R029598/1.