# Lattice Boltzmann investigation of droplet inertial spreading on various porous surfaces

Xavier Frank, Patrick Perre, Huai-Zhi Li

► To cite this version:

Xavier Frank, Patrick Perre, Huai-Zhi Li. Lattice Boltzmann investigation of droplet inertial spreading on various porous surfaces. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2015, 91 (5), 10.1103/PhysRevE.91.052405. hal-01186414

HAL Id: hal-01186414

https://hal.science/hal-01186414

Submitted on 14 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Lattice Boltzmann investigation of droplet inertial spreading on various porous surfaces

Xavier Frank*

IATE, INRA-CIRAD-UMII-SupAgro, 2 place Pierre Viala, 34060 Montpellier, France

Patrick Perré

École Centrale Paris, LGPM, Grande Voie des Vignes, 92290 Châtenay-Malabry, France

Huai-Zhi Li

Laboratoire Réactions et Génie des Procédés, Université de Lorraine, CNRS, 1 rue Grandville, BP 20451, 54001 Nancy Cedex, France

(Received 12 June 2014; published 26 May 2015)

The spreading of liquid drops on solid surfaces is a wide-spread phenomenon of both fundamental and industrial interest. In many applications, surfaces are porous and spreading patterns are very complex with respect to the case on smooth surfaces. Focusing on the inertial spreading just before the Tanner-like viscous regime, this work investigates the spreading of a low-viscosity droplet on a porous surface using lattice Boltzmann numerical simulations. The case of a flat surface is first considered, and it reveals a dependence on the solid equilibrium contact angle $\theta_s^{eq}$, which is in good agreement with published experimental data. We conducted numerical experiments with various surfaces perforated by a regular pattern of holes of infinite length. The results show that the global spreading dynamics is independent of the porosity morphology. Through the assumption that, for wetting, the pores can be regarded as surface patches with a contact angle of $\theta_{pore}^{eq} = 180^{\circ}$, we deduce an effective equilibrium contact angle $\theta_{eff}^{eq}$ on the porous surface from the Cassie-Baxter law. A spreading model is then proposed to describe both a prefactor and an exponent that are similar to a flat surface whose equilibrium contact angle is $\theta_{eff}^{eq}$. This model compares satisfactorily with a large number of numerical experiments under varying conditions.

DOI: 10.1103/PhysRevE.91.052405

PACS number(s): 68.08.–p, 47.55.dr, 47.55.nd, 47.56.+r

## I. INTRODUCTION

The spreading of liquid drops on solid surfaces is of both fundamental and industrial interest. In many applications, surfaces are porous or covered with a thin porous layer. The presence of a porous layer modifies the wettability of the substrate and, hence, the spreading patterns [1,2]. The spreading of liquid drops on porous media is a wide-spread phenomenon [3] involved in natural and industrial processes such as ink jet printing [4], adhesion, coating [5,6], lubrication, detergency, plant treatment, composite manufacturing, painting, and oil recovery. This topic has received much attention from researchers [7–13]. To the best of our knowledge, the understanding of the underlying mechanisms of drop spreading on a porous surface remains limited, in particular with respect to spreading on a smooth surface. Earlier investigations were mainly devoted to viscous wetting, and most modeling and simulation approaches were performed within the framework of a continuous porous material, pore-scale approaches being quite rare [14–16]. Long-term spreading dynamics, arising from a competition between capillarity and viscosity both inside and outside the porous material, are preceded by an initially faster spreading. The spherical cap-shaped droplet approximation fails to describe this short-term spreading.

Recently, both the inertial spreading dynamics and capillary invasion of pores were investigated by means of numerical simulations [15]. However, the dependence of spreading dynamics on the porosity remains quite qualitative in this study. The present paper aims to address this point with the

help of lattice Boltzmann (LB) numerical experiments. We will focus on simple porous substrate models, in which pore space is made of parallel holes of infinite length.

## II. NUMERICAL METHOD

