Mesoscale modeling of near-contact interactions for complex flowing interfaces

We present a mesoscale kinetic model for multicomponent flows, augmented with a short range forcing term, aimed at describing the combined effect of surface tension and near-contact interactions operating at the fluid interface level. Such mesoscale approach is shown to i) accurately capture the complex dynamics of bouncing colliding droplets for different values of the main governing parameters, ii) predict quantitatively the effective viscosity of dense emulsions in micro-channels and iii) simulate the formation of the so-called soft flowing crystals in microfluidic focusers.


Introduction
A thorough knowledge of the dynamic interactions between fluid interfaces is paramount to a deeper understanding of a variety of natural processes and engineering applications, such as combustion, microfluidic coating, food processing and many others. The coalescence and/or repulsion between droplets or bubbles can be traced to the hydrodynamic drag originating from the relative motion of two fluid interfaces in near contact ( [1,2,3,4,5]), and to the combined action of nanoscale attractive and/or repulsive forces, such as Van der Waals and electrostatic forces, steric interactions, hydration repulsion and depletion attraction ( [6,7]).
A wide body of theoretical and experimental work has elucidated the complex nature of the near contact interactions which develop within intervening liquid films: from the pioneering works of Gibbs and Marangoni on the thermodynamics of liquid thin films (see ( [7]) for a comprehensive review) to the separate work of Derjaguin and Overbeek ( [8,9]) which culminated in the joint DLVO theory (after its founders, Derjaguin, Landau, Verwey and Overbeek). These seminal works laid down the foundations for describing a broad variety of complex flowing systems, such as colloids, foams and emulsions, as well as flowing collections of droplets and bubbles characterized by highly ordered and uniform, crystal-like structures, now known as soft flowing crystals ( [10,11]). From a numerical standpoint, the direct introduction of near interaction forces at a molecular level reflects the need of solving simultaneously six spatial decades: from millimeters, i.e. the typical size of microfluidic devices, down to nanometers (and below), namely the relevant spatial scale of contact forces. This is far beyond the capability of any foreseeable computer ( [12]), hence placing a high premium on coarse-grained, mesoscale representations of near-contact forces, capable of retaining computational viability without compromising the essential physics. Of course, the success of such mesoscale strategy hinges crucially on the universality of the underlying physics, i.e. its dependence on suitable dimensionless parameters measuring the relative strength of the interactions, rather than on the details and strength of the interactions themselves ( [13]). Failing such universality, a genuine microscopic approach cannot be helped, thus hampering the possibility to reach up to the scales of the full device.
In the LB framework, immiscible fluids were first modelled by [14].These Authors developed a lattice Boltzmann model, augmented with a forcing term in accordance with the local color gradient, giving rise to the interface tension and a segregation step. In [15], the Authors corrected the segregation rule proposed by Gunstensen, to avoid the pinning problem that affected the Gunstensen model, allowing the fluids to moderately mix and to keep the color distribution symmetric with respect to the color gradient. The reason for lattice pinning is that at the sites where it happens, all of the particles of one kind are sent to the same direction and hence they cannot move from one site to another. More recently the model was further improved to simulate high density and viscosity ratio ( [16,17,18]). In addition to this, there have been several attempts to model interface interactions in amphiphilic fluids within the LB framework (see for example [19,20,21]). These models aim at directly simulating the effect of a surfactant by evolving two sets of distribution functions, allowing to take into account the presence of an amphiphilic fluid at the interface between two liquid phases. For a comprehensive review of multiphase and multicomponent LB models ( [22,23,24,25,26]) , the interested reader is referred to [27].
In this paper, we present a lattice Boltzmann (LB) based approach for multicomponent flows, based on the color-gradient model [16], augmented with an additional forcing term which is aimed at representing the effects of the near-contact forces operating at the fluid interface level. The proposed model, based is shown to accurately capture the collision outcomes between bouncing droplets for different values of the governing parameters, to predict the effective viscosity of dense emulsions in channels and to effectively simulate the evolution of soft flowing crystals in flow focuser devices. The paper is organized as follows. In Section 2 the lattice Boltzmann equation with the BGK collisional operator is described, together with the color gradient model and the regularization algorithm for simulating multicomponent fluids. The augmented LB model for repulsive near contact interactions is discussed in details in subsection 2.2. Section 3 collects the main results of the paper. Finally, a summary is reported in section 4.

