Gang Wang

Pierrette Guichardon

Umberto D'ortona

Improved pseudopotential lattice Boltzmann model to prevent droplet coalescence: from emulsion to dry foam

In this paper, we propose a droplet non-coalescence method based on the lattice Boltzmann (LB) pseudopotential model. Due to its simplicity, the pseudopotential model has been widely used to simulate multiphase flows. However, with the original pseudopotential model, droplets in contact tend to merge, which poses a limitation in simulating densely packed droplets without coalescence, as observed in emulsions and foams. To overcome this limitation, our aim is to develop a noncoalescence method that does not introduce additional sets of distribution functions to simulate a stable system of droplets in close contact. In this proposed method, all droplets must be labeled with unique indexes, enabling the identification of each droplet, and pseudopotential repulsions are applied at droplet contacting interfaces to prevent them from merging. Although this approach effectively prevents droplet coalescence, there is a small amount of mass transfer between droplets.

To correct this, we also incorporate a mass control scheme. It is found that the proposed method is able to model stable highly concentrated suspensions (HCS) with or without density and viscosity ratios. Besides, the proposed method can also simulate flowing HCS with arbitrarily high volume fractions of the dispersed phase allowing to simulate closely packed microfluidic soft flowing crystals or foams.

I. INTRODUCTION

The lattice Boltzmann method (LBM) has become an efficient and powerful tool to simulate fluid flows [START_REF] Krüger | The lattice Boltzmann method[END_REF][START_REF] Chen | [END_REF]. LBM can deal with complex geometries [3][4][5][6][7][8][9], flow in porous media [8][9][10][11][12][13][14][15], non-Newtonian fluids [15][16][17][18][19][20][21][22][23][24][25], is easy to implement in parallel computers or graphical processing units [26][27][28][29][30][31][32][33][34][35] and is extremely efficient to simulate multi-phase and multi-component flows [36][37][38][39][40][41][42][43][44][45][46][47][48][49]. As a consequence, LBM has been widely used to simulate microfluidic flows [36,[50][51][52][56][57][58][59][60][61].

Microfluidic emulsification has been considered as a promising technique because it can ensure the production of monodispersed droplets [62][63][64], contrarily to conventional emulsification in a batch that usually produces polydisperse droplets. The main drawback of microfluidics is the extremely low flow rates that are involved. Therefore, it is of the utmost importance to find the optimal operating conditions. In the case of the production of an emulsion, the objective is to find the optimal geometry and inlet flow rates such that at the outlet, the dispersed phase made of micro-droplets reaches the highest possible concentration in the continuous phase. Extremely high concentrations can be achieved in experiments such that droplets are in close packing and form "soft flowing crystals" [65][66][67][68]. With these high droplet concentrations, the produced emulsions are often metastable and the addition of a surfactant is necessary to prevent the merging of droplets, coarsening, and the loss of the monodisperse size distribution [65]. From a numerical point of view, the situation may be even more critical. Indeed, in multi-component LBM simulations, when two droplets of the same component come in contact, they immediately merge to form a larger droplet. This is inherent in LBM methods where interfaces have a finite width, typically a few lattice spacings, as opposed to interface-tracking Lagrangian techniques [53,54]. This phenomenon is unwanted to study numerically emulsions stabilized by a surfactant. As a consequence, numerical studies of T, X, and Y junctions in microfluidic devices usually restrict to low concentrations of the dispersed phase where droplets flow far away from each other.

A simple method has been proposed to prevent the coalescence of droplets of the same fluid and reach higher dispersed-phase concentrations [49,52,60]. Each time a new droplet is created, its color is changed such that it is different from the continuous phase and the neighboring droplet ones. This approach is interesting since it creates perfectly stable emulsions, but it suffers from several drawbacks. As pointed out by the authors, 5 to 6 different colors are necessary for 2D systems to prevent the neighboring of two droplets of the same color. In the LBM, this requires 5 to 6 different sets of the distribution function, which implies a larger amount of computer memory. In a 3D case, each droplet may have up to 12 first neighbors for a monodisperse system, and the required memory will be even larger. This problem can be solved by dynamically deallocating the memory associated to colors that are not present at the node [60], requiring more programming skills. A second problem arises if the droplets have relative displacements, for example being in a shear flow or crossing a channel width reduction. The colors should be dynamically assigned again to prevent any coalescence.

Another proposed approach adds to the short-range interaction (repulsion between different phases or attraction for each phase) that ensures immiscibility between the discrete (droplets) and the continuous phase, a mid-range repulsion between droplets that prevents them from touching and merging [52,54,55,69]. The main advantage is this method is that only two components are modeled which strongly reduces the memory require-ment. Essentially applied to the LB color-gradient model [41][42][43], this approach has successfully modeled complex emulsion structures in microchannel [55] that convincingly reproduces experiments [70].

The objective of this work is to propose a model that combines the advantages of the two previous approaches, i.e. a two-phase model that only uses two components to have a low memory requirement but where the newly created droplets are perfectly stable and will never merge with other droplets. In this work, the LB pseudopotential model is preferred since it is very popular due to its simplicity of implementation. In the pseudopotential model, the mesoscopic interactions of pseudopotentials (based on component densities) are used to represent the microscopic intermolecular interactions and the separation of different components is a result of the mesoscopic repulsions [47,48]. The extension to other LB multicomponent models, free-energy model [37][38][39], color-gradient model [41][42][43] and phase-field model [44][45][46] should cause any peculiar difficulties. The second objective of this model is to have a wide range of applications, from isolated drops to a highly concentrated suspension where the continuous phase occupies only a few percent of the volume, from isodensity phases to a high density difference to be closer to a liquid-gas system and to phases having a viscosity difference to model water-oil emulsions.

The article is organized as follows. In Sec. II, we introduce the LB pseudopotential multicomponent model and the proposed non-coalescence method. Section III tests the stability of contacting droplets with or without viscosity and density ratios and shows how the proposed method can model flowing HCS, foams and study flows in microfluidic devices. Finally, concluding remarks are given in Sec. IV.

II. NUMERICAL METHOD

A. The LB pseudopotential multicomponent model

The pseudopotential multicomponent model is widely used because of its mesoscopic nature and simplicity. This model can be used to simulate single-phase (liquidliquid) or multiphase (liquid-gas) flows. In this article, both flow systems are considered, The dispersed droplet (bubble) phase and the continuous fluid phase (suspending fluid) are represented by d and c, respectively. The lattice Boltzmann equation for each component σ can be given by [START_REF] Krüger | The lattice Boltzmann method[END_REF]:

f (σ) i (x+c i ∆t, t+∆t) = f (σ) i (x, t)+Ω (σ) i ∆t+S (σ) i ∆t (1)
where the Bhatnagar-Gross-Krook (BGK) collision operator Ω (σ) i is given by:

Ω (σ) i = - f (σ) i (x, t) -f eq(σ) i (x, t) τ m (x) (2)
The equilibrium distribution function f eq(σ) i (x, t) is defined by [71]:

f eq(σ) i (x, t) = ρ (σ) ω i 1 + c i • u c 2 s + (c i • u) 2 2c 4 s - u • u 2c 2 s (3)
where ω i are discrete weights for velocity sets c i . In this article, we utilize D2Q9 and D3Q19 velocity sets respectively for 2D and 3D simulations. The speed of sound is c s = 1/ √ 3. The fluid density for each component ρ (σ) , the total fluid density ρ, and the barycentric velocity u are defined by:

ρ (σ) = i f (σ) i (x, t), ρ = σ ρ (σ) , (4)
and

u = 1 ρ σ i f (σ) i c i + F (σ) tot ∆t 2 . ( 5
)
To simulate the flow with a large viscosity ratio, the system mixture relaxation time τ m is used in this paper. τ m is defined by [72,73]:

τ m (x) = σ ρ (σ) (x)ν (σ) ρ(x)c 2 s + 1 2 ∆t. (6)
When the viscosity ratio is 1, the system mixture relaxation time reduces to the traditional relaxation time with a constant value over the simulation domain. ν (σ) is the kinematic viscosity for the σ component, and it is given by:

ν (σ) = c 2 s (τ (σ) -0.5∆t) (7)
where τ (σ) is the relaxation time for the σ component. The Guo forcing scheme is utilized to compute the source term S (σ) i [START_REF] Krüger | The lattice Boltzmann method[END_REF]74]:

S (σ) i = ω i 1 - ∆t 2τ m (x) c i -u c 2 s + (c i • u)c i c 4 s F (σ) tot (8) where F (σ)
tot is the total force on each component σ, which consists of four forces. They are respectively the gravitational force F

F (σ) tot = F (σ) g + F (σ) inter + F (σ) intra + F (σ) s . (9)
The gravitational force is given by:

F (σ) g = ρ (σ) g. ( 10
)
Two situations are considered in this paper, a low gravity where droplets migrate, and get into contact but remain mostly spherical in a similar way to a wet foam, and a strong gravity where droplets become polygonal in a pattern analogous to a dry foam or a highly concentrated emulsion [75,[START_REF] Cantat | Foams: Structure and Dynamics[END_REF].

In the pseudopotential multicomponent model, the repulsive interaction is applied between different components to create the component separation: (11) where ψ (σ) (x) is the pseudopotential density of the component σ, and it can be given in many forms [START_REF] Krüger | The lattice Boltzmann method[END_REF]. One simple form of the pseudopotential density is ψ (σ) (x) = ρ (σ) (x), and this form is used in all the single-phase simulations throughout this paper. Another commonly used form is ψ (σ) (x) = ρ 0 (1 -exp(-ρ (σ) (x)/ρ 0 )), which bounds ψ (σ) between 0 and ρ 0 to avoid the appearance of too large interaction, and this form is used in all the multi-phase simulations throughout this paper. G σσ is the interaction strength between different components, and is chosen to be positive to ensure the repulsion between different components.

F (σ) inter (x) = -ψ (σ) (x)G σ σ i ω i ψ (σ) (x + c i ∆t)c i ∆t
For a multiphase system, the intra-component interaction force is added to induce a density ratio between different phases. (12) where ψ (σ) (x) has same forms as in Eq. (11). G σσ is the interaction strength and it is negative to mimic the attraction between elements of the same component.

F (σ) intra (x) = -ψ (σ) (x)G σσ i ω i ψ (σ) (x + c i ∆t)c i ∆t
For the boundary condition, the bounce-back method [START_REF] Ginzbourg | [END_REF][78][79] is usually used to achieve the no-slip boundary condition. The wetting boundary condition is realized by applying the pseudopotential interaction between the fluid and the wall. A virtual density is given to the wall, and the fluid-solid interaction force is given by:

F (σ) s (x) = -G σs ψ (σ) (x) i ω i s(x + c i ∆t)ψ s(σ) c i ∆t (13
) where s(x) is a switch function, and its value is 1 for solid nodes and 0 for fluid nodes. G σs is the adhesive strength between the fluid and the wall. ψ s(σ) is the virtual solid density for the component σ, which is a constant value. One can change both G σs and ψ s(σ) to adjust contact angles [80]. A recently developed curved boundary condition for the pseudopotential flows is also used in this paper, and details can be found in Wang et al. [7].

Originally, discrete weights ω i used in calculation of F are the same as that in Eq. (3), and it retains fourth-order isotropy of the interactions. To reduce spurious currents and increase the stability of the system, the calculation of the interactions can easily be extended to eighth-order isotropy by including the interactions with next-nearest fluid neighbors [81,82]. At the walls, next-nearest virtual solid neighbors are also required to apply the wetting wall boundary condition with eighth-order isotropy [7].

B. Non-coalescence of droplets in contact

A simple method is developed to prevent the coalescence of droplets in contact. In this method, all the droplets must be assigned unique indexes throughout the simulation. A pseudopotential-type repulsive force is applied at the region where droplets with different indexes are in contact. We will see that with this approach, the mass of each droplet is not perfectly constant, but the material is exchanged between droplets in an analogous way to Ostwald ripening occurring in foam or in suspensions [START_REF] Cantat | Foams: Structure and Dynamics[END_REF]. Indeed, in the pseudopotential model, the dispersed phase is partly miscible in the continuous phase allowing material exchange between droplets. In some situations this mass variation is unwanted and a simple method to compensate for mass transfer is proposed, by slightly adjusting the inter-component interaction strength G σ σ at the droplet region.

Droplet labeling

Here, we present algorithms to label all droplets with unique color indexes. In the pseudopotential model, since there is always miscibility between different fluids, we set a minimum droplet density ρ min . When the density of the droplet component at a node is above this minimum value (ρ (d) (x) > ρ min ), this node is considered inside the droplet region, otherwise, it is considered inside the continuous fluid region. In this paper, the threshold value is ρ min = 0.02 for all liquid-liquid simulations except for the microfluidic simulations where ρ min = 0.05. For all liquid-gas simulations, ρ min = 0.015. Two situations should be considered for the droplet labeling. If the number of droplets does not change during simulation, only the migration of the droplets should be considered. If the droplet is leaving a LB node, the color index is set to 0 when the discrete phase density goes below ρ min . Contrarily, if a droplet is arriving in the LB node (x), the discrete phase density ρ (d) (x) increases, and when the density becomes larger than ρ min , the color index is assigned to the one having the largest ρ (d) among the first neighbors (4 in 2D and 6 in 3D). This algorithm is simple, efficient, and easy to parallelize. As all droplets are surrounded by a layer of 3 to 4 nodes of the continuous phase, no inconsistency in the color indexing appears. The second situation is when new droplets may appear. A search algorithm should be implemented to check that all non-continuous regions have different color indexes, and if not, the color index of one of the regions must be reassigned. Several algorithms, recursive or not, exist for such a problem. One is proposed in appendix A. These algorithms are generally complex, slow, and difficult to parallelize.

Non-coalescence repulsion

To prevent droplets in contact from merging, the noncoalescence repulsion is applied at the contacting interface region. Here we only demonstrate how to apply the non-coalescence repulsion in 2D for simplicity, but it can be easily extended to 3D as well. For each node at droplet diffuse interfaces, we search for its midrange neighbors to find the nodes located inside neighboring droplets with different indexes, and then pseudopotential-type repulsions are applied between them. In this work, we do not search for all neighbors, instead, only one layer of midrange neighbors (third or fifth layer) is searched in order to save computational cost, and it is enough to simulate stable suspensions (see Sec. III A). Figure 1 shows the two midrange layers used in this paper. The third layer contains all nodes with a distance 2 √ 2, 3, or √ 10 to the center node. The distance for the fifth layer is 5, √ 26, or √ 29. Figure 2 illustrates the non-coalescence repulsion applied at the interface node x (filled red node) when two droplets are detected in contact. The red circles denote droplet nodes. The blue squares denote continuous fluid nodes. For the node x, its third layer of neighbors is searched, and three droplet nodes labeled with a different index are found. Thus, three repulsions are applied to the droplet phase at the node x. The total non-coalescence repulsion at the node x can be given as follows.

F (d) rep (x) = -ψ (d) (x)G rep i∈X l i ψ (d) (x + c i ∆t)c i ∆t (14)
where G rep is the strength of the repulsion. c i are discrete velocity sets for third layer of neighbors. Domain X contains all the nodes with a different droplet index on the third layer. l i are the repulsion weights for velocity sets c i , and they may be given in many forms. In this paper, it is given by:

l i = 1 |c i ∆t| 3 . ( 15
)
This choice is made to compensate for the difference in velocity length |c i |, where a larger |c i | corresponds to a smaller l i , ensuring that the non-coalescence repulsion from each neighbor is appropriately weighted. We have found that using a constant value for all l i also leads to a stable suspension but with a slightly higher mass transfer.

Control of mass transfer

Although the method described above can suppress the droplet coalescence, the mass of each droplet is not perfectly constant, but material is exchanged between droplets in an analogous way to Ostwald ripening occurring in foams or in suspensions [START_REF] Cantat | Foams: Structure and Dynamics[END_REF]. Indeed, in the LB pseudopotential model, the discrete phase is partly miscible in the continuous phase allowing material exchange between droplets. In some situations this mass variation is unwanted. To solve this, we propose to slightly adjust the inter-component interaction strength G σ σ at each time step: when a droplet loses mass, we slightly increase G σσ in the droplet region to acquire the mass of droplet component dissolved in the suspending fluid and vice versa. The slight change on G σ σ at node x is given by:

∆G(x) = k 1 - m(I(x)) m ref (I(x)) (16)
where m ref (I(x)) is the constant reference mass of the droplet labeled with index I(x) and m(I(x)) is the instantaneous mass of the droplet. k is the mass control strength, it is a simple scalar that controls the magnitude of ∆G.

The pseudopotential inter-component interaction is thus modified as follows.

F (σ) inter (x) = -ψ (σ) (x)(G σ σ + ∆G(x)) i ω i × ψ (σ) (x + c i ∆t)c i ∆t. ( 17
)

III. NUMERICAL RESULTS

A. Model stability validation

The objective of this test is to evaluate the stability of the proposed non-coalescence method by simulating confined highly concentrated suspensions (HCS) with size distributions under both low and high gravitational forces. We first test confined HCS with unit density and viscosity ratios and perform a parameter study (Sec. III A 1). Then, we show the effectiveness of the proposed non-coalescent method in simulating confined HCS with a viscosity ratio of around 41 (Sec. III A 2) and a density ratio of around 10 (Sec. III A 3). To launch a confined HCS simulation, we choose a lattice system with N x × N y = 560 × 560 (in lattice units). The bounce-back boundary condition [START_REF] Ginzbourg | [END_REF][78][79] is applied in all directions, and the wall is completely nonwetting to the droplet phase and completely wetting to the continuous phase. The interaction strength between different components is G σ σ = 4. The viscosity (ν (σ) ) is 1/6 for both components, resulting in a uniform τ m value of 1 throughout the simulation domain. Initially, 16 droplets are placed equidistant in the computational domain with 3 radii (50∆x, 60∆x, and 70∆x), and droplets with the same radius are represented by the same color (see Fig. 3). Inside the droplets, the initial fluid densities are ρ (d) = 1 and ρ (c) = 0.013, and the opposite densities for the outside of the droplets. We first let the droplets relax for 1 × 10 4 time steps, enabling the droplet masses to reach equilibrium with slight deviations from their initial values. Then a gravitational force is introduced in the upward vertical direction, acting on the dispersed droplet phase. Both high [g = 3 × 10 -5 (in lattice unit)] and low [g = 2 × 10 -6 (in lattice unit)] gravities are used in simulations. We also investigate the cases with (k > 0) and without (k = 0) mass control. In the case of mass control, the reference masses (m ref ) are chosen as the droplet masses at t = 10 4 ∆t. The third layer non-coalescence repulsion is employed to prevent droplets from merging. Figure 4 shows confined HCS obtained at t = 3×10 6 ∆t with the non-coalescence repulsion strength G rep = 5. In Fig. 4(a) and Fig. 4(c), droplets in contact remain spherical in shape since low gravity is applied, while in Fig. 4(b) and Fig. 4(d), droplets are strongly deformed and organized like a dry foam due to the high gravity. To quantitatively show the mass transfer between droplets without [Fig. 4(a) and Fig. 4(b)] and with [Fig. 4(c) and Fig. 4(d)] mass control, we plot the evolution of normalized masses of droplets A and B in Fig. 5, where droplets A and B (see Fig. 4) correspond to the droplets with the highest normalized mass (m max /m ref max ) and the lowest normalized mass (m min /m ref min ), respectively, after running simulations for 3 × 10 6 time steps, which is sufficiently long for the suspensions with mass control to reach a stable regime. Figure 5 shows that for both low gravity and high gravity cases with mass control (k = 20), droplet masses remain almost constant throughout simulations. However, for the low gravity case without mass control, there is mass change throughout the simulation and reaching a maximum of around 20% at the end. For the high gravity case without mass control, there is a sharp mass transfer at the initial stage when droplets are moving upwards rapidly due to the high gravity. Then the suspensions become stable and there is no mass exchange. We then study the influence of the non-coalescence repulsion strength G rep [Fig. 6(a)] and the mass control strength k [Fig. 6(b)] on mass transfer. The y-axes of both figures are kept identical to facilitate comparison. The mean absolute relative mass difference, defined as

E av = 1 16 16 i=1 (m i -m ref i )/m ref i ( 18
)
at t = 3 × 10 6 ∆t, is used as a measure of mass transfer. The non-coalescence repulsion strength G rep is tested in a confined HCS system without mass control, and its value ranges from 3 to 6. This range is determined to ensure that G rep is large enough to prevent coalescence in the high gravity case while avoiding system instability due to excessively large values. Figure 6(a) shows that mass transfer is not very sensitive to the value of G rep . The errors in low gravity cases are higher compared to those in high gravity cases. This is because for low gravity cases, there is mass transfer throughout the simulation and the suspensions do not have a stable regime, while for all high gravity cases, the suspensions rapidly reach a stable regime and mass transfer stops quickly. The mass control strength k is studied in a system with G rep = 5.

The range of k is set from 0 to 50, considering that excessively large k values lead to suspension instability, while negative k values accelerate mass transfer. Figure 6(b) shows that in the absence of mass control (k = 0), the mean error of mass is around 5% for the low gravity case and 1.5% for the high gravity case. The error decreases dramatically to around 0.03% for the low gravity case and 0.2% for the high gravity case when k is increased to 10. Further increase of k continues to reduce the error.

Remarkably, when k is increased to 50, the error is reduced to 0.01% for the low gravity case and 0.06% for the high gravity case ensuring an almost null mass transfer between droplets.

Viscosity ratio

In order to demonstrate the ability of the noncoalescence method in simulating stable confined HCS with a high viscosity ratio, we conduct a new test case using the same configurations as described in Sec. III A 1, except for the viscosity ratio. In this test, we set ν (d) = 1/600 and ν (c) = 1/6. By computing the mixture relaxation time τ m (x) through Eq. ( 6) at each time step, we obtain an effective viscosity ratio of approximately 1/41 between the dispersed phase and the continuous phase, which is large enough to model real microfluidic emulsions. The mass transfer, in this case, is equivalent to the case with unit viscosity and density ratios, and the implementation of the mass control scheme reduces it to almost zero. Figure 7 shows the viscosity contours of confined HCS under low gravity [Fig. 7 All units are in lattice units. See Videos 1 and 2 in Supplementary Material [START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF].

dispersed phase region. At droplet interfaces, there is a smooth transition of viscosity.

Density ratio

The non-coalescence method is also tested in a liquidgas system. Due to the stability issue of the pseudopotential multicomponent model with a very high density ratio, we use a density ratio of around 10 in this test. We adopt similar configurations as Sec. III A 1 with a few modifications. In this test, we set G σσ = 3.5 to induce the separation of different components and G σσ = -6 within the continuous-phase component, which leads to the phase separation. The initial fluid densities are ρ (d) = 0.27 and ρ (c) = 0.0041 inside the bubble, and ρ (d) = 0.0087 and ρ (c) = 2.91 outside the bubble. A low gravity [g = 2 × 10 -5 (in lattice unit)] and a high gravity [g = 2×10 -4 (in lattice unit)] are applied to the dispersed phase. We find that the third layer non-coalescence repulsion is not sufficient to prevent bubbles from merging, even though with a large repulsion strength G rep . To solve this problem, we employ the fifth layer repulsion, which allows the repulsion to be applied to more nodes at contacting interfaces. The repulsion strength G rep is set to 15. The mass transfer in this case is faster than in the case of unit viscosity and density ratios, but the mass control scheme still reduces it to almost zero. Figure 8 shows the density contours confined HCS (the density ratio is approximately 10) obtained at the stationary regime with mass control (k = 20).

B. From highly concentrated suspensions toward foams

The objective of this work is not to perform an extensive study of foams, but to show how the newly proposed model can simulate realistically foams and HCS. 2D and 3D foams and HCS are studied and typical properties of like yield stress and visco-elasticity are put in evidence.

As discussed in Sec. III A 3, with this simple LB model, a modest density ratio of around 10 is obtained, far from the density ratio of a real foam which is around 1000.

For simplicity, we nevertheless call these systems foams, while for a density ratio of 1, the name Highly Concentrated Suspensions is also used.

To create a realistic polydisperse foam before generating any flow, the following procedure is used. First, a simple regular pattern made of cells filled with the discontinuous phase and walls with the continuous phase is generated (see Fig. 9, t = 0). The thickness of the cell walls is chosen to control the volume ratio of the discrete to the continuous phases, allowing a continuous transition from a wet foam to a dry foam (see Fig. 10). Each cell is numbered and associated with a volume that is randomly generated. The sizes of the cells are uniformly distributed between d min and d max such that d max =2d min and such that the total surface fraction in 2D and the volume fraction is equal to F d , the wanted fraction of the discrete phase. The initial volume of each cell is then allowed to relax to its final volume. Two relaxation procedures have been tested. The first procedure uses the mass control method proposed in Sec. II B 3. This method is rather slow, typically 10 6 LB time steps, since it uses the natural miscibility of each component of the LB pseudopotential model and cells can only exchange material that has to diffuse through the continuous phase. In the second method, the node densities of each cell are multiplied by a factor (1 + ε (V (i) -V t (i)) /V (i)) where V t (i) and V (i) are the instantaneous and final volume of the i th cell respectively. ϵ is a factor chosen small enough such that the simulation is stable during the volume relaxation. With this procedure, the polydisperse foam is obtained in typically 15000 time steps (Fig. 9 and Video 5 in Supplemental Material [START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF]). After relaxation, the mass of each cell is controlled to keep all volumes un-t=0 t=2000 t=4000 t=15000

FIG. 9: Relaxation of a HCS to obtain a droplet size distribution. Time is expressed in the LB time step. The continuous phase is colored in blue, the discrete phase is colored in a random mixture of red and green.

The domain is a lattice of 1000 × 1000 and 100 cells (bubbles) are generated. See Video 5 in Supplementary Material [START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF].

changed using the mass control method of Sec. II B 3.

To allow the foam to flow, the periodic boundary condition is chosen in the horizontal directions (in 2D and 3D), and no-slip boundary conditions for the top and bottom walls. The partially saturated method is applied at the top and bottom walls, but as they are parallel to the lattice and placed equally spaced between two LB nodes, this simplifies to fullway bounce-back. During the relaxation of the foam, the wall surface properties are such that they equally attract each phase and as a consequence, the contact angle of all interfaces is nearly 90 • .

When the initial relaxation is finished, gravity is applied horizontally to induce a flow in the foam. To prevent experimentally a slip of the foam at the wall, a usual technique is to use grooved walls with grooves having a size comparable to that of the bubbles [START_REF] Katgert | [END_REF]. Numerically, this implies complex boundary conditions, and a simpler method is used. The wall surface properties of each solid node are modified such that the densities of both phases are equal to the ones of the closer fluid node. As can be seen from Videos 6 to 9 (see Supplemental Material [START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF]) or by comparing Figs. 9 and 10, the layers of cells in contact with the wall do not slip at all. A gravity equal to g = 1.5 10 -6 (in lattice unit) is applied horizontally towards the right on both phases in Fig. 10. For the 3 lowest volume fractions of the discrete phase, a flow is obtained, while for F d = 0.92, the applied force is not strong enough and after a few oscillations, the system remains frozen (see Videos 6 and 7 in Supplemental Material [START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF]).

Figure 11 shows the time evolution of the HCS mean velocity for various continuous phase volume fractions F c . Time and velocity are normalized using the equivalent single-phase flow whose mean velocity evolves like v 0 (1 -exp(-t/t 0 ), where v 0 is the stationary mean velocity and t 0 is the typical time of flow establishment. As expected, higher is the continuous phase volume fraction F c , higher is the mean velocity. Nevertheless, even for volume fraction as large as F c ≃ 30%, the mean velocity remains small, around 10% of that of the equivalent Figure 12 reports the stationary mean velocity as a function of the continuous phase volume fraction. The curve shows a threshold around F c = 0.08 where no flow is obtained for the applied gravity (g = 1. 5 10 -6 ). This confirms the yield stress nature of HCS and foam [START_REF] Cantat | Foams: Structure and Dynamics[END_REF].

The yield stress nature of the HCS can also be apprehended from the velocity profiles presented in Fig. 13. All velocity profiles present a plug flow at the center of the channel. In a Poiseuille flow, the highest shear is at the side walls while the shear decreases towards the center of the channel to be exactly null at the center-line. The widths of the plug flows increase with decreasing F c , indicating an increase of the yield stress with the volume fraction of the discontinuous phase. This can be seen in Fig. 10 where the cells near the top and bottom walls are stretched and elongated while in the center of the channel, the cells are more circular for a wet case, and polygons are more "regular" for a dry HCS. The effect of the applied driving force is now consid-ered. The previous HCS with a discontinuous phase volume fraction of F d = 0.9 is used with increasing gravity force ranging from g = 10 -6 to 2 10 -6 (expressed in lattice units). Figure 14 shows the time evolution of the flow mean velocity for the various gravity forces that are applied instantaneously at t = 0. Time and velocity are reduced like in the previous section. Indeed, as the equivalent single-phase velocity depends on the applied gravity, the case with g = 1.5 10 -6 is used as a reference. Two different behaviors are obtained, above a threshold (g > 1.25 10 -6 ), a flow is obtained and stronger is the gravity force, larger is the mean velocity. For lower applied forces, after a transient phase due to a reorganization of the HCS, the mean velocity oscillates and decays.

For g = 1.25 10 -6 , a nice sinusoidal oscillation is obtained indicating that the applied force is very near the threshold and confirming the visco-elastic nature of HCS and foams. Larger gravity forces have been tested, but they rapidly induce the breakage of cells near the endwalls. For the case g = 2 10 -6 in Fig. 14, the simulation was stopped at t/t 0 = 0.5 since for longer times, cell breakages happen and the mean velocity goes on increasing since a HCS with smaller cells flows more easily. Figure 15 shows the velocity profiles for the various applied gravity forces that induce a flowing HCS. As expected, stronger is the gravity force, larger is the velocity. The yield stress nature of the HCS can again be apprehended from the plug flow at the center of the flow whose width reduces when the applied gravity increases. More surprising is the negative velocities obtained near the side walls. To our knowledge, such negative velocities have never been reported in flowing HCS nor in flowing foams. To understand the origin of these negative veloc-ities, the velocity vector field near a side wall (y = 0) is presented (Fig. 16) and shows a recirculation flow inside the next wall cells. This LBM model allows to capture the flow inside each cell but also the flow in the continuous phase that separate the cells. FIG. 17: Examples of 3D HCS flowing in a channel with an applied gravity g = 1.5 10 -6 for two lattice sizes: (a) a 600 3 lattice and (b) a 300 × 300 × 600 lattice. See Videos 8 and 9 in Supplementary Material [START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF].

Figure 17 presents the equivalent 3D flowing HCS for two lattice sizes, (a) a 600 3 lattice and (b) a 300 × 300 × 600 lattice. We recognize the usual structure of 3D HCS and foams including bubbles, faces, edges (plateau borders) and vertices [START_REF] Cantat | Foams: Structure and Dynamics[END_REF]. To generate a flowing 3D HCS, the same procedure as in 2D is applied: relaxation of the bubble volumes such that d max = 2d min (see Video X in S.M.), coating of the top and bottom wall to prevent any slip and application of a horizontal gravity force g = 1.5 10 -6 on both phases at t = 0 that induces a flow thanks to the x-y periodic boundary conditions (see Videos Y and Z in S.M.). The mean cell size is around 60 lattice nodes giving 1000 cells for the large lattice (600 3 ) and 250 cells for the small one (600×300×300).

Figure 18 shows the time evolution of the flow mean velocity for various gravity forces and the two studied lattice sizes. Time and velocity are reduced using the equivalent single-phase flow. The large lattice (600 3 ) is only simulated for a reduced time (t/t 0 ≃ 0.1) since it requires about 300 Gbytes of memory and has run for two weeks on a 48-core computer. It nevertheless shows that the smaller lattice gives very close results compared to the large one. Like in the 2D case, a threshold is obtained in the applied gravity force below which no flow occurs and only decaying oscillations are observed, confirming the yield stress viscoelastic nature of a HCS. When gravity is increased, a flow occurs and as expected, larger is the gravity, larger is the flow. Velocity profiles are not reported since they are very similar to 2D cases with a central plug flow and negative velocities near the walls.

Changing density, toward a real foam

As already stated, the used multiphase LBM model is too simple to obtain a density ratio comparable to those observed experimentally in real foams. A modest density ratio of 10 is obtained allowing to show some differences between a HCS and a foam. Flows of foams with a density ratio of 10 have been simulated and they give velocity profiles and video extremely similar to HCS with no density ratio. As a consequence, they are not reported here.

Differences are observed when comparing T1 transformations during the relaxation process in the foam and in the HCS [Figs. 19(a)-(f)]. The previous 2D 1000×1000 lattice with 100 cells is used again, but with 4 walls that are non-wetting to the discrete phase. The continuous phase volume fraction is F c = 0.1. The initial square lattice of cells rapidly relaxes to an irregular lattice (See Video XX and YY in Supplementary Materials). The sizes of the cells are modified extremely slowly like in Sec. III B to reach d max = 2d min , but instead of previously 15000 time steps, the relaxation is performed in t 0 = 5 10 5 time steps. This relaxation period is chosen long enough such that T1 transformations are distinguishable over time. A relaxation rate is computed (g) Time evolution of the relaxation rate for the foam and the HCS.

according to

R = 1 N δt N i f d i (t + δt) -f d i (t) (19)
where f d i (t) is the fraction of the discrete phase at the node i at time t. The sum is made over all the N lattice nodes and δt, the time between two pictures, is taken equal to δt = 1000 time steps. R measures the rate of modification between two systems distant of 1000 time steps. Figure 19(g) reports the time evolution of the relaxation rate where the relaxation time t 0 is used to rescale the time axis. After a rapid relaxation occurring before t/t 0 ≃ 0.2, R decays exponentially but with T1 transformations that are visible as peaks in the relaxation rate curve. The peaks indicated by arrows correspond to the T1 transformations illustrated in Fig. 19(a)-(f). T1 transformations occur two to three times more rapidly for a foam with a density ratio of 10 compared to the HCS without a density ratio. Indeed, for T1 transformation occurring in foam, liquid films are reorganized while the gas in each bubble plays no role. In a HCS, the liquid in the adjacent cells is moving as well, slowing down the process due to inertia and the liquid viscosity. Even though a modest density ratio of 10 is achieved with this LB multiphase model, and differences in T1 transformation in foams and HCS can be put in evidence. These few examples demonstrate how efficient is this LB pseudopotential model with the proposed non-coalescent scheme to study foams and HCS. The proposed non-coalescence method is also applied to model concentrated oil/water emulsions through a microfluidic flow-focusing device, whose geometry is shown in Fig. 20. This geometry is the same as the one used by Du et al. [64], and it has a dispersed inlet, two continuous inlets that are symmetric, and an enlarged outlet. Droplet formation in this channel occurs through the periodic pinch-off of the dispersed phase by the continuousphase flow. This pinch-off process is influenced by the combined effects of interfacial tension, viscous forces, and inertial forces. In the simulations, we adopt a 3D lattice system with N x × N y × N z = 400 × 90 × 30 (in lattice units). Since the channel has a curved geometry, the recently proposed improved partially saturated method [7], which is a curved boundary condition for the LB pseudopotential flows, is used to simulate the wetting wall boundary condition. At the dispersed inlet, the wall is completely wetting to the dispersed phase and completely non-wetting to the continuous phase. At other parts of the channel, the wetting condition is the reverse. The diameter of the inlet channels is 30∆x, and each lattice spacing ∆x responds to 10/3 µm in the experiment.

The interaction strength between different components is G σσ = 4. To achieve a viscosity ratio ϕ ν consistent with the experiment, we set the viscosity ν (d) = 0.479 and ν (c) = 0.087, which gives an effective viscosity ratio ϕ ν of around 5.1 between the dispersed and continuous phases [64]. At the three inlet layers, we impose the Poiseuille velocity profile. Outlet boundary conditions are complex for a two-phase flow and there is a lot of ongoing work on that. For simplicity, we impose at the outlet the velocity profile obtained for a single phase flow and simply adjust the total flow rate to be exactly equal to the inlet flow rate. Two layers of ghost nodes next to the outlet layer are updated at each time step in order to calculate the pseudopotential interaction at the outlet. The ghost densities are given by solving the convective boundary condition with the average velocity at the outlet layer [85]. With this procedure, we notice that the movement of droplets leaving the numerical domain is almost not perturbed and the previous droplets are not perturbed at all (see Video 12-15 in Supplementary Material [START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF]). We set non-coalescence repulsion strength G rep = 5, and there is no need to apply mass control since droplets leave the short outlet channel quickly and mass transfer between droplets is almost negligible. f) to better show the spatial arrangement of droplets inside the channel. To numerically reproduce the experimental results, we ensure that Reynolds numbers (Re = U D/ν) of the dispersed and continuous inlet flows in our simulation match the values used in the experiments. This ensures that the flow characteristics, such as the velocity and viscosity, are consistent between the numerical and experimental setups. Additionally, we also match the capillary numbers (Ca = ρνU/γ) of both the dispersed and continuous phases to those used in the experiments. This ensures that the interfacial tension effects and the relative importance of viscous and surface tension forces are accurately represented in the simulation, allowing for a meaningful comparison with the experimental results. In Fig. 21 As a result, the droplet size decreases, which is attributed to the higher viscous shear force exerted by the continuous-phase flow, which makes it easier to pinch off the droplets.

To further demonstrate the effectiveness of the proposed non-coalescence method in simulating microfluidic concentrated emulsions, we adopt the same microfluidic channel to simulate a system where the continuous phase has a significantly higher viscosity compared to the dispersed phase (ϕ ν ≈ 1/13). This is equivalent to a system where oil is the continuous phase. This large viscosity contrast enables the continuous phase to readily pinch off the dispersed phase even at relatively low flow rates. Conversely, the dispersed phase exhibits a higher flow rate, causing the droplets to occupy a larger portion of the channel and become closely packed, even with deformations in their shapes. Figure 22 shows 3D snapshots of obtained highly concentrated emulsions. In Fig. 22(a), the flow rate ratio between the dispersed and continuous phases is ϕ Q = 20/3. The droplets are closely packed in a single row, almost completely filling the channel. They exhibit strong deformations, resembling cubic-like shapes. We then maintain the flow rate of the dispersed phase and decrease the flow rate of the continuous phase by half (ϕ Q = 10/3), resulting in the formation of smaller droplets self-assembling into two rows [Fig. 22(b)]. The droplets deform into triangular shapes since they are in close contact.

IV. CONCLUSION

To conclude, we have proposed a droplet noncoalescence method based on the LB pseudopotential model. By assigning unique indexes to each droplet and applying pseudopotential repulsion forces, we have effectively prevented droplet coalescence. In situations where mass exchange occurs, we have introduced a mass control method to compensate for such transfers. We conducted simulations of confined HCS with different viscosity and density ratios under low and high gravitational forces to evaluate the stability of the proposed non-coalescence method. We have shown that the method can prevent droplet coalescence, and stable suspensions can be obtained when the mass control scheme is implemented. The proposed non-coalescence method has been successfully applied to model flows of highly concentrated suspensions, foams, and concentrated oil/water emulsions in a microfluidic flow-focusing device. To reproduce the experimental results numerically, we matched Re and Ca numbers in our simulations to that of the experimental setup. The simulation results show agreement with experimental observations. Additionally, the method is further validated by simulating a system with a large viscosity contrast, where the continuous phase readily pinches off the dispersed phase even at a low flow rate to form highly concentrated emulsions. The obtained emulsions exhibit close packing of droplets and deformations in their shapes. These findings demonstrate the effectiveness of the non-coalescence method in simulating various types of microfluidic concentrated emulsions.

FIG. 1 :FIG. 2 :

 12 FIG.1:The third and fifth layers of lattice neighbors of the red node at the center for a 2D case.

1 .

 1 FIG. 3: The initialization of confined HCS where droplets with the same radius have the same color. All units are in lattice units.

20 FIG. 4 :

 204 FIG. 4: Confined HCS with the unit viscosity and density ratio between the dispersed and continuous phases. In all cases, G rep = 5. In (a) and (c), g = 2 × 10 -6 (in lattice unit). In (b) and (d), g = 3 × 10 -5 (in lattice unit). In (a) and (b), there is no mass control (k = 0). In (c) and (d), the mass control strength k = 20.

FIG. 5 :

 5 FIG. 5: Normalized masses of droplet A (m max /m ref max ) and droplet B (m min /m ref min ) with respect to t.

FIG. 6 :

 6 FIG. 6: The mean absolute relative error of the droplet masses at t = 3 × 10 6 ∆t with respect to (a) the non-coalescence repulsion strength G rep and (b) the mass control strength k.

  FIG. 7: The viscosity contours of confined HCS obtained at the stationary regime. The viscosity ratio between the dispersed and continuous phases is approximately 1/41, and the density ratio is unity. Both cases have mass control with k = 20. (a) Low gravity g = 2 × 10 -6 and (b) high gravity g = 3 × 10 -5 .All units are in lattice units. See Videos 1 and 2 in Supplementary Material[START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF].

FIG. 8 :

 8 FIG. 8: The density contours of confined HCS obtained at the stationary regime. The density ratio between the dispersed and continuous phases is approximately 1/10, and the viscosity ratio is unity. (a) Low gravity g = 2 × 10 -5 and (b) high gravity g = 2 × 10 -4 . All units are in lattice units. See Videos 3 and 4 in Supplementary Material [83].

1 .

 1 From wet to dryIn this first series of 2D simulations, the influence of the volume fractions of the discrete F d and continuous F c phases is highlighted. F d is in the range [0.71; 0.92] and F c = 1 -F d in the range [0.08; 0.29]. As the initial density and viscosity ratios are equal to one, the system is a Highly Concentrated Suspension. Figure10shows 4 typical HCS of increasing discrete volume fraction where the geometry of the discontinuous phase evolves from round shapes to polygonal cells.

F d = 0. 71 F d = 0. 8 FFIG. 10 :

 71810 FIG. 10: Examples of HCS flowing in a channel, four dispersed phase volume fractions are considered, from F d = 0.71 to 0.92. Pictures are taken when the flows have reached a stationary regime t/t 0 ≃ 0.3. See Videos 6 and 7 in Supplementary Material [83].

FIG. 11 :FIG. 12 :

 1112 FIG.11: Time evolution of the mean velocity for different continuous-phase volume fractions F c and a gravity g = 1.5 10 -6 . v 0 is the mean velocity of the equivalent single-phase flow.

FIG. 13 :

 13 FIG.13: Velocity profiles for different continuous phase volume fractions F c . The applied gravity is g = 1.5 10 -6 and v 0 is the center-line velocity of the equivalent single-phase flow.

6 g = 2 10 - 6 v/v 0 FIG. 14 :

 66014 FIG.14: Time evolution of the mean velocity for different applied gravity g and a continuous phase volume fraction of F c = 0.1. v 0 is the mean velocity of the equivalent single-phase flow for g = 1.5 10 -6 .

- 6 FIG. 15 :

 615 FIG. 15: Velocity profile for different applied gravity. The continuous phase volume fraction is F c = 0.1. v 0 is the center-line velocity of the equivalent single-phase flow for g = 1.5 10 -6 .

FIG. 16 :

 16 FIG. 16: Velocity field near the side wall for a HCS with F c = 0.1 and a gravity g = 1.75 10 -6 .

6 W 6 W 6 v/v 0 FIG. 18 :

 666018 FIG.18: Time evolution of the mean velocity for different applied gravity g and for the two studied lattice sizes. v 0 is the mean velocity of the equivalent single-phase flow for g = 2.5 10 -6 .

10 FIG. 19 :

 1019 FIG. 19: Example of T1 transformations (indicated by thick arrows) in a foam (a) to (c) and in a HCS (d) to (f), see Videos XX and YY in Supplementary Materials.(g) Time evolution of the relaxation rate for the foam and the HCS.

FIG. 20 :

 20 FIG. 20: The geometry of the microfluidic flow-focusing device. All units are in lattice units

Figure 21 presents

 21 concentrated emulsions obtained from experiments and 3D simulations. Snapshots at the x-y midplane [Fig. 21(b) and Fig. 21(e)] are compared with experimental results [Fig. 21(a) and Fig. 21(d)]. Additionally, 3D snapshots of dispersed droplets are shown in Fig. 21(c) and Fig. 21

  (

  (a)-(c), the concentrated emulsion is obtained with Re d = 0.21, Ca d = 1.3 × 10 -2 , Re c = 0.45, and Ca c = 1.08 × 10 -3 . The droplets maintain spherical shapes and have a relatively large size. They are in close contact with each other and organized in a staggered way. In Fig. 21(d)-(f), Re c and Ca c are, respectively, increased to 1.51 and 3.59 × 10 -3 while Re d and Ca d remain the same.

FIG. 21 :FIG. 22 :

 2122 FIG. 21: Comparison of the concentrated emulsions obtained by the experiments [(a) and (d)] and the 3D simulations [(b), (c), (d), and (f)]. In (a)-(c), Re d = 0.21, Ca d = 1.3 × 10 -2 , Re c = 0.45, and Ca c = 1.08 × 10 -3 . In (d)-(f), Re d = 0.21, Ca d = 1.3 × 10 -2 , Re c = 1.51, and Ca c = 3.59 × 10 -3 . See Videos 12 and 13 in SupplementaryMaterial[START_REF]= 2 × 10 -5 towards a stationary state. Video 4: Same for g = 3 × 10 -4 . Video 5: Relaxation of a HCS to obtain the desired size distribution. Video 6: Flow of a HCS having a continuous phase volume fraction Fc = 0.29. Video 7: Same for Fc = 0.1. Video 8: Relaxation and flow of a 3D HCS in a 600 3 lattice and g = 2.5 10 -6 . Video 9: Same for a 600 × 300 × 300 lattice. Video 10: Slow relaxation of a foam with a density ratio ρc/ρ d = 10. Video 11: Same for a HCS with no density[END_REF].

ACKNOWLEDGMENTS

We would like to thank Dr. Jiupeng DU for performing microfluidic experiments. Centre de Calcul Intensif d'Aix-Marseille University is acknowledged for granting access to its high performance computing resources.

Appendix A: A droplet labeling algorithm

In the algorithm, D is the simulation domain. I old is the droplet index field for the previous time step t -∆t, and I is the field of droplet indexes for the current time step t. At the initial time step, unique indexes for all droplets should be assigned to I old , and then, I old is updated from I at the beginning of the labeling process. Starting from I old , we use the same indexes to label droplets in I. As a result, the indexes for each droplet remain constant throughout the simulation. However, when new droplets are generated from the breakup, new indexes will be assigned to them.

Algorithm 1 An algorithm to label each droplet with a unique index for a 2D case

if (I old (x, y) > 0 & ρ (r) (x, y) ≤ ρmin) then 9:

if (indexused(I(x, y)) = 0) then 13:

index ← I old (x, y) if

if (I(m + 1, n) = 0) then 33:

isearch ← isearch + 1


Enrico Mingo Hoffman

Silvio Traversaro
email: silvio.traversaro@iit.it

Alessio Rocchi
email: alessio.rocchi@iit.it

Mirko Ferrati
email: mirko.ferrati@centropiaggio.unipi.it

Alessandro Settimi
email: alessandro.settimi@centropiaggio.unipi.it

Francesco Romano
email: francesco.romano@iit.it

Lorenzo Natale
email: lorenzo.natale@iit.it

Antonio Bicchi
email: antonio.bicchi@centropiaggio.unipi.it

Francesco Nori
email: francesco.nori@iit.it

Yarp Based Plugins for Gazebo Simulator

Keywords: Simulator, Robotics, YARP, Gazebo, Open-Source

This paper presents a set of plugins for the Gazebo simulator that enables the interoperability between a robot, controlled using the YARP framework, and Gazebo itself. Gazebo is an open-source simulator that can handle different Dynamic Engines (ODE, DART, Bullet, SimBody), backed up by the Open Source Robotics Foundation (OSRF) and supported by a very large community. Since our plugins conform with the YARP layer used on the real robot, applications written for our robots, COMAN and iCub, can be run on the simulator with no changes. Our plugins have two main components: a YARP interface with the same API as the real robot interface, and a Gazebo plugin which handles simulated joints, encoders, IMUs, force/torque sensors and synchronization. The robot model is provided to the simulator by means of an SDF file, which describes all the geometric, dynamic and visual characteristics of a robot. Once the SDF is read from Gazebo, our plugins are loaded and associated with the simulated robot model and the simulated world. Different modules for COMAN and iCub have been developed using Gazebo and our plugins as a testbed: joint impedance control plus gravity compensation, dual arm Cartesian control for manipulation tasks, dynamic walking, etc. This work has been developed as part of a joint effort between three different European Projects "WALKMAN", "CoDyCo" and "SoftHands" aiming at implementing a common simulation platform to develop and test algorithms for our robotic platforms. This work is available as open-source to all the researchers in the YARP community (https://github.com/robotology/gazebo_yarp_plugins).

Introduction

In the past years, robotics researchers have been developing many robotics frameworks such as OpenRDK [START_REF] Calisi | Open-RDK: a modular framework for robotic software development[END_REF]), YARP [START_REF] Metta | Yarp: Yet another robot platform[END_REF]) or ROS [START_REF] Quigley | Ros: an open-source robot operating system[END_REF]) in order to ease the creation of generic applications for robots and encourage code reuse. The performance overhead introduced by these frameworks is balanced by the architectural benefits, for example they allow to build modular systems to execute one or more assigned tasks.

In these frameworks, the simulator is a module that represents the real robot at the interface level (Figure 1). Such simulator module accepts control input (desired joint torques, desired joint position, ...) and outputs sensory feedback (cameras, joint positions, ...) from the simulated world. These simulators usually allow to have the human in-the-loop permitting to train a human operator. The most important aspect is that they allow to develop modules that directly will work in the real robot without any need to rewrite code. In fact, when the real robot is used, there is a module that replaces the simulator by providing the same hardware interfaces.

By accurately simulating robots and environments, code designed to operate on a real robot can be executed and validated on the simulated equivalent system. This avoids common problems associated with With these concepts in mind we decided to extend one of the most known robotics simulator, Gazebo [START_REF] Koenig | Design and use paradigms for gazebo, an open-source multirobot simulator[END_REF]), to be compatible with one of the most used robotics framework, YARP (Figure 2), developed in the Italian Institute of Technology. YARP is supported by the iCub simulator (iCubSim, [START_REF] Tikhanoff | An open-source simulator for cognitive robotics research: The prototype of the icub humanoid robot simulator[END_REF]) that is dedicated to a specific platform. The needs of a more generic tool for simulating different robots rise up. Gazebo, which has been recently chosen as the simulator for the DARPA Virtual Robotic Challenge (VRC, DARPA [2013]), allows the use of different dynamic engines, it is easily expandable through plugins and it has a strong and active community. Gazebo is maintained by the Open Source Robotics Foundation [START_REF] Osrf | Open source robotics foundation[END_REF]).