The LB method is a computational fluid dynamics approach derived from nonequilibrium statistical physics [17–20], and it has been successfully applied to multiphase flows and wetting problems since the early 1990s [21–25]. In the LB framework, the fluid is described at the particle kinetics level. Particle probability density functions (PPDFs) are defined, from which hydrodynamic quantities can be deduced. The Navier-Stokes equation can be recovered from the LB equation [18,26,27]. A no-slip boundary condition is imposed using the well-known bouncing-back method, in which incoming particles are reflected back toward the direction from which they came. Among multiphase LB approaches, pseudopotential methods are appealing due to their simplicity, and they have provided valuable results in various fields, such as droplet dynamics [28,29], capillary filling [25,30,31], porous media [32], microfluidics [33,34], colloidal fluids [35], contact line dynamics [36], complex fluid-fluid interfaces [37], and wetting phenomena [38–41].

In the present study, gas-liquid flows are simulated with the help of the well-known Shan-Chen pseudopotential approach [42]. The velocity used to compute equilibrium distributions is modified according to a local fluid-fluid force that emerges from interactions between neighboring fluid particles. A liquid phase and a vapor phase emerge spontaneously within this framework. Liquid density $\rho_L$, vapor density $\rho_V$, and surface tension $\sigma$ depend on the magnitude of fluid attraction

<footnote>* xavier.frank@supagro.inra.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Scheme of a droplet spreading numerical experiment. A nonporous surface is used as an example. A: perspective view of the initial state. The solid surface is in the xy plane, and both the xz and yz symmetry planes are shown. B1: definition of the droplet radius R. B2: definition of the spreading radius r.</figcaption></figure>

in the Chen-Shan model. An interface is located at the middle of the transition from $ρ_G$ to $ρ_L$, where density $ρ$ reaches the value $ρ_M = (ρ_G + ρ_L)/2$.

In this study, Martys and Chen's fluid-solid pseudopotential [43] describes interactions between fluid particles and the solid. The equilibrium contact angle $θ_s^{eq}$ depends on the magnitude of both fluid-fluid and fluid-solid local interactions in the Chen-Shan model. Please see the previous works [24,42,43] for more details about the numerical approach.

### III. NUMERICAL SETUP

We employ symmetries through the xz and yz planes to reduce the demand for numerical resources (Fig. 1, panel A). As a consequence, the initial position of the droplet center has to be placed on the z axis. The initial z position of the droplet center of mass is fixed to place the bottom of the liquid phase at the middle of z range. Simulation box dimensions are fixed to 180 × 180 × 450, and the initial droplet diameter is set to 190. Simulations are conducted without gravity, and the wettability of the solid is defined using the solid equilibrium contact angle $θ_s^{eq}$. First, the droplet is equilibrated during 10 000 LB iterations without any solid node in the simulation box. The equilibrated droplet radius R, which is slightly different from the initial radius, has to be measured at the final state of this equilibrating phase (Fig. 1, panel B1). Then, the solid nodes are introduced in the simulation box, with a solid position adjusted in such a way that the droplet interface slightly touches the solid phase. To compute the radius of the droplet's footprint r (Fig. 1, panel B2), a density profile is extracted from a z,x = const. line, just above the solid surface. This specific position is chosen to extract the density profile only above the solid, and not above the pores. Wetting simulation runs during 15 000 iterations, and the spreading radius r is computed and saved at each time step.

### IV. SPREADING ON FLAT SURFACES

First, numerical spreading experiments are conducted with flat surfaces. As spreading dynamics can depend significantly on liquid viscosity $μ_L$ [44], it is crucial to compare the magnitude of inertial and capillary forces with that of viscous