Method
Lattice Boltzmann models for non-ideal fluids come mainly in two families. The first one is based on heuristic assumptions ( [16,17,28,29,30]) while the second one, builds on the projection of the kinetic equation on a discrete set of microscopic velocities ( [22,23,24,25,31,26]). For a more exhaustive review, see ( [27,32]). In the following we provide a brief introduction to the one which we found most suitable for the description of flowing crystals, namely the regularized color gradient method ( [29]).

Regularized color gradient lattice Boltzmann model
In the color gradient LB for multicomponent flows, two sets of distribution functions are needed to track the evolution of the two fluid components, which occurs via a streaming-collision algorithm (for a comprehensive review on the lattice Boltzmann method, please refer to ( [32,33])): where f k i is the discrete distribution function, representing the probability of finding a particle of the k th component at position x and time t with discrete velocity c i .
where and i is the index running over the lattice discrete directions i = 0, ..., b, where b = 26 for a three dimensional 27 speed lattice (D3Q27). The lattice time step ∆t has been taken as 1 (in lattice units) for convenience, which is a common practice in LB literature (see [32]). The density ρ k of the k th component is given by the zeroth moment of the distribution functions: The total fluid density is given by ρ = k ρ k , while the total momentum of the mixture is defined as the sum of the linear momentum of the two components: The collision operator can be split into three parts ( [14,16,17]): In the above, Ω k i (1) stands for the standard collisional relaxation ( [34]) which reads Ω k where ω ef f = 2/(6ν − 1) is the effective relaxation parameter beingν the viscosity at the interface between the two fluids which is computed as 1 ν = ρ1 ν2 (ν 1 and ν 2 are the kinematic viscosities of the two fluids in the bulk). The equilibrium distribution function of the k th component f k,eq i is given by a low-Mach, second-order, expansion of a local Maxwellian, namely f k,eq (2) is the perturbation step ( [14]), which contributes to the build up of an interfacial tension. Finally, is the recoloring step ( [14,15]), which promotes the segregation between species, so as to minimize their mutual diffusion.
In order to reproduce the correct form of the stress tensor ( [35]), the perturbation operator can be constructed by exploiting the concept of the continuum surface force ( [36]). Firstly, the perturbation operator must satisfy the following conservation constraints: By performing a Chapman-Enskog expansion, it can be shown that the hydrodynamic limit of Eq.1 is represented by a set of equations for the conservation of mass and linear momentum: where p = k p k is the pressure and ν = c 2 s (τ − 1/2) is the kinematic viscosity of the mixture, being τ the single relaxation time and c s = 1/ √ 3 the sound speed of the model ( [34,33]). Note that the divergence of the stress tensor (last term in equation eq. 2.8), which is responsible for the build up of surface tension, acts only at the interface between the fluids (see eq. 11) The stress tensor in the momentum equation is given by: The surface stress boundary condition at the interface between two fluids can be expressed as follows ( [35,36]): where, I is the identity tensor, σ is the surface tension coefficient (with dimensions of force per unit area), n is the unit normal to the interface, T = −pI + ρν(∇ u + ∇ u T ) is the stress tensor of the k th component and ∇ · n is the local curvature of the fluid interface. The local stress jump at the interface can be induced by adding an interfacial volume force F (x, t) ( [37]): In the above, δ I = 1 2 |∇Θ| is an index function which explicitly localizes the force on the interface and Θ = ρ 1 −ρ 2 ρ 1 +ρ 2 is the phase field ( [37]). The normal to the interface can be approximated by the gradient of the phase field, n = ∇Θ/|∇Θ|.
Since the perturbation operator is responsible for generating interfacial tension, the following relation must hold: By choosing ( [16]) Ω k |∇Θ| 2 − B i , substituting it into 5 and 12 and by imposing that the set B i must satisfy the following isotropy constraints: we obtain an equation for the surface tension of the model: The above relation shows a direct link between the surface tension and the parameter A = A 1 + A 2 (A 1 = A 2 ). In actual practice, after choosing the viscosity of the two components and the surface tension of the model, at each time step, one locally computes the A = A 1 + A 2 coefficient by using the formula reported in eq. (14). It is worth noting that, in this work, a fourth order isotropic discrete gradient operator on a 27 points stencil is employed (for details please refer to [17]).
As pointed out above, the perturbation operator generates an interfacial tension in compliance with the capillary-stress tensor of the Navier-Stokes equations for a multicomponent fluid system.
Nonetheless, the perturbation operator alone does not guarantee the immiscibility of different fluid components. For this reason, a further step is needed (i.e. the recoloring step) to minimize the mutual diffusion between components.
Following the work of Latva-Kokko and Rothman ( [15]), the recoloring operator for the two sets of distributions takes the following form: denotes the set of post-perturbation distributions, ρ = ρ 1 + ρ 2 , cos φ i is the angle between the phase field gradient and the i th lattice vector and f eq,0 Note that the coefficient β in the above expressions is a free parameter which can be used to tune the interface width, thus playing the role of an inverse diffusion length scale ( [15]). We wish to point out that the present model is employed to simulate droplet-based microfluidic applications which are often characterized by very small W e, Re and Ca numbers (i.e. much smaller than one). Hence, by considering typical flow speeds of 10 −3 − 10 −2 lu/step the u 3 error contribution is of the order O(10 −12 − 10 −9 ), which reflects in a very negligible compressibility errors. It is important to note that, following the work of [16],in this work we perform the entire collision step (collision+perturbation+recoloring steps) on two separate distibutions, this at variance with the works of [14] and [15], in which the collision and perturbation are written in terms of the blind distribution The color gradient LB scheme is further regularized by filtering out the high-order non-hydrodynamic (ghost) modes, emerging after the streaming step (see ([38, 39, 40, 41, 42, 43, 44]) for further details).
Indeed, it was noted that sizeable non-isotropic effects arise in the model ( [29]), whenever the LB scheme is under-relaxed (τ ≥ 1). As a consequence, we exploit the regularization procedure in order to recover the loss of isotropy by suppressing the non-hydrodynamic modes ( [45,29,46]). For the sake of clarity, here we report a pseudocode of the regularization procedure employed in our simulations: In the next subsection, we show how to include the effect of repulsive near contact interactions, directly within the LB framework, by augmenting the regularized color-gradient model with a forcing term aimed at coarsegraining the near-contact forces at the fluid surface.