This paper is organized in the following sections: State of Art, discusses some of the most popular simulation environments in robotics, Structure introduces in detail the Gazebo plugins developed in this work, Conclusions summarize the outcome of this effort and finally Future Works discusses the follow up activities.

State of Art

A large number of simulators have been developed in the past two decades [START_REF] Ivaldi | Tools for dynamics simulation of robots: a survey based on user feedback[END_REF]). Such simulators range from dynamic solver libraries to complex simulation environments/systems. The latter are usually large projects that provide both rigid body dynamic simulations and tools such as graphical editors, planner libraries, visualization tools, controllers and so on.

The Open Dynamics Engine (ODE, [START_REF] Smith | Open dynamic engine[END_REF]) is one of the most widely used rigid body dynamics engine in robotics simulation. ODE simulates chains of rigid bodies connected and constrained by different types of joints. It has a built-in collision detection system and implements hard contacts using nonpenetration constraint whenever two bodies collide. Beside the large number of project that use it, at the moment the development has been paused. Bullet [START_REF] Erwin | [END_REF]) is another dynamic engine. It implements different direct/inverse rigid body dynamic algorithms (eg. Featherstone articulated body algorithm, Featherstone [2007]) as well as different solvers (eg. Mixed Linear Complementarity Problem, MLCP) and contact models. Bullet is used for a wide range of projects and its community is active and continues to improve it constantly.

