Electrokinetic droplet transport from electroosmosis to electrophoresis

Abstract

Droplet transport in microfluidic channels by electrically induced flows often entails the simultaneous presence of electroosmosis and electrophoresis. Here we make use of coupled lattice-Boltzmann/molecular dynamics simulations to compute the mobility of a droplet in a microchannel under the effect of an external electric field. By varying the droplet solvation free energy of the counterions released at the channel walls, we observe the continuous transition between the electroosmotic and electrophoretic regime. We show that it is possible to describe the mobility of a droplet in a unified, consistent way, by combining the theoretical description of the electroosmotic flow with, in this case, the Hückel limit of electrophoresis, modified in order to take into account the Hadamard–Rybczynski droplet drag.

---

1 Introduction

The ability to manipulate and transport droplets in a controlled fashion is one of the central technological assets in modern microfluidics. Droplet-based microfluidic devices¹ are using small amounts of liquids, typically in the range from micro- to picoliters, in the form of a binary immiscible liquid mixture at low Reynolds numbers, forming a dispersed component (the droplets) carried in a continuum component (the medium). Dealing with such small amounts of liquids allows fast and controlled mixing thanks to the advantageous dimensional scaling. Individual control over the droplets makes them cheap and viable microreactors² that can be transported and analysed selectively³ or in parallel,⁴ allowing to achieve higher throughput than in continuous phase microfluidics.

Droplets can be moved around in microfluidics devices using a variety of techniques,⁵ among which electrokinetic approaches like electroosmosis,^6,7 electrophoresis,^8–10 dielectrophoresis,^11,12 or induced-charge electrokinetic effects¹³ have proven to be very popular. Electrokinetic phenomena, as opposed to electrohydrodynamics ones,¹⁴ are characterised by the presence of local charge imbalances in the fluid due to the presence of a counter-ions cloud (the so-called diffuse layer) that screens surface charges. A flow can be then generated by acting with an external electric field on the diffuse layer of the continuum phase or of the droplet, in the case of electroosmosis and in the case of electrophoresis respectively.

The distinction between electrokinetics and electrohydrodynamics, similarly to the one between electroosmosis and electrophoresis, has mostly a historical, rather than a physical justification. In practice, electrophoresis and electroosmosis often appear simultaneously. In the case of colloidal electrophoresis, for example, this complicates the correct evaluation of different contributions to the transport properties and, in turn, the measurement of the colloidal charge itself.^15,16 Droplet electrophoresis represents an even more complex problem as in principle counterions can permeate the droplet.

From the computational point of view, following Rotenberg and Pagonabarraga,¹⁷ one could classify different computational approaches based on whether the solvent and/or the ions are treated explicitly or implicitly. Complex flows in microfluidic channels have been often performed at the mesoscale^18–27 using numerical approaches like the lattice Boltzmann method,²⁸ kinetic theory approaches,^29,30 dissipative particle dynamics^31,32 or the multi-particle collision dynamics^33–36 methods, which allow, at the expense of microscopic detail, to overcome the size limitation of atomistic investigations, which are still limited to the nanoscale.^37–40 Regarding the electrophoresis of charged droplets in an electrolyte solution, the theoretical predictions of Ohshima and coworkers⁴¹ have been compared to numerical approaches, for example, using control volume techniques.⁴²

The problem of the electroosmotic contribution to the droplet motion in nanochannels, to the best of our knowledge, has not been addressed so far. In principle, a recently proposed approach for the solution of the Nernst–Plank equations using the lattice-Boltzmann method²⁷ could represent a viable and efficient method to investigate this issue using continuum ionic distributions.

In this study, we employ a molecular dynamics/lattice-Boltzmann coupling method to investigate the continuous transition between the purely electroosmotic transport regime and the electrophoretic-dominated one, for a droplet confined in a slit pore. We explore the transition between these two regimes by systematically varying the solvation free energy of counterions in the dispersed phase.

2 Methods

We use an implementation of the bicomponent lattice-Boltzmann method of Shan and Chen,⁴³ which solves the fluctuating hydrodynamics of the fluid, extending also the Guo forcing term⁴⁴ to the bicomponent case.⁴⁵ The lattice-Boltzmann equations are coupled to the molecular dynamics simulation of point-like ions.^45,46 All simulations are performed using an in-house modified version of the ESPResSo^47,48 simulation package that includes a CPU version of the bicomponent lattice Boltzmann code.