Augmented LB model for repulsive near contact interactions
The stress-jump condition across a fluid interface is augmented with a repulsive term aimed at providing a mesoscale representation of all the repulsive near-contact forces (i.e., Van der Waals, electrostatic, steric and hydration repulsion) acting on much smaller scales (∼ O(1 nm)) than those resolved on the lattice (typically well above hundreds of nanometers). It takes the following form: where π[h( x)] is responsible for the repulsion between neighboring fluid interfaces, h( x) being the distance along the normal n, between locations x and y = x + h n at the two interfaces, respectively (see Fig. 1).
The above expression can be rewritten in the following form ( [47]): in which ∇ s is used to identify the gradient tangent to the interface. By neglecting any variation of the surface tension along the interface, we can approximate T = −pI ([36]) and the above equation takes the following form: By projecting the equation along the normal to the surface, we obtain the extended Young-Laplace equation ( [48,49]): The additional term can be readily included within the LB framework, by adding a forcing term acting only on the fluid interfaces in near contact, namely: In the above, A h [h( x)] the parameter controlling the strength (force per unit volume) of the near contact interactions (please refer to the sketch in figure 1). In this work, A h is set to a constant value (A hM ) if h < h min and then decreases as ∼ h −3 , as shown in figure  1, although other choices are certainly possible. The near contact force has been defined solely as a function of the distance between two fluid interfaces. Nonetheless, it could be easily extended to take into account local variations due to the effect of spontaneous migrations of the surfactant along the fluid interface ( [50]). The addition of the repulsive force, naturally leads to the following (extended) conservation law for the momentum equation: This is the Navier-Stokes equation for a multicomponent system augmented with a surface-localized repulsive term, expressed through the potential function ∇π.
There have been other attempts to model interface interactions in amphiphilic fluids in the literature (see for example [19,20,21]). These models aim at directly simulating the effect of a surfactant by evolving two sets of distribution fucntions, allowing to take into account the presence of an amphiphilic fluid at the interface between two liquid phases. Indeed, the propagation of the amphiphilic molecules is described a set of LB equations, one for the distribution function and one for the relaxation of the average dipole vector to its local equilibrium orientation. [51] and [52] extended the color gradient model to any number of components. This method makes use of a set of distributions for each immiscible fluid species. More recently, [53] proposed a thermodynamically consistent free energy model for fluid flows comprised of one gas and two liquid components using the entropic lattice Boltzmann scheme. The standpoint of our work is quite different in that, the repulsive action of a surfactant, arising when two interfaces come in close contact, is taken into account by introducing a repulsive forcing term localized at the interface of the two fluids. Importantly, this just requires two distributions functions regardless of the number of droplets.
To conclude, the extended approach still holds to a continuum description of the interface dynamics, being the governing equations modified only by the presence of a distributed body force, which can heuristically be interpreted as a coarse-grained version of the short-range molecular forces acting at the nanometers and sub-nanometer scales.