OpenRAVE [START_REF] Diankov | Automated Construction of Robotic Manipulation Programs[END_REF]) provides an environment for testing, developing, and deploying motion planning algorithms in real-world robotics applications. The main focus is on simulation and analysis of kinematic and geometric information related to motion planning. It provides many command line tools to work with robots and planners, and the run-time core is small enough to be used inside controllers and bigger frameworks. Industrial robotics automation is an important target application.

Webots [START_REF] Michel | Cyberbotics ltd. webots tm : Professional mobile robot simulation[END_REF]) is a development environment used to model, program and simulate mobile robots. With Webots the user can design complex robotic setups, with one or several, similar or different robots, in a shared environment. A large choice of simulated sensors and actuators is available. The robot controllers can be programmed with the built-in IDE or with third party development environments. The robot behavior can be tested in dynamic simulated worlds (ODE based). The controller programs can optionally be transferred to commercially available real robots.

V-REP [START_REF] Rohmer | V-rep: A versatile and scalable robot simulation framework[END_REF]), similarly to Webots, embeds different tools that permit fast developing of algorithms, the code can be transferred inside real robotic hardware.

Gazebo [START_REF] Koenig | Design and use paradigms for gazebo, an open-source multirobot simulator[END_REF]) is a multi-robot simulator for outdoor environments. As Stage (part of the Player project, [START_REF] Gerkey | The player/stage project: Tools for multirobot and distributed sensor systems[END_REF]), it is capable of simulating a population of robots, sensors and objects. It generates realistic sensor feedback and physically consistent interactions between objects. It includes an accurate simulation of rigid-body physics and allows the user to select between multiple dynamics engines (ODE, Bullet, SimBody [START_REF] Sherman | Simbody: multibody dynamics for biomedical research[END_REF] and DART Tech [2013]). Gazebo has been used to compare algorithms for navigation and grasping in a controlled environment.