forces. In the present work, the value of the Laplace number $La = \frac{ρ_L σ R}{μ_L^2}$ is $La = 135$ (in other words, the Ohnesorge number is $Oh = 1/\sqrt{La} = 8.6 \times 10^{-2}$), compared with the La values in the Bird et al. study [45], which ranged from 130 to 72 000. As $La \gg 1$ we can then consider that the droplet is of low viscosity in our simulations. In such a case, the three following steps can be identified within droplet spreading dynamics [46]: a first step, independent of wettability, in which the spreading radius scales as $r \approx t^{1/2}$, a second step, in which the spreading dynamics depend on the solid equilibrium contact angle [45–47], and, finally, a slow viscous step [48]. The first step is not observed in our numerical simulations, as a consequence of grid resolution limits, and only the intermediate regime is studied in the present work. Under such conditions, the inertial dynamics of droplet spreading is dominated by the propagation of a capillary wave from the bottom to the top of the droplet interface [45]. As a consequence, the time scale $τ$ of inertial spreading is $τ = \sqrt{ρ_L R^3}/σ$. The normalized spreading radius $r/R$ is deduced from r using R as the length scale, and it is plotted as a function of $t/τ$ for 5 values of the equilibrium contact angle $θ_s^{eq}$ (Fig. 2).

Strictly speaking, only the first inertial regime ($t/τ \lesssim 0.04$) can be characterized by a power law [46], and the intermediate regime ($0.04 \lesssim t/τ \lesssim 2$) is a crossover between the early inertial spreading and the viscous regime. However, as shown by Bird et al., a power law is a good approximation for low viscosity

<figure><img src="image_3.png" /><figcaption>FIG. 2. Normalized spreading radius on a flat surface as a function of dimensionless time for various surface equilibrium contact angles. Each curve is fitted using a power law [Eq. (1)]. Insets are power law parameters as functions of the contact angle. Left inset: power law prefactor, right inset: power law exponent. In each inset, empty markers are experimental results from Bird et al. [45], full markers are values deduced from our numerical simulations, and lines are models of the power law prefactor and exponent, $C(θ_s^{eq}) = 1.468 - 0.00631 \times θ_s^{eq}$ and $α(θ_s^{eq}) = 0.485 - 0.00102 \times θ_s^{eq}$, respectively. The theoretical duration of inertial spreading is underlined with a vertical line.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Pore lattices used. C-S: square lattice of circle-sectioned pores, C-CS: centered square lattice of circle-sectioned pores, S-S: square lattice of square-sectioned pores. A unit cell is underlined with dashed lines. Both the pore period $λ$ and pore size $d_{pore}$ are shown for each case.</figcaption></figure>

droplet spreading dynamics in this intermediate phase [45]. This power law approximation is expressed as

$$ \frac{r}{R} = C \left(\frac{t}{\tau}\right)^{\alpha}, \qquad (1) $$

where $α$ is the power-law exponent and $C$ is the prefactor, which is in good agreement with the experimental results [45]. The evolution of the power law prefactor $C(\theta_s^{eq})$ and exponent $α(\theta_s^{eq})$ with $\theta_s^{eq}$, shown in Fig. 2 as insets, has to be underlined and compared with earlier experimental results. Both decrease with $\theta_s^{eq}$, which is in agreement with the experiments from Bird et al. [45]. The duration of inertial spreading [49] can be easily calculated from $τ$ as $T \approx (\rho_L \sigma R / \mu_L^2)^{1/8} τ$. Bird et al. showed that the duration of the inertial regime does not depend significantly on the equilibrium contact angle [45]. In the present study, the evaluation leads to $T \approx 1.84τ$, which is in good agreement with our numerical simulations (Fig. 2).

## V. SPREADING ON POROUS SURFACES

We consider now surfaces perforated by a regular pattern (square or centered square lattice) of holes of infinite length. Between neighboring cylindrical pores, solid walls are effective to avoid any liquid communication. Surface porosity is defined as $ε = S_{pores}/λ^2$; $S_{pores}$ is the sum of the transversal sections of all pores belonging to one lattice cell, and $λ$ is the lattice period. The higher is the porosity, the lower is the driving force of spreading. As a consequence, the spreading dynamics are expected to slow down as $ε$ increases. However, the porosity of a surface can be distributed in various ways: small pores with high numerical density or larger pores with lower numerical density, or various pore lattices and various pore section shapes. Do these details matter or not? If this is not the case, the spreading dynamics will depend only on $ε$ regardless of the pore section shape, the pore lattice unit cell pattern, or the pore lattice period. To clarify this point, porous surfaces are built up using various geometrical parameters for both square and centered square lattices of circle-sectioned pores, and the square lattice of square-sectioned pores (Fig. 3).