Numerical Results
In this section we test the extended LB model on three applications namely, head-on and off-axis collision between two bouncing droplets, pressure driven flow of a dense emulsion in a channel and the formation of soft flowing crystals in a microfluidic flow focuser.

Droplets coalescence
Here we first test the ability of the color-gradient LB model to capture the physics of the coalescence process between two equally sized droplets. The simulation setup consists of a three-dimensional fully periodic box (160 × 100 × 121) in which two liquid droplets of radius R = 30 , surrounded by a dispersed fluid of the same density, are placed at close distance. The surface tension of the mixture was fixed at a constant value, σ = 0.01 while the surface tension of the two fluids was varied from ν = 0.05 to ν = 0.15. The viscosity ratio between the droplets and the surrounding phase has been set to unity, and all the physical quantities are reported in lattice units. As shown in fig.2 (panel a-c), a liquid bridge forms between the two droplets. As pointed out in ( [54]) the coalescence process in the early stage is so fast that only a fraction of the fluid caught in the narrow gap between the two spheres is able to escape, while the rest accumulates in a "bubble" which forms at the meniscus. The simulations predict the formation of such a bubble, which remains trapped between the two coalescing droplets. To be more quantitative, we measured the normalized radius of the liquid bridge, r b /R, as a function of the square root of the non-dimensional time t/τ (panel d), being r b the bridge radius, t the simulation time and τ a characteristic time scale defined as τ = R 3 /σ. As evidenced in the figure, the growth of the liquid bridge follows a scaling law r b ∝ (t/τ ) 0.5 , with the data collapsing on a single master curve r b /R = 1.6 t/τ . It is also worth noting that, the prefactor predicted by the simulations ( 1.6), is in very close agreement with the value 1.62 reported in ( [55]).