Our model consists of a cubic simulation box of size 48 × 48 × 48 lattice units Δx with two planar walls of surface area S = 48 × 48 Δx² each positioned so that the distance between the hydrodynamic no-slip planes is d = 42Δx. This grid resolution has been previously shown to be enough to reproduce the analytical results of the ion distribution and of the velocity field for the electroosmotic flow of a single-component fluid.

The fluid has two components, a and b, with total density ρ = ρ_a + ρ_b = 5.0/Δx³ and barycentric velocity u. We solve the fluctuating hydrodynamic equations of two uncharged fluids with same dielectric permittivity (dielectrophoretic forces are not taken into account)

where the forcing term f_ζ that governs the interaction between the fluid components is

In the equations above, the coupling parameter g_c = 0.8Δx³/Δt², where ζ and ζ′ indices for the fluid components, δ_ζζ′ is the Kronecker symbol, r identifies the position of a node and r′ loops over adjacent nodes, so that the term (r′ − r)ρ(r′) represents a discretised gradient of the density. In addition, p is the pressure, and Π and D_ζ are, respectively, the stress tensor and diffusion current, whose fluctuating parts obey the corresponding fluctuation-dissipation theorem.⁴⁵ In the present simulation scheme all non-conserved modes are relaxed independently and, in particular, the shear stress and interspecies diffusion relaxations are governed by the kinematic viscosity of the two fluid species, ν_a = ν_b = 40Δx²/Δt, and by the interspecies mobility M = Δt/6, respectively. The dispersed phase has the volume fraction ϕ ≃ 0.04 of the total fluid and the form of a spherical droplet with radius R = 11.2 Δx, calculated as the distance from the centre of the droplet to the shear surface. The thickness ξ of the interface is ξ ≃ 5Δx and is therefore not within the sharp interface limit.⁴⁹ The no-slip hydrodynamic boundary condition is imposed at the surface of the walls and is implemented as a bounce-back⁵⁰ with tuneable wettability feature.^51,52

Each pair (i,j) of charges interacts via the Coulomb potential U_c = k_BT[small script l]_Bq_iq_j/r with q_i being the valency of the i-th particle. The characteristic ratio of the energy of thermal fluctuations and the electrostatic energy between two particles is given by the Bjerrum length [small script l]_B = e²/(4πε₀ε_rk_BT), where ε₀ is the dielectric permittivity of vacuum, ε_r is the relative dielectric constant, e the elementary charge and k_BT is the thermal energy. We set l_B = Δx with no dielectric contrast between the fluid components. The two walls are located at z_B1 = 2Δx and z_B2 = 44Δx and are decorated with 576 immobile, pointlike ions each with q = −1, thus bearing a surface charge density σ_e = −0.125e/Δx², corresponding to the Gouy–Chapman length [small script l]_GC = e/(2πl_Bσ_e) ≃ 1.27Δx.

The same number of pointlike counterions with opposite valency q = 1 are free to access the slit pore volume and therein confined by a Weeks–Chandler–Anderson potential βU_WCA = 4ε[(σ/z)¹² − (σ/z)⁶ + 1/4] for z < σ2^1/6, and zero otherwise, which depends on the distance z of the ions to each of the walls. The wall-ions interaction (ε = 50k_BT, σ = Δx) prevents positive and negative ions to collapse on each other. The initial distribution of counterions is taken from an equilibrated single-component simulation.

The electrostatic coupling parameter Ξ = (eq)²l_b/l_GC = 2πl_B²σ_e ≃ 0.8 and the ratio of the wall distance to the average ion distance on the plate is about 15; therefore, the system is in the weak coupling regime.^53,54 Strictly speaking, since no salt is present in the system, the Debye screening length cannot be defined; instead, a screening constant κ from the Poisson–Boltzmann solution of the charged walls can be defined as , where ρ₀ is the charge density in the middle of the channel. Using the approximation ,⁵⁵ one can estimate κ ≃ 0.07/Δx. Hence, one can define a reduced screening length with respect to the droplet radius as δ = 1/(κR) ≃ 1.28. We apply an electric field of strength E, ranging from 0.05 to 1.0k_BT/(eΔx); This quantity should be compared to the potential drop equivalent of the thermal energy (k_BT ≃ 25 meV at room temperature) across the droplet, which defines a reduced electric field strength E* = eER/k_BT.