For each kind of lattice, the geometry of the porous surface is completely defined by the pore size $d_{pore}$ (diameter for circle-sectioned ones and side for square-sectioned ones) and spatial period $λ$. Porosity $ε$ is easily computed from $d_{pore}$ and $λ$. Surface parameters are summarized in Table I. As an example, a perspective view of one solid phase is shown in Fig. 4.

Numerical spreading experiments are conducted in the same way as for flat surfaces. The surface porosity leads to

<table><caption>TABLE I. Parameters of pore lattices: pore lattice type, pore size $d_{pore}$, and pore period $λ$, which is deduced the surface porosity $ε$. Both $d_{pore}$ and $λ$ in lattice units (l.u.).</caption><thead><tr><th>Pore lattice</th><th>d<sub>pore</sub></th><th>λ</th><th>ε</th></tr></thead><tbody><tr><td>C-S</td><td>10</td><td>20</td><td>0.20</td></tr><tr><td>C-S</td><td>10</td><td>30</td><td>0.087</td></tr><tr><td>C-S</td><td>10</td><td>40</td><td>0.049</td></tr><tr><td>C-S</td><td>20</td><td>30</td><td>0.35</td></tr><tr><td>C-S</td><td>20</td><td>40</td><td>0.20</td></tr><tr><td>C-S</td><td>30</td><td>40</td><td>0.44</td></tr><tr><td>C-CS</td><td>10</td><td>30</td><td>0.18</td></tr><tr><td>C-CS</td><td>20</td><td>40</td><td>0.39</td></tr><tr><td>S-S</td><td>10</td><td>20</td><td>0.25</td></tr><tr><td>S-S</td><td>20</td><td>40</td><td>0.25</td></tr><tr><td>S-S</td><td>20</td><td>30</td><td>0.44</td></tr><tr><td>S-S</td><td>30</td><td>40</td><td>0.56</td></tr></tbody></table>

a greater complexity of the gas-liquid interface close to the solid surface: the solid-liquid-vapor triple line is no longer circular, but exhibits oscillations (Fig. 5). Liquid spreads on the solid surface and crosses the pores. Imbibition can start after the pore opening has been covered by the liquid interface. As expected, pores are invaded when $\theta_s^{eq} < 90°$, not when $\theta_s^{eq} > 90°$ (Fig. 6). Both droplet height evolution and droplet volume variation emerge from a competing mechanism between spreading on the solid phase of the surface and capillary imbibition inside the porous medium. Such a point was detailed in an earlier paper [15] and is not developed here.

The spreading dynamics are normalized as for flat surfaces, by means of $R$ and $τ$ (Fig. 7). Compared with the case of

<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) Perspective view of the solid phase. The square lattice of circular section pores; the pore diameter is $d_{pore} = 20$ l.u., and the pore lattice spatial period is $λ = 30$ l.u. A small part of the solid phase is removed to facilitate the visualization.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) Snapshot of a droplet spreading simulation on a porous surface. The square lattice of circular section pores; the pore diameter is $d_{pore} = 20$ l.u., the pore lattice spatial period is $\lambda = 30$ l.u., and the intrinsic solid equilibrium contact angle is $\theta_s^{eq} = 77^{\circ}$.</figcaption></figure>

a flat surface, successive accelerations and decelerations can be observed as the three-phase contact line crosses the pores. Clearly, the spreading radius dynamics are similar to a power law, despite these oscillating phenomena.

<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) Snapshot of droplet spreading simulation on a porous surface with two values of the intrinsic solid equilibrium contact angle, the solid phase being removed. The square lattice of circular section pores; the pore diameter is $d_{pore} = 20$ l.u., and the pore lattice spatial period is $\lambda = 30$ l.u.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Normalized radius of the wetted zone as a function of dimensionless time during droplet spreading on two porous surface pairs of porosity, $\epsilon = 0.20$ and $\epsilon = 0.44$, respectively, and with the solid equilibrium contact angles $\theta_s^{eq} = 41^{\circ}$ and $\theta_s^{eq} = 112^{\circ}$. Both $d_{pore}$ and $\lambda$ in lattice units (l.u.).</figcaption></figure>