Bouncing colliding droplets
In this subsection, we show the capability of the extended LB model to accurately reproduce the correct dynamics of head-on and off-axis collisions between two bouncing droplets ( [56]).
The characteristic nondimensional parameters, governing the collision outcome, are the Weber and the Reynolds numbers, defined as W e = ρU 2 rel D/σ and Re = U rel D/ν, respectively. In the above, U rel the relative impact velocity, D the droplet diameter, σ the surface tension coefficient and ν the kinematic viscosity, as well as the impact number b = χ/D, namely, the distance between the collision trajectories in units of the droplet diameter. In Table 1, we report the main simulation parameters (expressed in lattice units). Also in this case, the viscosity ratio between the droplets and the surrounding phase has been set to unity. Figure  3 shows three collision sequences for different impact, Weber and Reynolds numbers. The collision outcomes are compared with those reported in ( [56]). The experiments were performed with near millimetric droplets of immiscible fluids, with diameters ranging between 700 − 800µm and impact velocities in the range of 1 − 3m/s. The other relevant parameters can be found in ( [56]). By taking a droplet diameter of 700µm, discretized with 30 lattice units, we obtain a lattice spacing ∆x ≈ 20µm = 1lu, which is the minimum interaction distance between the simulated droplets. It is worth noting that, the thin film between two interfaces in near contact has characteristic lengthscales of the order of nanometers, far below the spatial scales accessible to our simulations. Indeed, in our simulations, the smallest spatial scale is around 20 µm (∆x), while the characteristic distance between two non coalescing impacting droplets is of the order of one to ten nanometers, i.e. three orders of magnitude smaller than the grid size. Nothwithstanding this gap, by matching the main governing parameters (Weber and Reynolds), the Figure 4: Sequence of the flow field during the impact between the two droplets (mid-plane slice). As shown in panel (b), during the first stage of the collision the dispersed phase flows outwards, allowing the two droplets to approach. Afterwards, when the droplets come in close contact, the fluid within the thin film begins to recirculate, thereby preventing film rupture hence the coalescence between the droplets. A remark is in order: in real interacting systems the thin film between two interfaces in near contact develops on characteristic lenghscales of the order of the nanometers, far below the spatial scales accessible to our simulations. Indeed, in panel (b), the smallest spatial scale is around 20µm while the characteristic distance between two non coalescing impacting droplets is of the order of the tenth of nanometers, i.e. three orders of magnitude smaller than the grid size.  overall impact dynamics is correctly captured by the simulations. This, again, calls into cause the universality of the underlying physical processes, i.e., at the spatial scale at hand, the interaction physics depends upon the dimensionless numbers, measuring the relative strength of the interactions, rather than on the strength of the interactions themselves. Panel (a) reports the sequence of an head-on collision between two equally sized droplets.
As expected, at these Weber and Reynolds numbers, the two droplets bounce off without coalescing, because the coalescence is frustrated by the effect of the near-contact repulsive forces. The bouncing collision also occurs by increasing the impact parameter (panel (b) and (c)). Indeed, in both cases, the impact velocity is not sufficient to break the intervening thin film and a kissing-like collision is finally observed. We then inspected the evolution of the thin liquid film during the collision process. As reported in fig. 4, in a first stage, the fluid between the two droplets flows outwards, thus allowing the droplets to approach each other (panel b, left). Afterwards, as they get closer, the liquid begins to recirculate inwards within the intervening film, stabilizing it and preventing the coalescence between the colliding droplets (panel b, right). This phenomenon resembles the so-called Marangoni flow in liquid films, namely a fluid recirculation occurring in liquid thin films in the presence of shear and temperature gradients, which was observed to prevent the coalescence between bodies of liquids ( [57]) . It is straightforward to note that, by varying the magnitude of the repulsive force, it is possible to promote or inhibit the coalescence of the impacting droplets. From this standpoint, it proves expedient to introduce two further non-dimensional parameters, which measure the relative strength of inertia and surface tension versus repulsive forces. The first, (S σ ) is defined as the ratio between the repulsive force (∼ A hM ) which frustrates the coalescence between droplets, and the surface tension (σ), which promotes it. The second (S κ ), is defined as the ratio between the local impact kinetic energy and the work done by the repulsive forces to prevent the two interfaces from coalescing.
The two non-dimensional parameters read as follows:  Table  1, except for the repulsive parameter, which is four times smaller. It is evident that the repulsive force is not strong enough to frustrate the coalescence between the two impacting droplets, as instead occurs in the case reported in the lower panel (S σ = 1, S κ = 0.36). In the above ∆x, the lattice spacing, is the characteristic length scale of near-contact interactions between two droplets on the lattice. In our simulations S σ = 1 and S κ = 0.36, meaning that the inertial and the surface forces (both promoting coalescence) are not strong enough to withstand the effect of the local repulsion. By lowering the intensity of the repulsive forces of a factor four ( S κ = 1.5 and S σ = 0.25), we finally observe the coalescence between the impacting droplets (see fig. 5). It will be shown below (subsection 3.4) that the balance between inertia, surface forces and repulsive actions crucially affects the formation and the overall dynamics of crystal-like structures in foams and emulsions.