Finally, two notable softwares are the OpenHRP project used in Japan for the HRP series [START_REF] Kanehiro | Openhrp: Open architecture humanoid robotics platform[END_REF]) and MuJoCo [START_REF] Todorov | Mujoco: A physics engine for model-based control[END_REF]) used for model-based control.

Our decision to add a YARP interface to Gazebo is motivated by the following considerations. We want to switch between fast, not accurate simulations and slow, accurate ones, thus we need the capability of choosing among different dynamic engines. We also want a simulator which is both easy to use and to expand in order to add new robot models. Finally, we prefer an open-source software with an active community and money investments. Gazebo fulfills our requirements, in particular it is expandable with a plugin structure: in this work our YARP interface is a collection of Gazebo plugins.

Structure

It is useful to understand Gazebo plugins and YARP device drivers before describing the structure of our plugins (from now on gazebo yarp plugins).

Gazebo plugins are C++ classes that extend the functionalities of Gazebo, while YARP device drivers are C++ classes used in YARP for abstracting the functionality of robot devices. Usually, each class of gazebo yarp plugins embeds a YARP device driver in a Gazebo plugin.

Gazebo Plugins

A plugin is a piece of code compiled as a shared library and inserted into the simulator. A plugin has direct access to all the functionalities of Gazebo from the physics engine to the simulated world. Furthermore, plugins are self-contained routines that are easily shared and can be inserted and removed from a running system. There are 4 types of plugins in Gazebo: world, model and sensor plugins are attached to and control a specific simulated world/model/sensor respectively, while system plugin is specified on the command line and loads during the Gazebo startup.