As expected, spreading slows down when $\theta_s^{eq}$ increases for a given porosity, or when porosity increases, $\theta_s^{eq}$ being fixed (Fig. 7). Surprisingly, the global spreading dynamics do not depend on how the porosity is distributed over the surface. Increases in $\theta_s^{eq}$ and $\epsilon$ induce similar evolutions of spreading dynamics, and a global model, including both parameters, should be attainable.

The simplest approach to consider both the contact angle and morphology of a complex surface in a single law is to assume that the surface behaves as an equivalent smooth surface. Such a description leads to the famous Cassie-Baxter law [50]. Considering a composite surface involving two materials with two different equilibrium contact angles, the surface is supposed to exhibit an effective homogeneous equilibrium contact angle $\theta_{eff}^{eq}$, depending on both the intrinsic equilibrium contact angle and the proportion of basic materials. In the present case, we have a solid phase with an equilibrium contact angle $\theta_s^{eq}$ and a surface fraction $(1 - \epsilon)$, and a pore phase with equilibrium contact angle $\theta_{pore}^{eq}$ and surface fraction $\epsilon$. Using these assumptions, an effective equilibrium contact angle $\theta_{eff}^{eq}(\theta_s^{eq}, \epsilon)$ could be deduced from the Cassie-Baxter law as follows:

$$ \cos \theta_{eff}^{eq} = (1 - \epsilon) \cos \theta_s^{eq} + \epsilon \cos \theta_{pore}^{eq}. \quad (2) $$

Effective power-law parameters can stem from the computed values of $\theta_{eff}^{eq}$, as if the porous surface is an equivalent smooth solid surface with equilibrium contact angle $\theta_{eff}^{eq}$:

$$ C_{eff}(\theta_s^{eq}, \epsilon) = C(\theta_{eff}^{eq}), \quad (3) $$

$$ \alpha_{eff}(\theta_s^{eq}, \epsilon) = \alpha(\theta_{eff}^{eq}). \quad (4) $$

The application of the Cassie-Baxter law to our case requires $\theta_{pore}^{eq}$ to be fixed. A simple hypothesis can then be proposed:
<figure><img src="image_9.png" /></figure>

FIG. 8. Normalized spreading radius on various porous surfaces with various intrinsic solid equilibrium contact angles vs the prediction of the proposed effective spreading law with $\theta_{\text{pore}}^{\text{eq}} = 180^{\circ}$. Pore lattices are square lattice of circle-sectioned pores (C-S), centered square lattice of circle-sectioned pores (C-CS), and square lattice of square-sectioned pores (S-S). Both $d_{\text{pore}}$ and $\lambda$ in lattice units (l.u.). For each of the 12 pore lattices, 5 values of $\theta_s^{\text{eq}}$ were used (cf. Fig. 2), which means a total of 60 curves.

as long as a pore is not invaded by liquid, it remains an obstacle to global spreading. Thus, $\theta_{\text{pore}}^{\text{eq}} = 180^{\circ}$ seems to be a straightforward hypothesis.

To compare our effective law with the numerical results, we plot the ratio $r/R$ as a function of $C_{\text{eff}}(t/\tau)^{\alpha_{\text{eff}}}$, $t/\tau$ being a curve parameter (Fig. 8). Concretely, for each simulation, $\theta_{\text{eff}}^{\text{eq}}$ is computed as a function of both $\theta_s^{\text{eq}}$ and $\epsilon$, and effective power-law parameters $C_{\text{eff}}$ and $\alpha_{\text{eff}}$ are deduced from linear fits of values extracted from the numerical simulation of droplet spreading on a smooth surface (Fig. 2). Despite oscillations, it is clear that the spreading radius dynamics from numerical simulations are similar to a power-law whose prefactor and

exponent are predicted by the Cassie-Baxter law. Moreover, various curves fit this simple effective law well, although the proposed spreading law does not explicitly take into account the distributed pattern of porosity on the surface. Only the pore fraction matters whatever its distribution. The global behavior of the droplet is very close to that on an equivalent smooth surface.