We apply periodic boundary conditions along the x and y directions (parallel to the slit pore), taking into account the long-range electrostatic interaction between periodic copies in the xy plane, using the electrostatic layer correction (ELC)^56,57 modification of the P3M algorithm.⁵⁸ The lattice-Boltzmann simulation is coupled to the molecular dynamics simulation according to the scheme of Ahlrichs and Dünweg⁵⁹ by integrating with timestep Δt = 0.01 a modified Langevin equation:

ma_i = F_i,ext − γ[v_i − u(r_i)] + F_i,R + F_i,ps, (3)

where γ = 10/Δt is a bare friction coefficient, F_i,R is a stochastic term with zero mean and variance 〈F_i,R(t)·F_j,R(t′)〉 = 6k_BTγδ_ijδ(t − t′). The momentum transferred from the fluid to the particles via the viscous coupling in eqn (3) is then redistributed back to the fluid to ensure momentum conservation using a trilinear interpolation.⁴⁶F_i,ext represents the external forces acting on the i-th ion, namely, the electric force eq_iE and the force from the repulsive walls. F_ps is the particle solvation force, which models the solvation free energy of counterions by acting on particles depending on the gradient of the dispersed phase density

F_ps = −ΔG∇ρ_a. (4)

We use ΔG to denote the work needed to bring one particle from the middle of the droplet component deep into the other phase. Far from the interface the density is basically constant, and there the ions do not feel any attraction/repulsion by the solvation forces.

3 Results

In Fig. 1 we report two simulation snapshots at different values of ΔG, showing the isodensity surface at half of the maximum density. Fixed ions on the channel surface and mobile counterions are depicted using blue and red spheres respectively. We apply the electric field along the y axis, parallel to the channel surface.

Fig. 1 Two snapshots of the system. The grey outline marks the surface of the droplet, while the orange and the blue spheres represent mobile and surface ions, respectively. Top: ΔG = 0.0k_BT; bottom: ΔG = 100k_BT.

The mobility μ of the droplet is defined by the terminal velocity v_t reached by the droplet in stationary conditions under the effect of an applied electric field E, as μ = v_t/E. Since what we want to address is the problem of linear electrokinetic transport, we first checked, in which range of the reduced field E* the terminal velocity depends linearly on E* itself. In Fig. 2 we report the reduced terminal velocity v_t*, which corresponds to the square root of the Weber number

where σ_s is the surface tension (for the present choice of parameters, σ_s ≃ 0.41k_BT/Δx², computed via the Laplace pressure jump), as a function of E* for different values of the solvation free energy ΔG. The terminal velocity is clearly linear at least up to values E* = 4. By using the result of the best fit to a linear function in the range E* = [0, 2.5], as shown in Fig. 2 using dashed lines, we report, in Fig. 3, the reduced mobilityas orange squares.

Fig. 2 Reduced droplet terminal velocity as a function of the reduced electric field for different values of the solvation free energy ΔG. Dashed lines are the result of a fit to linear functions in the E* ∈ [0, 2.5] interval.

Fig. 3 Reduced total droplet mobility μ* as a function of the solvation free energy ΔG of the counterions (orange squares); electroosmotic contribution, μ_eof*, calculated using eqn (9) (light blue line); electrophoretic mobility, μ_ep*, calculated from the droplet charge using eqn (8) (blue circles); electrophoretic mobility, μ_ep*, calculated by subtracting μ_eof* from the total mobility μ* (orange bars, no points).

When ΔG is negative, the droplet phase repels the ions and the mobility is slightly lower than in the neutral case, ΔG = 0. After we increase the solvation free energy to positive values, the mobility raises sharply, until the droplet transport enters into what seems to be a saturated regime. Intuitively, it seems straightforward to interpret the mobility values at ΔG ≤ 0 as mainly caused by the electroosmotic flow induced by the counterions. In fact, the droplet is in the middle of the channel, where the ion density is the lowest, and the electroosmotic flow is mainly generated in the high ion density regions next to the wall surface.

As soon as the droplet acquires some charge, when ΔG > 0, one could expect the onset of an electrophoretic behaviour. More precisely, for the droplet to move as a single charged object, the maximum solvation force needs to be larger than the electric force acting on each ion, or, ΔG/k_BT > E*ξ/R, which, in our case, is valid for ΔG/k_BT ≥ 1. In other words, given the parameters we have chosen, a free energy barrier large enough to prevent the thermal escape of the counterions will also prevent the electric field from stripping counterions from the droplet. Therefore, it makes sense to compute the charge of the droplet and to test whether its mobility at high solvation free energies is proportional to its charge Q, after removing the electroosmotic contribution.