YARP Device Drivers

In YARP, a device driver is a class that implements one or more interfaces. There are three separate concerns related to devices in YARP:

• Implementing specific drivers for particular devices • Defining interfaces for device families • Implementing network wrappers for interfaces For example the Control Board device driver implements a set of interfaces that are used to control the robot (IPositionControl, ITorqueControl, etc.) and another set of interfaces to read data from the motors (IEncoders, etc).

Gazebo-YARP Plugins

The gazebo yarp plugins is made of:

• Gazebo plugins that instantiate YARP device drivers, • YARP device drivers that wrap Gazebo functionalities inside the YARP device interfaces.

The plugins/devices already implemented are the Control Board, 6-axis Force Torque sensor, Inertial Measurement Unit (IMU) and a Clock plugin used for synchronization. The first three plugins are directly related to the simulated objects and sensors, while the last one is a system plugin that synchronizes all the other YARP modules with the simulation time.

Control Board

The Control Board plugin allows to control the robot using YARP Interfaces, it is implemented as a Gazebo Model plugin. Every control board allows the user to control one or more joints (a kinematic chain such as the arm or leg, etc.) as specified in a configuration file. For each controlled joint the control board opens different interfaces, permitting the use of different type of controllers for each joint. Such interfaces include position control, torque control, encoders reading, torque The values read from the simulator are broadcasted through YARP interfaces in the YARP network, in a similar way the desired joint values come from YARP interfaces (Figure 4). The following YARP interfaces are used to control the robot.