Recently, Stapelbroek et al. conducted droplet spreading experiments with various complex solid surfaces, such as microtextured and chemically striped substrates [51]. In this study, the data collapse of the crossover time between the early inertial spreading and intermediate regime could be achieved using only a single parameter, the effective contact angle, regardless of the details of the substrate. Such result corroborates our findings. To overcome the energy barriers in the presence of topographic or chemical patterns that lead to a spatial variation of the local surface energy, an additional source of energy is required in the vicinity of the contact line. This energy must come from the kinetic energy of the flow inside the droplet. As a consequence, it could be reasonably argued that the balance between the kinetic energy of the spreading drop and the surface energies is a relevant mechanism, and that the details of a complex surface play a minor role in the inertial spreading global dynamics. Of course, such a picture is no longer valid in the viscous regime [47].

## VI. CONCLUSIONS

In summary, we have proposed a droplet spreading law on porous surfaces, including both the intrinsic solid contact angle and surface porosity. Lattice Boltzmann numerical simulations were performed for smooth surfaces with various equilibrium contact angles and the predicted power law for spreading dynamics compares favorably with published experimental results. Simulations with various porous surfaces and solid equilibrium contact angles show that increasing porosity slows down spreading in the same way as increasing the equilibrium contact angle for smooth surfaces. Assuming that pores are perfectly nonwetting patches, we deduced an effective equilibrium contact angle of the surface, and both the power law prefactor and exponent from an effective equilibrium contact angle are in satisfactory agreement with the simulation.

[1] H. Koivula, M. Toivakka, and P. Gane, J. Colloid Interface Sci. **369**, 426 (2012).

[2] V. M. Starov, S. R. Kostvintsev, V. D. Sobolev, M. G. Velarde, and S. A. Zhdanov, J. Colloid Interface Sci. **252**, 397 (2002).

[3] D. Quéré, Annu. Rev. Mater. Res. **38**, 71 (2008).

[4] A. Määttänen, P. Ihalainen, R. Bollström, M. Toivakka, and J. Peltonen, Colloids Surf., A **367**, 76 (2010).

[5] A. Zadražil, F. Stepanek, and O. Matar, J. Fluid Mech. **562**, 1 (2006).

[6] P. Colombo, G. Mera, R. Riedel, and G. D. Sorarù, J. Am. Ceram. Soc. **93**, 1805 (2010).

[7] P. Neogi and C. A. Miller, J. Colloid Interface Sci. **92**, 338 (1983).

[8] M. Denesuk, B. J. J. Zelinski, N. J. Kreidl, and D. R. Uhlmann, J. Colloid Interface Sci. **168**, 142 (1994).

[9] S. H. Davis and L. M. Hocking, Phys. Fluids **11**, 48 (1999).

[10] S. H. Davis and L. M. Hocking, Phys. Fluids **12**, 1646 (2000).

[11] V. M. Starov, S. A. Zhdanov, S. R. Kosvintsev, V. D. Sobolev, and M. G. Velarde, Adv. Colloid Interface Sci. **104**, 123 (2003).

[12] L. Cueto-Felgueroso and R. Juanes, Phys. Rev. Lett. **101**, 244504 (2008).

[13] T. G. D'Onofrio, H. K. Navaz, B. Markicevic, B. A. Mantooth, and K. B. Sumpter, Langmuir **26**, 3317 (2010).
[14] J. C. R. Neyval, R. F. Griffiths, and J. M. Santos, Appl. Math. Model. 32, 341 (2008).

[15] X. Frank and P. Perré, Phys. Fluids 24, 042101 (2012).

[16] S. Meng, R. Yang, J. S. Wu, and H. Zhang, Int. J. Heat Mass Transf. 77, 828 (2014).

[17] R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222, 0130008 (1992).

[18] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Clarendon Press, Oxford, 2001).

[19] D. d'Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, and L. S. Luo, Philos. Trans. R. Soc. London A 360, 437 (2002).

[20] C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. 42, 439 (2010).

[21] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, and F. Toschi, Phys. Rev. Lett. 97, 204503 (2006).