If our droplet were a solid colloid, we would expect an electrophoretic behaviour of the Hückel type, since the ratio of the droplet radius to the screening length κR < 1. In this case, the mobility will be μ = Q/(6πRνρ), if the slip surface is located at R. However, the droplet is not solid and it is therefore not correct to use the Stokes friction formula in the derivation of the mobility. In the limit of low droplet deformations,⁶⁰ however, there is a simple solution for the friction coefficient D_visc = F/v_t of a droplet subject to a force F, as found by Hadamard and Rybczynski,^61,62 and later rederived by Booth,⁶³

where λ is the viscosity ratio of the inner/outer fluids. In our case λ = 1, hence, D_visc = 5πνρR, where one readily recognises a Stokes-like expression with a coefficient of 5 instead of 6 for the solid sphere, which leads to a Hückel electrophoretic mobility

μ_ep = Q/(5πηR) (8)

for a charge distributed homogeneously in the droplet (we neglect conductivity contributions). However, not all charges in the droplet will contribute to the electrophoretic mobility. As a rough estimate, we consider the number of ions within the droplet radius R in the absence of solvation force to be always freely moving. Using this value of the charge to calculate the electrophoretic mobility, eqn (8), we report its rescaled value in Fig. 3 as blue circles.

The behaviour of the electrophoretic mobility contribution calculated this way resembles closely that of the total mobility, apart from a roughly constant offset. It is reasonable to believe that the origin of this offset is the electroosmotic flow, which drags the droplet along. From this perspective, the fact that the droplet mobility at values ΔG < 0 is lower than at ΔG = 0 can be interpreted as the droplet behaving as a charge hole in the local charge background, since the repelling solvation force is depleting the droplet of ions. Therefore, in order to interpret the observed mobility, we need not only the expression for the mobility, eqn (8), but also an estimate for the electroosmotic flow contribution. Here our Ansatz is that we can approximate the electroosmotic contribution as that in the middle of an identical channel in absence of the droplet. We take into account finite size effects by considering only the charges that are not trapped in the droplet as a source of the electroosmotic flow. Then an effective surface charge can be defined as σ_e′ = σ_e − Q(ΔG)/S and used to compute the Gouy–Chapman length, and, in turn, the approximation for the electroosmotic contribution to the mobility⁵⁵

The electroosmotic contribution μ_eof is reported in reduced units as a continuous line in Fig. 3 and over the whole range of ΔG simulated, the contribution changed by 17% of its maximum value. With all contributions at hand, we are now ready to check that it is possible to express the droplet mobility as

μ ≃ μ_eof + μ_ep (10)

with the expressions for the electrophoretic and electroosmotic contributions given by eqn (8) and (9). In Fig. 3 we report the difference between the reduced mobility μ* and the reduced electrophoretic contribution μ_eof*, which is found to be in good agreement with the reduced electrophoretic contribution μ_ep*.

4 Conclusions

The problem of transport of droplets under microscopic confinement is relevant for a large number of microfluidic applications. The presence of counterions released in solution from the confining surfaces, by charging the fluids, opens the possibility of controlling droplets using electric fields. Even in the linear electrokinetic regime, where the potential drop across the droplet is smaller than the thermal energy, understanding the droplet transport is far from a trivial task, because of the superposition of the electroosmotic flow contribution and the electrophoretic mobility. The two contributions cannot be easily separated experimentally, a well-known issue in the context of colloidal electrophoresis, which has never been tackled in the case of fluid droplets. Computer simulations provide the possibility to disentangle these two phenomena by giving direct access to the droplet charge. Using a combination of on- and off-lattice simulation methods we modelled the transport of a droplet in a microfluidic channel as a function of the ions' solvation free energy difference between the dispersed phase and the medium. The solvation free energy of the ions is, in the present approach, an independent parameter. This allowed us to study the transition from the electroosmotic regime, where the counterions are dragging the medium, to the electrophoretic regime, in which the droplet moves as a charged object.

The present simulation does not take into account dielectric mismatch between the two fluids (nor with the confining medium). The electrophoretic mobility, however, as O'Brien and White pointed out,⁶⁴ should not depend on the dielectric mismatch between solvent and, in their case, the colloid, but only on the zeta-potential. In other words, dielectric boundary forces should not influence the mobility in the linear regime. The same applies to droplets or bubbles.⁹ In our case, a change in dielectric mismatch would alter the Bjerrum and the Gouy–Chapman length. In practical terms, one would need to take care of using the correct values for the Bjerrum length (which depends on the dielectric permittivity) and of the surface charge (which should include polarisation charges).