• IPositionControl: a position control with a linear trajectory generator considering a max joint speed • IPositionDirect: a position control using Gazebo position PIDs • ITorqueControl: a perfect torque follower • IImpedanceControl: a joint impedance control with the following law

τ d = -P d (q -q d ) -D d q + τ of f set (1)
where q d is the desired equilibrium position, P d is the desired joint stiffness and D d is the desired joint damping. τ of f set is an extra signal that can be used for gravity compensation or inverse dynamics control.

Furthermore, the Control Board implements the ICon-trolMode interface that allows to change the type of controller online. All these interfaces are also available on the robot and they have the same behaviour.

6-axis Force/Torque sensor

A Force/Torque sensor measures a wrench in the robot structure (Figure 3). The sensor, at the time of writing, is simulated in Gazebo as if it was attached to the reference frame associated to a joint. On the YARP side, the reading of a generic sensor is implemented as a IAnalogSensor interface (Figure 5). The broadcasted data is a vector of six numbers representing the forces and the torques applied on that reference frame.

IMU sensor

An IMU measures velocity, orientation, and gravitational forces, using a combination of accelerometers and gyroscopes, of the link where it is placed. It is also possible to add white Gaussian noise on the measurement (Figure 7). Similar to the Force/Torque sensor, it is implemented as a IAnalogSensor interface.