[22] M. L. Blow and J. M. Yeomans, Philos. Trans. R. Soc. A 369, 2519 (2011).

[23] X. Frank and H.-Z. Li, Phys. Rev. E 71, 036309 (2005).

[24] X. Frank, D. Funfschilling, N. Midoux, and H.-Z. Li, J. Fluid Mech. 546, 113 (2006).

[25] X. Frank, G. Almeida, and P. Perré, Int. J. Multiphase Flow 36, 599 (2010).

[26] Y. H. Qian, D. D'Humières, and P. Lallemand, Europhys. Lett. 17, 479 (1992).

[27] H. Chen, S. Chen, and W. H. Matthaeus, Phys. Rev. A 45, R5339 (1992).

[28] Q. Kang, D. Zhang, and S. Chen, J. Fluid Mech. 545, 41 (2005).

[29] J. Hyvälöma, A. Koponen, P. Raiskinmäki, and J. Timonen, Eur. Phys. J. E 23, 289 (2007).

[30] P. Raiskinmäki, A. Shakib-Manesh, A. Jäsberg, A. Koponen, J. Merikoski, and J. Timonen, J. Stat. Phys. 107, 143 (2002).

[31] F. Diotallevi, L. Biferale, S. Chibbaro, A. Lamura, G. Pontrelli, M. Sbragaglia, S. Succi, and F. Toschi, Eur. Phys. J. Spec. Top. 166, 111 (2009).

[32] A. Parmigiani, C. Huber, O. Bachmann, and B. Chopard, J. Fluid Mech. 686, 40 (2011).

[33] Z. Yu, O. Hemminger, and L. S. Fan, Chem. Eng. Sci. 62, 7172 (2007).

[34] L. Clime, D. Brassard, J. P. Pezacki, and T. Veres, Microfluid. Nanofluid. 12, 371 (2012).

[35] A. S. Joshi and Y. Sun, Phys. Rev. E 82, 041401 (2010).

[36] M. Sbragaglia, K. Sugiyama, and L. Biferale, J. Fluid Mech. 614, 471 (2008).

[37] M. Sega, M. Sbragaglia, S. S. Kantorovich, and A. O. Ivanov, Soft Matter 9, 10092 (2013).

[38] Q. Li, Z. Chai, B. Shi, and H. Liang, Phys. Rev. E 90, 043015 (2014).

[39] M. Sbragaglia, L. Biferale, G. Amati, S. Varagnolo, D. Fer-raro, G. Mistura, and M. Pierno, Phys. Rev. E 89, 012406 (2014).

[40] H. Farhat, S. Kondaraju, S. K. Na, and J. S. Lee, Phys. Rev. E 88, 013013 (2013).

[41] H. P. Jansen, K. Sotthewes, J. van Swigchem, H. J. W. Zandvliet, and E. S. Kooij, Phys. Rev. E 88, 013008 (2013).

[42] X. Shan and H. Chen, Phys. Rev. E 47, 1815 (1993).

[43] N. S. Martys and H. Chen, Phys. Rev. E 53, 743 (1996).

[44] A. Eddi, K. G. Winkels, and J. H. Snoeijer, Phys. Fluids 25, 013102 (2013).

[45] J. C. Bird, S. Mandre, and H. A. Stone, Phys. Rev. Lett. 100, 234501 (2008).

[46] K. G. Winkels, J. H. Weijs, A. Eddi, and J. H. Snoeijer, Phys. Rev. E 85, 055301 (2012).

[47] L. Courbin, J. C. Bird, M. Reyssat, and H. A. Stone, J. Phys.: Condens. Matter 21, 464127 (2009).

[48] L. H. Tanner, J. Phys. D: Appl. Phys. 12, 1473 (1979).

[49] A. L. Biance, C. Clanet, and D. Quéré, Phys. Rev. E 69, 016301 (2004).

[50] A. B. D. Cassie and S. Baxter, Trans. Faraday Soc. 40, 546 (1944).

[51] B. B. J. Stapelbroek, H. P. Jansen, E. S. Kooij, J. H. Snoeijer, and A. Eddi, Soft Matter 10, 2641 (2014).