Dense emulsion in a planar channel
We now discuss the pressure-driven flow of a dense emulsion within a narrow channel, made of a regular arrangement of equal sized droplets As φ increases, the velocity profiles flatten in the central region of the channel. In order to quantify the effect of the packing fraction, we measured the effective (or relative) dynamic viscosity, defined as the flow rate ratio µ ef f = Q(φ = 0)/Q(φ) and compared it against the model proposed by [59] and [58]. Taylor theory, the effective dynamic viscosity (for ν d /ν c = 1) can be expressed as a linear function of the volume fraction (in the limit of small droplets and volume fractions) as follows: The effective dynamic viscosity predicted by Bullard, which is based on the differential effective medium theory, read as follows: being [η] the intrinsic viscosity, whose value is restricted between the undeformable and the freely deformable limit ( [60]). As reported in figure 7 (left panel), the simulation results are in a very close match with respect to the theoretical prediction of Taylor (inset in left panel) and Bullard (with intrinsic viscosity set to [η] = 1).

Soft flowing crystals in a microfluidic focuser
As an application, we simulated the formation of oil/water emulsions in a microfluidics flow focusing device, whose sketch is reported in fig. 8. The micro-device is made of three channels supplying the dispersed (A) and the continuous (B) phase, plus an orifice (C) placed downstream of the three coaxial inlet streams. The mechanism of droplet formation follows from the periodic pinch-off of the dispersed jet by the continuous stream and the pinch-off mechanism takes place in the small orifice.
Flow focusers are nowadays widely employed for the production of mono-dispersed emulsions, due to the precise control over the emulsion monodispersity and the droplet size ( [61,62,63]). The high degree of flow reproducibility is due to the dominance of the viscous forces over inertia which smoothes flows and tames hydrodynamic instabilities. The pinch-off process is thus metronomic, allowing to produce droplets with measured standard deviations in size as little as the 0.1% ( [10,64,65]). Here, we show that the mesoscale approach proposed in this paper is able to reproduce different packing configurations in the outlet channel of the flow focuser, which are obtained by varying the dispersed-to-continuous flow rate ratio. Moreover, we also show that, by tuning the magnitude of the near force interaction, different structures of the flowing crystals can be achieved. The simulation setup is sketched in figure 8, while the main simulation parameters are reported in Table 2. The prototypical focuser, used in the simulations, consists of three main channels (A and B) of width H = 200µm, a nozzle (C) h = 100µm and an outlet channel of width H c = 400µm, while the height of the focuser is W = 100µm. By taking an interfacial tension of an oil-water mixture (∼ 50mN/m), the dynamic viscosity of the water (dispersed phase) µ ∼ 10 −3 P a · s and an inlet velocity of the dispersed phase ∼ 0.1m/s we obtain   Figure 9 reports two different flow configurations which can be obtained by varying the flow rate ratios between the dispersed and the continuous phase. In both cases, droplets readily self-assemble in ordered patterns described as soft flowing crystals ([11, 10, 66]. In panel (a) (φ = 1/2 ), the reported sequence shows a typical wet foam-like configuration in which the droplets are circular (cylinders in three dimensions) and automatically arrange on three rows, as evidenced in fig.9 panel (c), which also provides a visual comparison with the experimental data reported in ( [67]). Panel (b) ( φ = 1/1) shows a more ordered crystal structure, made of larger cylindrical droplets disposed along two staggered rows. We further investigated the effect of the magnitude of the near force on the formation of the crystal pattern. Figure 10 panel (a) shows a time sequence of the flow field during the breakup stage occurring downstream the striction of the focuser, with a magnitude of the near contact force lowered by a factor 5 with respect to the base case of fig.9 panel (b). It is evident that the magnitude of the repulsive force is no longer sufficient to prevent the coalescence between the droplet downstream the striction and the upstream jet, due to the high impact velocity. Thus, as evidenced in panel (c), the droplets in the outlet channel are about twice the size with respect to the previous case, and proceed in single-file motion, i.e. aligned along the horizontal axes of the focuser. It is also interesting to note that the action of the repulsive force keeps frustrating the coalescence in the outlet channel, where the typical velocity are much lower than at the exit of the nozzle. Once again, we can compute the non-dimensional parameter S κ , downstream of the nozzle, where the coalescence occurs. By taking a characteristic jet velocity U J ∼ 0.1, we estimate S κ ∼ 1.4 ( fig. 10 panel (a)) and S κ ∼ 0.18 ( fig. 10  panel (b)). Thus, S κ ∼ O(1) may be interpreted as a threshold above which the repulsive force is no longer able to balance the inertia, so as to frustrate the coalescence between the jet and the droplet. These preliminary simulations clearly highlight the pivotal role of the near contact interactions on the structure of the resulting flowing crystal. From this standpoint, the proposed approach may open up new chances to investigate the complex dynamics of flowing microfluidics crystals, helping in identifying the optimal operational regimes required to precisely control the production of mono-dispersed emulsions.  In the first case the strenght of the repulsive force is not sufficient to prevent the coalescence between the droplet and the jet, due to the high speed of the jet at the outlet of the nozzle. Nonetheless, the repulsive force is strong enough to prevent the coalescence between droplets moving in the outlet channel (see panel (c), left figure). In the second case, the repulsive force is strong enough to completely suppress the coalescence between the jet and the newly formed droplets, thus allowing the formation of an ordered soft flowing crystal in the outlet channel. experimental researchers endeavored to explain the basic physics behind near contact interactions. Notwithstanding the surge of experimental and theoretical activity, the numerical description of the soft flowing matter is still in its early stage. One reason is that, the concurrent solution of the macroscopic equations, needed to evolve the fluid interface, together with a direct description of the near contact interactions at the nano-scales stands out as a prohibitive multiscale problem.
In this paper, we have proposed a coarse-grained approach to include the effect of the near contact interactions within the lattice Boltzmann computational framework, and shown that such extended LB model is able to accurately describe a number of relevant physical effects. That is, i) to capture the evolution of two bouncing impacting droplets for different values of the main governing parameters namely, Weber Reynolds and impact number; ii) to predict the effective viscosity of a dense emulsion flowing in a micro-channel, in agreement with the theoretical model of Bullard ([58]). Moreover, the extended LB approach is also able to reproduce different droplets arrangements at the outlet channel of a microfluidic focuser, thus permitting to simulate soft flowing crystals at the scale of the actual microfluidic device.
To this purpose, two additional non-dimensional parameters (S σ and S κ ) have been introduced which measure the strength of inertia and surface tension versus the repulsive near-contact interactions. We found that, S κ = 1 acts as a natural threshold, above which the repulsive near-contact forces are no longer able to withstand the impact kinetic energy and prevent the coalescence between colliding fluid interfaces.
Even though in this paper we have discussed the specific case of a flow-focuser microfluidic device, the method presented in this work is expected to apply to a much broader variety of engineering and biomedical problems.
An application where our methods can be specifically useful (and in future will be developed for) is the design of systems where droplets would undergo an internal transition from viscous solutions to elastic materials. This in-droplet gelation "freezes" the shape of microparticles as it is at the outlet channel. The resulting microparticle systems (also produced in microfluidic devices) are employed e.g. for the encapsulation of living cells in hydrogels (biological reactors or sensors for toxicological screening) or of pharmaceutically active compounds in polymer matrices (controlled drug release). By applying our methods to such materials, it will be possible to rationally design complex, micro particle-based flowing crystals, where morphology/aspect ratio and any ensuing properties are precisely and topologically controlled. This would allow e.g. a) a permanent and tailored modulation of optical properties orthogonally to the channel direction, producing soft waveguides, or b) a very fine control of the conditions of dynamic arrest (macroscopic 'gelation') of the emulsion or foam obtained through the microfluidic device, which can be particularly useful in applications of 3D printing.
We would like to stress that the possibility of employing a mesoscale instead of a full scale description crucially relies upon the universality of the underlying physics or, in other words, its dependence on suitable dimensionless parameters measuring the relative strength of the interactions, rather than on the microscopic details of the interactions themselves. The aim of this work was precisely to present a coarse-grained approach encompassing the basic physical features of near contact interactions. In this regard, the proposed model represents an upscale of the interactions occurring at the interface level. Whilst being a dramatic simplification of the underlying physics at the molecular level, the results obtained in this paper suggest that, at least at the spatial scale at hand, a coarse grained description is appropriate to describe the mesoscale evolution of an interacting multidroplet system. To conclude, it is hoped that the numerical approach presented here may open the way to an experimental-scale modeling of soft flowing crystals, promising new chances to decode the complexity which characterize this fascinating state of flowing matter.