Clock

A fundamental aspect in simulations is the synchronization between YARP modules and the simulated robot. A YARP module is a process in which one or more threads are started. When such modules are used in the real robot, the thread rate is timed by the machine (system) clock, also called the wall clock.

When the simulation is running we want the rate of such modules to be synchronized with the simulated time, otherwise the control loop could run faster or slower with respect to the simulated robot dynamics.

The real-time factor (RTF ) of the simulation is given by

RT F = update f requency × step time (2)
and is kept to one when the desired update frequency is the inverse of the time increased at each step in the simulation. For instance if the simulation runs with a real time factor of 0.1, 10 seconds are needed to simulate 1 real second. Within this situation, the controller process should also be slowed down 10 times to be coherent with the simulation. To solve this issue we developed a clock plugin that synchronizes modules with the simulated time. The clock plugin is implemented as a System plugin and publishes on a YARP port the time information from the simulator. For every simulation step, the simulation time is incremented and the timestamp is sent via socket. All the YARP threads implemented as control threads (which need to be run at a desired rate, i.e. YARP's RFModule or RateThread ) are automatically synchronized using the simulation clock if the YARP CLOCK environment variables is set, or if the module explicitly asks to. The yarp::os::Time functionalities are also transparently working using the wall-clock or the simulation clock depending on the environment variable. Thread sleeps are performed using the right wall or simulated time. When synchronized with the simulation clock the yarp::os::Time delay does not explicitly sleep on a wall clock, rather a scheduler is synchronized with the simulation clock by performing blocking reads on the YARP CLOCK port. This scheduler wakes up the threads that required a delay just once, when they have slept for the desired duration. Compared to the ROS::Time implementation which uses small sleeps on wall clock to check synchronization with the simulated clock, this allows to run simulations both slower and faster than real time and still have synchronization between threads and controls. In any case, when accessing the simulated clock Experiments showed the approach to be successful in synchronizing 1kHz control loops against simulations running 1kHz, thus having a 1ms clock granularity.

Simulation Description Format (SDF)

Gazebo uses an XML-style format, Simulation Description Format (SDF), to save and load information about a simulated world or model. An SDF encapsulates all the necessary information for a simulation such as:

• Scene: ambient lighting, sky properties, shadows.

• Physics: gravity, time step, physics engine.

• Models: collection of links, collision objects, joints, and sensors. • Lights: point, spot, and directional light sources.

• Plugins: world, model, sensor, and system plugins.

Our Control Board, Force Torque sensor and IMU plugins are included inside the SDF file that describes our robots. For our humanoid bipedal robot, COMAN, we have five Control Board plugins (one for each kinematic chain), four Force/Torque sensor plugins (two in the legs and two in the arms) and one IMU sensor plugin (placed on the back of the waist). The SDF descriptions of COMAN and iCub are available in https://github.com/EnricoMingo/ iit-coman-ros-pkg and in https://github.com/ robotology-playground/icub_gazebo respectively. The clock plugin is loaded trough a command line parameter when the simulator is started.

Conclusions

In this work we have presented a set of Gazebo plugins, named gazebo yarp plugins, that allow to connect the robotics framework YARP to Gazebo itself. Gazebo was chosen since it is easy to use, it has the possibility to switch between different rigid multi-body dynamics engines, it is Open-Source and has an active community. Our plugins are based on YARP device drivers in order to have exactly the same interfaces in the real and simulated robot. This allows to write modules that will work both in the simulator and in the real robot without the need to change the code. This is a very important paradigm in robotics research and develop since it minimizes the presence of errors due to code porting. Furthermore the simulator becomes a tool that helps the developer in testing and validation before using the real platform. Such plugins consist in: a Control Board plugin to control the robot, a Force Torque sensor plugin and an IMU plugin. A special plugin dedicated to synchronization between modules and simulator was also implemented. The plugins were tested to simulate two humanoid bipedal robots, the COMAN and the iCub, both from the Italian Institute of Technology.

Future Works

Gazebo yarp plugins is a project at an early stage that is gaining more and more interest inside the YARP community. Beside the good results obtained up to now, some works are still missing in order to be able to have 100% compatibility with all the Gazebo functions. Furthermore we still need to implement plugins to connect YARP device drivers dedicated to cameras and RGB-D sensors to the simulated ones in Gazebo.

We are also interested in multi-robot and human-robot simulation: we already have the possibility to easily simulate different robot models but it is still difficult to simulate multiple instances of the same robot. Furthermore, since our robots in IIT have flexible joints, we are investigating on how to simulate flexible joints without specifying extra joints/links inside the SDF of the robot. Finally we are planning an official release of our plugins inside the Gazebo community.

  (a) Module A connected to the simulator (b) Module A connected to the robot

Figure 1 :

 1 Figure 1: Module A write desired joint position and read actual joint position without knowing if it is interfaced with the simulator or the robot since they expose the same interface. hardware such as short battery life, hardware failures, and unexpected and dangerous behaviors, particularly during the initial stages of development and tuning of new modules and controllers. It is also much faster to have a simulation engine up and running than using a real robot, especially when the simulation engine can run faster than real-time. In this way the simulator becomes a fundamental part of the framework and the robot software development cycle as the first step to validate algorithms, thus minimizing the risks of hardware breaks.

Figure 2 :

 2 Figure 2: COMAN and iCub interacting inside a Gazebo simulation of a kitchen. Blue dots represent contact points.

  Plot of forces measured at the Force/Torque sensor placed on the right arm (data logged during simulation). Forces along x,y and z are respectively in red, green and blue Plot of torques measured at the Force/Torque sensor placed on the right arm (data logged during simulation). Torques along x,y and z are respectively in red, green and blue

Figure 3 :

 3 Figure 3: A Gazebo simulation running COMAN interacting with a debris

Figure 4 :

 4 Figure 4: Control Board plugin for the left arm kinematic chain. yarp::IPositionControl interface has a method positionMove() that can be used to set joint values inside a YARP module. The plugin implements such interface by calling the Publish() method inside the Gazebo API to move the simulated joints at each OnUpdate().

Figure 5 :

 5 Figure 5: The Force/Torque sensor in the left arm is implemented as a YARP IAnalogSensor interface. At every step the internal state of the plugin is updated with the last readings of forces and torques from the simulation.

Figure 6 :

 6 Figure 6: Time elapsed between each execution of the control loop, measured in simulation clock time and in wall clock time. Desired thread rate is 1kHz and simulation time step is 1ms

  Plot of accelerations measured by the IMU (data logged during simulation). Accelerations along x,y and z are respectively in red, green and blue

Figure 7 :

 7 Figure 7: iCub inside Gazebo.

The research leading to these results has received funding from the European Union Seventh Framework Programme [FP7-ICT-2013-10] under grant agreements n.611832 WALKMAN, ERC Advanced Grant no. 291166 SoftHands and the CoDyCo project (FP7- ICT-2011-9, No. 600716).