In order to interpret the simulation results, we formulated a model for the total droplet mobility, which combines the electrophoretic contribution in the Hückel limit and an analytical expression for the electroosmotic flow. It is worth noticing that the well-known expression for the electrophoretic mobility of droplets by Baygent and Saville,⁹ based on the solution proposed by Booth⁶³ reduces, in the zero salt concentration limit, to the Stokes case Q/(6πηR). This expression, however, applies only to uncharged droplets,⁶³ therefore, not to the present case, where the Hadamard–Rybczynski solution holds.

In summary, we have applied a mesoscopic simulation technique to the study of droplet transport in microfluidic channels. The simulation results support the interpretation of the total mobility as the superposition of an electroosmotic and an electrophoretic terms. The expression we proposed relates the total droplet mobility to its charge, as a function of known parameters such as fluid viscosity and channel surface charge density. This expression could be of practical relevance for the determination of individual droplet charge in microfluidic devices.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work has been supported by the Marie-Skłodowska-Curie European Training Network COLLDENSE (H2020-MSCA-ITN-2014 Grant No. 642774). The computational results presented have been achieved using the Vienna Scientific Cluster (VSC). This work was partially funded by the Austrian Research Fund (FWF): START-Projekt Y 627-N27.

References

1. S.-Y. Teh, R. Lin, L.-H. Hung and A. P. Lee, Lab Chip, 2008, 8, 198–220.
2. H. Song, J. D. Tice and R. F. Ismagilov, Angew. Chem., Int. Ed., 2003, 42, 768–772.
3. H. Song, D. L. Chen and R. F. Ismagilov, Angew. Chem., Int. Ed., 2006, 45, 7336–7356.
4. E. Z. Macosko, A. Basu, R. Satija, J. Nemesh, K. Shekhar, M. Goldman, I. Tirosh, A. R. Bialas, N. Kamitaki and E. M. Martersteck, et al. , Cell, 2015, 161, 1202–1214.
5. T. M. Squires and S. R. Quake, Rev. Mod. Phys., 2005, 77, 977.
6. J. N. Israelachvili, Intermolecular and Surface Forces, Academic Press, Amsterdam, 3rd edn, 2011.
7. P. N. Nge, C. I. Rogers and A. T. Woolley, Chem. Rev., 2013, 113, 2550–2583.
8. A. Tiselius, Trans. Faraday Soc., 1937, 33, 524–531.
9. J. Baygents and D. Saville, J. Chem. Soc., Faraday Trans., 1991, 87, 1883–1898.
10. S. A. Nespolo, M. A. Bevan, D. Y. Chan, F. Grieser and G. W. Stevens, Langmuir, 2001, 17, 7210–7218.
11. H. A. Pohl, J. Appl. Phys., 1951, 22, 869–871.
12. K. Ahn, C. Kerbage, T. P. Hunt, R. Westervelt, D. R. Link and D. A. Weitz, Appl. Phys. Lett., 2006, 88, 024104.
13. M. Z. Bazant and T. M. Squires, Phys. Rev. Lett., 2004, 92, 066101.
14. M. Z. Bazant, J. Fluid Mech., 2015, 782, 1–4.
15. I. Semenov, O. Otto, G. Stober, P. Papadopoulos, U. Keyser and F. Kremer, J. Colloid Interface Sci., 2009, 337, 260–264.
16. I. Semenov, S. Raafatnia, M. Sega, V. Lobaskin, C. Holm and F. Kremer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 022302.
17. B. Rotenberg and I. Pagonabarraga, Mol. Phys., 2013, 111, 827–842.
18. C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech., 2010, 42, 439–472.
19. J. Zhang, Microfluid. Nanofluid., 2011, 10, 1–28.
20. E. Chiarello, A. Gupta, G. Mistura, M. Sbragaglia and M. Pierno, Phys. Rev. Fluids, 2017, 2, 123602.
21. A. Gupta and M. Sbragaglia, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 023305.
22. J. Smiatek and F. Schmid, Advances in Microfluidics, InTech, London, 2012, pp. 97–126.
23. J. Smiatek and F. Schmid, Comput. Phys. Commun., 2011, 182, 1941–1944.
24. J. Harting, C. Kunert and J. Hyväluoma, Microfluid. Nanofluid., 2010, 8, 1.
25. S. Van der Graaf, T. Nisisako, C. Schroen, R. Van Der Sman and R. Boom, Langmuir, 2006, 22, 4144–4152.
26. J. Smiatek, M. Sega, C. Holm, U. D. Schiller and F. Schmid, J. Chem. Phys., 2009, 130, 244702.
27. N. Rivas, S. Frijters, I. Pagonabarraga and J. Harting, J. Chem. Phys., 2018, 148, 144101.
28. R. Benzi, S. Succi and M. Vergassola, Phys. Rep., 1992, 222, 145–197.
29. U. M. B. Marconi and S. Melchionna, Langmuir, 2012, 28, 13727–13740.
30. S. Melchionna and U. M. B. Marconi, Europhys. Lett., 2011, 95, 44002.
31. P. Hoogerbrugge and J. Koelman, Europhys. Lett., 1992, 19, 155.
32. P. Español and P. Warren, Europhys. Lett., 1995, 30, 191.
33. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605–8613.
34. A. Malevanets and R. Kapral, J. Chem. Phys., 2000, 112, 7260–7269.
35. M. Ripoll, K. Mussawisade, R. Winkler and G. Gompper, Europhys. Lett., 2004, 68, 106.
36. R. Kapral, Adv. Chem. Phys., 2008, 140, 89–146.
37. J. B. Freund, J. Chem. Phys., 2002, 116, 2194–2200.
38. V. Knecht, H. Risselada, A. Mark and S. Marrink, J. Colloid Interface Sci., 2008, 318, 477–486.
39. D. M. Huang, C. Sendner, D. Horinek, R. R. Netz and L. Bocquet, Phys. Rev. Lett., 2008, 101, 226101.
40. M. Sega, M. Sbragaglia, L. Biferale and S. Succi, Soft Matter, 2013, 9, 8526–8531.
41. H. Ohshima, T. W. Healy and L. R. White, J. Chem. Soc., Faraday Trans. 2, 1984, 80, 1643–1667.
42. S. Bhattacharyya and P. S. Majee, Phys. Fluids, 2018, 30, 082008.
43. X. Shan and H. Chen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 1815–1819.
44. Z. Guo, C. Zheng and B. Shi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 046308.
45. M. Sega, M. Sbragaglia, S. S. Kantorovich and A. O. Ivanov, Soft Matter, 2013, 9, 10092–10107.
46. P. Ahlrichs and B. Dünweg, J. Chem. Phys., 1999, 111, 8225–8239.
47. H.-J. Limbach, A. Arnold, B. A. Mann and C. Holm, Comput. Phys. Commun., 2006, 174, 704–727.
48. A. Arnold, O. Lenz, S. Kesselheim, R. Weeber, F. Fahrenberger, D. Roehm, P. Košovan and C. Holm, Meshfree methods for partial differential equations VI, Springer, 2013, pp. 1–23.
49. F. Magaletti, F. Picano, M. Chinappi, L. Marino and C. M. Casciola, J. Fluid Mech., 2013, 714, 95–126.
50. S. Succi, The lattice Boltzmann equation: for fluid dynamics and beyond, Oxford University Press, Oxford, 2001.
51. N. S. Martys and H. Chen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 53, 743.
52. H. Huang, D. T. Thorne Jr, M. G. Schaap and M. C. Sukop, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 066701.
53. M. Gouy, J. Phys. Theor. Appl., 1910, 9, 457–468.
54. D. L. Chapman, Lond. Edinb. Dubl. Phil. Mag., 1913, 25, 475–481.
55. A. Bazarenko and M. Sega, J. Mol. Liq., 2018, 271, 301–304.
56. A. Arnold, J. de Joannis and C. Holm, J. Chem. Phys., 2002, 117, 2496–2502.
57. J. de Joannis, A. Arnold and C. Holm, J. Chem. Phys., 2002, 117, 2503–2512.
58. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
59. P. Ahlrichs and B. Dünweg, J. Chem. Phys., 1999, 111, 8225–8239.
60. T. Taylor and A. Acrivos, J. Fluid Mech., 1964, 18, 466–476.
61. J. Hadamard, Acad. Sci., Paris, C. R., 1911, 152, 1735–1738.
62. W. Rybczynski, Bull. Acad. Sci. Cracovie Ser. A, 1911, 1, 40–46.
63. F. Booth, J. Chem. Phys., 1951, 19, 1331–1336.
64. R. W. O'Brien and L. R. White, J. Chem. Soc., Faraday Trans. 2, 1978, 74, 1607–1626.

---