Lattice-Boltzmann Simulations of Colloidal Particles
at Fluid Interfaces
Jens Harting
Department of Applied Physics, Eindhoven University of Technology P.O. Box 513, NL-5600MB Eindhoven, The Netherlands
and
Faculty of Science and Technology, Mesa+ Institute, University of Twente P.O. Box 217, NL-7500AE Enschede, The Netherlands
E-mail: j.harting@tue.nl
Particle-stabilized fluid interfaces are very common in industrial applications as they can be found for example in the food, cosmetics, or oil industries. However, until recently our un-derstanding of these systems was mostly based on experiments, highly simplified theoretical calculations and only very few numerical approaches. The lack of well established and widely applied simulation codes is not surprising since computer simulations of particle stabilized fluid interfaces are a highy complex task. Suitable algorithms need to be able to treat the hydrody-namics of the involved solvents, the dyhydrody-namics of the suspended colloidal particles, as well as the interactions between all those constituents at the same time. Furthermore, the simulation of large scale 3D emulsions requires highly efficient and massively parallel implementations and access to state of the art supercomputers. In these lecture notes we summarize the relevant details of our own simulation method for particle stabilized interfaces which is based on a com-bined Lattice-Boltzmann and molecular dynamics solver. We provide an overview on important implementation details and review a number of recent applications from our group with the aim to demonstrate the specific features of particle stabilized fluid interfaces.
1
Introduction
Particle stabilized emulsions have a high potential for various purposes with industrial ap-plications, such as cosmetics, improved low-fat food products, ice cream, drug delivery, or tertiary oil recovery1–3. While amphiphilic surfactant molecules are traditionally
em-ployed as emulsification agents, their effects can be mimicked or supplemented by the use of colloidal particles. These may be a cheaper or less toxic alternative to surfactants, but most importantly they may be customized to include additional desirable properties. Examples include ferromagnetic particles4, 5, particles with different interfacial properties
on different parts of their surface (for example Janus particles)6–10, or nonspherical
parti-cles11–18, where the geometric anisotropy has an impact on their stabilization properties. Those chemically or geometrically anisotropic particles will under certain conditions de-form the surrounding fluid-fluid interface. This leads to capillary interactions which can be used to tune the attraction or repulsion of adsorbed colloidal particles and as such might find applications such as the formation of new soft and highly tunable materials17, 10.
The energy differences involved in the adsorption of colloidal particles at a fluid in-terface are generally orders of magnitude larger than thermal fluctuations. Therefore, this adsorption process is practically irreversible19, 18 and Ostwald ripening can be fully
blocked20–24. In this manner, particles allow for long-term stabilization of an emulsion15, 25.
However, colloidal particles stabilize fluid interfaces kinetically and not thermodynam-ically. They reduce the interfacial free energy by their presence: maintaining a fluid-fluid
interface is energetically expensive, and it is favourable to replace it with particle-fluid in-terface. This is exactly what happens when a particle adsorps to such a fluid-fluid inin-terface. In contrast, due to their amphiphilic ineractions with the involved fluids, surfactants reduce the interfacial tension directly, which also reduces the total interfacial free energy. These effects are described in some detail in Frijters at al.26. Note that emulsions are in general not necessarily thermodynamically stable, as the energy cost of retaining the remaining in-terfacial area is still larger than the gain in entropy when the fluids mix completely. There are several ways to improve the stabilization properties, for example by using anisotropic particles with an additional rotational degree of freedom15.
Particle-stabilized emulsions can present in various forms, classified by the shapes and sizes of their fluid domains. The form that most resembles a traditional, surfactant-stabilized emulsion is the “Pickering emulsion”27, 28, which consists of particle-covered
droplets of one fluid suspended in another fluid. A more recent discovery is the bicon-tinuous interfacially jammed emulsion gel, or “bijel”, which is characterized by two large continuous fluid domains that are intertwined and are only stable because of the particles present at their interfaces. This form was first predicted by numerical simulations29, and
shortly thereafter confirmed experimentally30, 31. Parameters that affect the final state of an emulsion include the ratio between the two fluid components, the volume fraction of the particles, and their wettability. These parameters have been studied numerically by various authors15, 25, 32–35.
Computer simulations of particle stabilized emulsions require a solver for the dynam-ics of the involved fluid species combined with an algorithm to track the suspended parti-cles and their interactions. Here we use the Lattice-Boltzmann (LB) method for the fluid components, coupled to solid particles whose inter-particle interactions are simulated by molecular dynamics. This approach has proven very successful during recent years due to its ease of implementing multiphase flows, inherent parallelism allowing to harness the power of state of the art supercomputers, and straightforward methods to couple the fluid-and particle solvers. The simulation method is briefly explained in the following section 2. Within the scope of this lecture we then review three distinct applications of the method where we studied fundamental properties of particle stabilized fluid interfaces before we conclude in Sec. 4. The current lecture nodes represent a shortened collection of several original articles36, 26, 37, 38, 18, 15, 25and the interested reader is referred to those for more de-tails.
2
Simulation Method
2.1 The Lattice-Boltzmann method
For the simulation of the fluids we apply the Lattice-Boltzmann method which is based on the discrete form of the Boltzmann equation:39
fic(x + ci∆t, t + ∆t) = fic(x, t) + Ωci(x, t), (1)
wherefc
i(x, t) is the single-particle distribution function for a fluid component c with
discrete lattice velocity ciat timet located at lattice position x. Numerous discretizations
D3Q19 lattice. The distance between lattice nodes is given by the lattice constant∆x and nineteen directions are allowed for the velocity.∆t is the timestep and
Ωci(x, t) =− fc i(x, t)− f eq i (ρc(x, t), uc(x, t)) (τc/∆t) (2)
is the Bhatnagar-Gross-Krook (BGK) collision operator40. The fluid density is defined as
ρc(x, t) = ρ0
X
i
fic(x, t), (3)
whereρ0is a unit mass factor.τcdenotes the relaxation time for componentc and
fieq(ρc, uc) = ζiρc 1+ci·u c c2 s +(ci·u c)2 2c4 s − (uc ·uc) 2c2 s +(ci·u c)3 6c6 s − (uc ·uc) (c i·uc) 2c4 s (4) is a discretized third order expansion of the Maxwell-Boltzmann distribution function.
cs= 1 √ 3 ∆x ∆t (5)
is the speed of sound,
uc=X
i
fic(x, t)ci/ρc(x, t) (6)
is the fluid velocity andζidenotes a coefficient depending on the direction:ζ0= 1/3 for
the zero velocity,ζ1,...,6= 1/18 for the six nearest neighbors and ζ7,...,18= 1/36 for the
next nearest neighbors in diagonal direction. The kinematic viscosity can be calculated as νc= c2 s∆t τc ∆t− 1 2 . (7)
In the following we choose∆x = ∆t = ρ0= 1 for simplicity, but a conversion to SI units
is trivial41. In all simulations the relaxation time is set toτc
≡ 1. 2.2 A multicomponent Lattice-Boltzmann method
Several extensions for the Lattice-Boltzmann method have been developed to simulate multicomponent and multiphase fluids42–46. The article Liu et al.47 provides an extensive overview. Here, we use the multicomponent pseudopotential method introduced by Shan and Chen42. Every species has its own set of distribution functions following Eq. (1). A mean field force
Fc(x, t) =−Ψc(x, t)X c0 gcc0 X x0 Ψc0(x0, t)(x0− x) (8)
is calculated locally and acts between the components. The summation includes the dif-ferent fluid speciesc0 and x0, the nearest neighbors of lattice positions x. g
cc0 is a phe-nomenological coupling constant between the species andΨc(x, t) is a monotonous weight
function representing an effective mass. For the results presented here, the form
is used. To incorporate Fc(x, t) in fieqwe define ∆uc(x, t) =τcFc(x, t)
ρc(x, t) . (10)
The macroscopic velocity included infieqis shifted by∆ucas
uc(x, t) = P
ific(x, t)ci
ρc(x, t) − ∆u
c(x, t). (11)
As we are interested in immiscible fluids we choose a positive value forgcc0 which leads to a repulsive interaction and thus to the emergence of surface tension. A typical range of the coupling parameter is0.08≤ gcc0 ≤ 0.14.
2.3 Colloidal particles
We couple a traditional molecular dynamics algorithm to the Lattice-Boltzmann solver in order to compute the interactions and trajectories of the suspended particles. The particle trajectories follow Newton’s equations of motion
F= m ˙upar, D= J ˙ωpar, (12)
which are integrated using a classical leap frog integrator. F and D are the force and torque acting on the particle with massm and moment of inertia J. uparandωparare the velocity
and the rotation vector of the particle.
The particles are discretized on the Lattice-Boltzmann lattice. They are coupled to both fluid species by a modified bounce-back boundary condition as pioneered by Ladd and Aidun33, 48–52. The lattice Boltzmann equation then becomes
fic(x + ci, t + 1) = f¯ic(x + ci, t) + Ω¯ic(x + ci, t) +C. (13)
C depends linearly on the local particle velocity, ¯i is defined in a way that ci = −c¯i
is fulfilled. A change of the fluid momentum due to a particle leads to a change of the particle momentum in order to keep the total momentum conserved:
F(x, t) = 2f¯ic(x + ci, t) +Cc¯i. (14)
If the particle moves, some lattice nodes become free and others become occupied. The fluid on the newly occupied nodes is deleted and its momentum is transferred to the particle as
F(x, t) =−X
c
ρc(x, t)uc(x, t). (15)
A newly freed node (located at x) is filled with the average density of theNFNneighboring
fluid lattice nodes xiFNfor each componentc
33, ρc(x, t) ≡ N1 FN X iFN ρc(x + c iFN, t). (16)
Hydrodynamics leads to a lubrication force between the particles. This force is reproduced automatically by the simulation for sufficiently large particle separations. If the distance between the particles so small that there is no lattice node between them anymore, this is
supposed to fail. Therefore, if the smallest distance between two identical spheres with radiusR is smaller than a critical value ∆c= 23a lubrication correction is introduced as52
Fij = 3πµR2 2 ˆrij(ˆrij(ui− uj)) 1 rij− 2R− 1 ∆c . (17)
Here,µ is the dynamic viscosity, ˆrij a unit vector pointing from one particle center to
the other one and ui is the velocity of particlei. These corrections can be generalized
to nonspherical particles in several ways. A simple, but less accurate approach is based on a scaling of the potentials following Berne and Pechukas53, 54, 35. By appropriately de-termining the distance between particle surfaces and also including tangential forces, the accuracy can be improved substantially55. For very small particle distances, the diverging nature of the lubrication correction can cause the simulation to become unstable. Since this case only happens a very few times even in very dense systems, we introduce a Hertz potential56which has the following shape for two identical spheres with radiusR:
φH= KH(2rp− r)5/2forr < 2rp. (18)
Here,r is the distance between particle centers. For larger distances φH vanishes.KH is
a force constant and is chosen to beKH = 100 in all our simulations.
The Shan-Chen forces also act between a node in the outer shell of a particle and its neighboring node outside of the particle. This leads to an unphysical increase of the fluid density around the particle. In order to avoid this effect, the lattice nodes in the outer shell of the particle are filled with a virtual fluid with a fluid density corresponding to the average of the value in the neighboring free nodes for each fluid component:
ρcvirt(x, t) = ρc(x, t). (19)
This can be used to control the wettability of the particle surface. We define the parameter ∆ρ named “particle color”. Positive values of ∆ρ are added to the “red” fluid component,
ρrvirt= ρr+ ∆ρ, (20)
and negative values are added to the “blue” component,
ρbvirt= ρb+|∆ρ|. (21)
This leads to an approximately linear relation between∆ρ and the three-phase contact angleθp35.
2.4 Implementation
The locality of the Lattice-Boltzmann algorithm makes the implementation of massively parallel codes straightforward57. Generally, a regular, orthogonal grid is used, and the collision operator and boundary implementations are based on local operations so that at each lattice node only information from its own location is required. The computationally most demanding parts of a typical Lattice-Boltzmann code are “streaming” and “collision” routines. In practice, the number of floating point operations required at every timestep is not the limiting factor, but the random access in memory required for the streaming and eventual non-local force computations makes the algorithm a memory-bound numerical method. This implies a number of issues for the implementation of highly efficient codes,
1 2 4 8 16 32 64 128 256 1 2 4 8 16 32 64 128 256
speedup normalized to 1024 cores
number of cores [1024] LB3D v6.5.4 ideal (a) 1 2 4 8 16 32 64 128 256 1 2 4 8 16 32 64 128 256
speedup normalized to 1024 cores
number of cores [1024] LB3D v6.5.4 ideal
(b)
Figure 1: Strong scaling of LB3D on JUQUEEN. (a) relates to a system with only one fluid component, (b) to one with two fluid species and suspended particles37, 35.
be it on shared-memory multicore nodes, distributed memory clusters, or accelerator cards such as graphical processing units (GPUs) or Xeon Phi cards.
Several well established and highly scalable Lattice-Boltzmann implementations exist, which have all demonstrated excellent scaling on hundreds of thousand CPUs. To name a few, Ludwig58, LB3D59, walBerla60, MUPHY61, and Taxila LBM62are all codes which
have been used for numerous high quality scientific publications. Most of these codes are able to handle solid objects suspended in fluids. The first three can even combine this with multiple fluid components or phases.
Our own implementation, LB3D, is written in Fortran 90 and parallelized using MPI. Long-running simulations on massively parallel architectures require parallel I/O strategies and checkpoint and restart facilities. LB3D uses the parallel HDF5 library for I/O which has proven to be highly robust and performant on all supercomputers we had access to. LB3D has been shown to scale almost linearly on up to 262,144 cores on the European Blue Gene systems Jugene and Juqueen based at the J¨ulich Supercomputing Centre in Germany36, 47.
However, such excellent scaling required some optimizations of the code: Initially, LB3D showed only low efficiency in strong scaling beyond65 536, which could be related to a mismatch of the network topology of the domain decomposition in the code and the network actually employed for point-to-point communication. The Blue Gene/P provides direct links only between direct neighbors in a three-dimensional torus, so a mismatch can cause severe performance losses. Allowing MPI to reorder process ranks and manually choose a domain decomposition based on the known hardware topology, efficiency can be brought close to ideal. The importance of these implementation details is depicted by strong scaling measurements based on a system of1 0242×2 048 lattice sites carrying only
one fluid species (Fig. 1a) and a similarly sized system containing two fluid species and 4 112 895 uniformly distributed particles with a diameter of ten lattice units (Fig. 1b).
Wall
Wall
Fluid 1
Fluid 2
F < 0
F > 0
R
h
θ
Figure 2: Equilibrium state of a spherical particle at a fluid-fluid interface. The contact angle,θ = cos−1(h/R), whereh is the height of the particle centre of mass above the interface and R the particle radius, is defined with respect to Fluid 1. The force,F , on the particle acts to detach it into the wetting (positive force) or non-wetting (negative force) fluid (Figure reprinted from Davies et al.18).
3
Applications
3.1 Detachment energies of spheroidal particles from fluid-fluid interfaces
The detachment energy of a single particle from a fluid-fluid interface plays a crucial role in our understanding of particle-stabilised emulsions. Several authors have studied the detachment of particles in the past, but to the best of our knowledge, none of these extended their treatment to the case of anisotropic spheroidal particles, which we focused on in a recent paper by Davies et al.18. This section gives a summary of the results presented in
that publication.
Previous detachment energy studies focussing on free energy differences between an equilibrated particle at an interface63–70 and in the bulk revealed a crucial dependence on
particle shape: prolate and oblate spheroidal particles attach to interfaces more strongly because they reduce the interface area more than spherical particles for a given particle volume.71–75 For a particle already adsorbed at an interface to detach itself, the particle
must deform the interface and overcome the interface’s resistive force: there is a free-energy barrier and an associated activation free-energy.
We develop a simple thermodynamic model for the detachment energy of spheroidal particles from fluid-fluid interfaces as a function of contact angle and aspect ratio only, and highlight the implications of our simplifications. Assuming that the effect of gravity can be neglected, the surface free energy of a particle at an interface (Fig. 2) is given by
E = σ12A12+ σp1Ap1+ σp2Ap2, (22)
where Aij is the area and σij is the surface-energy of the i, j interface where
i, j ={1: fluid 1, 2: fluid 2, p: particle}.76 We neglect line-tension since it is relevant
only for nano-sized particles.72 The surface area of the particle isA
p = Ap1+ Ap2. The
free energy of a system in which the particle is fully immersed in either fluid1 or fluid 2 is given by Ei = σ12A12+ σpiApi wherei = 1, 2. Taking the free energy difference
between a spherical particle at an interface (Fig. 2) and a spherical particle immersed in the bulk fluid yields the detachment energy73, 74
Equilibrium Resistance Detachment
Figure 3: Snapshots of a prolate spheroidal particle of aspect ratioα = 2 under the influence of an external force detaching from an interface. For each snapshot, we run a simulation with the particle fixed and measure the resistive force from the interface on the particle (Figures reprinted from Davies et al.18).
For neutrally wetting micron-sized particles at an interface with surface-tensionσ12=
50 mN m−1, the detachment energy is much larger than the thermal energy,E/k BT ∼
107, and particles irreversibly attach to the interface. For nano-sized particles at the same
interface with very large or small contact angles, E ∼ kBT , and particles may freely
adsorb to and desorb from the interface19.
For our simulations, we initialise a system volume of size1283 in lattice units which
is half filled with liquid1 and half filled with liquid 2 (ρ(1) = ρ(2) = 0.7), such that an
interface forms atx = 64. The top and bottom are closed with solid walls. The particle density is chosen asρp = 2 which is an arbitrary choice. The particle is placed at the
interface and is not under the influence of any external forces. We first equilibrate the system until the interface diffuses and the particle establishes its equilibrium position, and hence its contact angle, on the interface.33 Then, we apply a constant external force to the particle. As stated above, the wettability of the particles can be tuned and we determine the contact angle by subtracting the height of the particle centre of mass above the interface (we linearly interpolate the interface position) and dividing by the particle radius,cos θ =
h R.
33, 26
To obtain the minimum detachment force, we employ a binary search algorithm: we start the algorithm by using the fact that the particle remains attached at the interface for a zero-force,Fatt = 0, and guessing a force which detaches the particle, Fdet. We then run
a new simulation with a force,Fnew = 12(Fatt+ Fdet) and repeat this procedure until we
determineFdetto the desired accuracy. As a next step, we run a single simulation with the
minimum detachment force, saving the simulation state frequently. We then run several simulations from the saved simulation snapshots (Fig. 3) but now with the particles fixed so that drag, buoyancy and gravity forces can be neglected. We let the systems from each snapshot equilibrate, and we measure the resulting force on the particle which is exactly the resistive force of the interface. This allows us to build a force-distance curveF (x). We fitF (x) with a fourth-order polynomial, which allows us to capture the linear regime and the detachment break-off regime accurately, and integrate the fitted function numerically to obtain the detachment energy. The detachment distance is the minimum distance at which the resistive force exerted by the interface on the particle is zero. As discussed shortly, the resistive force decreases discontinuously to zero at the point of detachment.
The restoring force provided by the interface to the particle as a function of displace-ment from equilibrium is shown in Fig. 4a. The corresponding fourth-order polynomial
(a) (b)
Figure 4: a) The normalized resistive force for a spherical particle of radiusR = 10 is approximately linear for small displacements as predicted by de Gennes et al.77, 78 and O’Brien.79 The symbols are simulation data and the dashed lines represent fourth-order polynomial fits to that data. The fourth-order fits are integrated in order to obtain the detachment energy,E
b) Dependence of the detachment energy on the aspect ratio,α, for several different contact angles, θ. Each set of coloured data points represents different wettabilities of the particles. Stars are theoretical calculations from our thermodynamic model in eq. (24) and eq. (25) and symbols are numerical data (Figures reprinted from Davies et al.18).
fits are also shown.
Davies et al.18 introduce a simplified thermodynamic model describing the detachment
energies of prolate and oblate spheroids in their equilibrium positions from interfaces.
∆E⊥= ¯h 2 2 1− α 2G(α) −¯h2 + α 4G(α), (24) ∆Ek= ¯h 2 2 1− α 2G(α) −¯h2 + 1 4G(α). (25)
These approximate expressions are simple quadratic functions of the dimensionless height ¯
h and the aspect ratio α. G(α) is a geometry factor. In Fig. 4b we compare the analytical results from Eq. (24) and Eq. (25) with our simulation data.
We find good agreement between our thermodynamic model and the numerical simula-tions. For neutrally wetting particles, the measured detachment energies from simulations are larger than those predicted by the thermodynamic model; this is expected since the particle has to deform the interface to overcome its resistive force. The differences are of the order of10%, suggesting that the thermodynamic model, which does not take into account interface deformations, is fairly accurate for the particle aspect ratios we investi-gated. Similarly, we find good agreement forθ = 68◦where the numerical data show a
higher detachment energy than predicted by the thermodynamic model, as expected. For contact anglesθ = 52◦ andθ = 32◦, we still find good qualitative agreement between
thermodynamic theory and numerical simulations for both prolate and oblate spheroids. However, for oblate spheroids the numerical detachment energy is less than the analytical predictions, though within errors of the order of the symbol size.
a)
b)
c)
Caeff = 0.04 Caeff = 0.08 Caeff = 0.12
Figure 5: Side-view of deformed droplets ( a)χ = 0.00, b) χ = 0.27 and c) χ = 0.55). The particles prefer to stay in the center of the channel where the shear flow is weakest. The particles further exhibit tank-treading-like behaviour: they move around the interface following the shear flow (Figure reprinted from Frijters et al.26).
3.2 Shear induced deformation of particle-armored droplets
Droplets covered by colloidal particles appear as a building block of emulsions, or might be produced as individual containers for example for drug delivery purposes. In this sec-tion we review some simulasec-tion results on how such droplets behave under an externally imposed shear flow. The section is based on an article by Frijters et al.26 and the reader is
referred to this article for more details.
A fractionχ of the interface of a droplet immersed in a second fluid is covered with particles and a shear flow is imposed by moving boundaries at the top and bottom of the simulation domain, which causes the droplet to deform. The resulting deformation can be quantified using Taylor’s dimensionless deformation parameter80, 81
D≡ L− B
L + B. (26)
HereL and B are the length and the breadth of the droplet. For small deformations, Taylor predicts a linear dependence of the deformation of a droplet on the capillary number (see Eq. 27), which becomesD = 35/32Ca for equiviscous fluids. Due to the finite system size, we define the capillary and Reynolds numbers in terms of an effective shear rate ˙γeff,
which is measured at the interface of the droplet during the simulation: Caeff ≡ µm˙γ effR d σ , Re eff ≡ ρm˙γ effR2 d µm , (27)
These are defined in terms os the dynamic viscosityµmof the medium, the radius of the
initial undeformed dropletRdand the surface tensionσ. The effective Reynolds number
0.04 0.08 0.12 Caeff 0.0 0.1 0.2 0.3 0.4 Dd 0.000 Caeff 0.12 8 16 24 R e eff (a) 0 40 80 120 160 ρp 0.2 0.3 0.4 Dd (b)
Figure 6: a) The deformation parameterD as a function of the effective capillary number Caeff forχ = 0 (red squares),χ = 0.27 (blue circles), χ = 0.41 (green diamonds), and χ = 0.55 (purple triangles). The effect of adsorped particles is very weak for lowχ, becomes noticeable at χ > 0.4. Taylor’s law is reproduced for smallCa (dashed line). The inset depicts that the Reynolds number scales linearly with the capillary number. b) The deformation parameterD as a function of the rescaled mass of the particles m∗
p, withχ = 0.55 and Caeff = 0.1. The inertia of the heavier particles causes additional deformation as they drag the droplet interface in the direction of the shear flow (Figures reprinted from Frijters et al.26).
Let us first understand the effect of nanoparticles on the deformation properties of the droplet, by discussing how they position themselves at and move over the droplet interface. The fluid-fluid interaction strength is held fixed atgbr= 0.10 and the particles with radius
rp = 5.0 are neutrally wetting (θp = 90◦) and have a mass density which is identical to
the fluid mass density. The simulation volume is chosen to benx= ny= 256, nz = 512,
with an initial droplet radius ofRinit
d = 0.3· nx = 76.8, while the number of particles
is varied asnp = 0, 128, 256, 320, 384, 446 and 512. This results in surface coverage
fractions betweenχ = 0 and χ = 0.55. The capillary number is changed by changing the shear rate. Some examples of the deformations thus realised are shown in Fig. 5, for Caeff = 0.04, 0.08, 0.12 and χ = 0.0 (a), χ = 0.27 (b) and χ = 0.55 (c).
The interface of a sheared and thus deformed droplet is increased as compared to its original spherical shape and more space is available for the particles to move freely over the interface (cf. Fig. 5). The particles are swept over the interface with increasing velocity as they move away from the center plane of the system and up the shear gradient. If the particles would not be affected by the shear flow, they would prefer to occupy interface with high local curvature as can be explained by a geometrical argument: the interface removed by a spherical particle at a curved interface is larger than the circular area removed from a flat interface. This explains why in this dynamic equilibrium, most particles can be found at the tips of the droplet (see Fig. 5 b)). Even though the overall structure of the particles on the droplet interface remains stable over time, individual particles move over the interface, performing a quasi-periodic motion.
The linear dependence of the deformation parameterD on the effective capillary num-ber introduced above is recovered for low capillary numnum-ber and low particle coverage (cf. Fig. 6a). When the coverage fraction grows beyondχ > 0.40 the deformations in this regime increase with increasingχ and constant capillary number. Combining Eqns. 27
results in a relation between the capillary and Reynolds number: Reeff = σ ρ mRd µ2 m Caeff (28)
Since we changeCaeff explicitly by changing the shear rate, Reeff is proportional to the capillary number for a fixed value of the surface tension and inertial effects increase the deformation. Furthermore, since the particles do not affect the surface tension, all curves of the Reynolds number versus capillary number have the same slope (cf. inset of Fig. 6a and Eq. (28)). This implies that the increased deformation in the case of added particles is not caused by changes in inertia of the fluids, but instead the inertia of the particles themselves plays a decisive role here: As shown in Fig. 6b, we have varied the mass of the particles over two orders of magnitude. We keptχ = 0.55 and Caeff = 0.1 constant and rescaled the mass scale with the reference mass: m∗p = mp/524. The particles are accelerated
as long as they are on the part of the droplet interface that experiences a shear flow at least partially parallel to the particle movement. Eventually, particles have to “round the corner”. The increased inertia of heavier particles makes it more difficult to change their movement, leading to a situation where the droplet interface is in fact initially dragged farther away in the direction of the shear flow instead. This then explains the increase of deformation with increasing particle mass. As our deformation is increased substantially, the system size limits the deformation we can induce. Therefore, the values presented here are underpredictions of the actual effect of increased mass at high deformations.
3.3 Timescales of emulsion formation caused by anisotropic particles
In this section different types of particle stabilized emulsions and the effect of the particle shape on some of their properties are discussed. The results presented here are a summary of an article recently published by G¨unther et al.35.
We find two different types of emulsions in our simulations, namely the Pickering emulsion (Fig. 7, left) and the bijel (Fig. 7, right). The choice of parameters (such as particle contact angle, particle concentration, fluid-fluid ratio, particle aspect ratio) deter-mines the type of emulsions. Parameter studies for emulsions have been discussed in the past.33, 35, 25 for spherical and ellipsoidal particles, respectively. Here, we limit ourselves to
anisotropy effects on the time dependence of the emulsion formation. We use prolate ellip-soids (m = 2; Fig. 7, top), spheres (m = 1; Fig. 7, center) and oblate ellipellip-soids (m = 1/2; Fig. 7, bottom). The interaction parameter between the two fluids (see Eq. (8)) is chosen asgbr = 0.08 which corresponds to a fluid-fluid interfacial tension of 0.0138. The
parti-cles are neutrally wetting (contact angleθp = 90◦) and the particle volume concentration
is chosen asΞ = 0.24. The simulated systems have periodic boundary conditions in all three directions and a side length ofLS = 256∆x. Initially, the particles are distributed
randomly. At each lattice node a random value for each fluid component is chosen so that the designed fluid-fluid ratio is kept (1:1 for the bijels and 5:2 for the Pickering emul-sions). When the simulation evolves in time, the fluids separate and droplets/domains with a majority of one of the fluids form.
The average size of droplets/domainsL(t) can be determined by measuring L(t) = 1 3 X i=x,y,z L(t)i, with L(t)i=p 2π hk2 i(t)i . (29)
Figure 7: Snapshots of typical simulated Pickering emulsions (left) and bijels (right) after 105timesteps. The emulsions are stabilized by prolate ellipsoids (m = 2, top), spheres (m = 1, center) and oblate ellipsoids (m = 1/2, bottom). The parameter determining if one obtains a bijel or a Pickering emulsion is the fluid ratio which is chosen as 1:1 for the bijels and 5:2 for Pickering emulsions (Figure reprinted from G¨unther et al.15).
hk2
i(t)i = Pkk 2
i(t)ς(k, t)/Pkk 2
i(t) is the second-order moment of the
three-dimensional structure functionς(k, t) = (1/ςn)|ϕ0k(t)|. ϕ0= ˜ϕ− h ˜ϕi is the fluctuation of
˜
ϕ which is the Fourier transform of the order parameter field ϕ = ρr
− ρb.
The time development ofL(t) for the three different particle types (prolate, spherical and oblate (m = 2, 1 and 1/2) and for Pickering emulsions and bijels is shown in Fig. 8a. We can identify three regimes: in the first few hundred timesteps the initial formation of the droplets/domains starts. Then, the growth of droplets/domains is being driven by Ostwald ripening. At even later times, droplets/domains grow due to coalescence. When two droplets unify, the area coverage fraction of the particles at the interface is increased because the surface area of the new droplet is smaller than that of the two smaller droplets before. At some point the area coverage fraction of the particles is sufficiently high to prevent further coalescence. The state which is reached at that time is (at least kinetically) stabilized and one obtains a stable emulsion. The values forL(t) are larger for bijels than for Pickering emulsions. This can be explained by the way we calculateL(t) (see Eq. (29) and related text) using a Fourier transformation of the order parameter field.
It can clearly be seen that anisotropic particles are more efficient in interface sta-bilization than spheres since they can cover larger interfacial areas leading to smaller
0.5 1.0 1.5 2.0 2.5 3.0 3.5 0 2 4 6
L [units of parallel particle radius]
t [105 timesteps] bijel, m=2 bijel, m=1 bijel, m=1/2 Pickering, m=2 Pickering, m=1 Pickering, m=1/2 (a) 0 5 10 15 20 25 30 35 40
L [units of parallel particle radius]
t [105 timesteps] 2.85 2.90 2.95 3.00 m=2 m=1/2 (b)
Figure 8: a) Pickering emulsion and bijel: Time development of the average domain sizeL(t) (see Eq. (29)) for m = 1 and m = 2. At first view, a steady state is reached after about 105timesteps. L(t) is larger for bijels than for Pickering emulsions, which is due to the measurement being based on the Fourier transform of the order parameter. Ellipsoids are able to stabilize larger interface areas than spheres leading to smallerL(t).
b) Bijel: Zoom form = 2 and m = 1/2: Time development of the average domain sizes depicting the impact of the additional timescales. The range of the variation ofL(t) is larger as compared to the Pickering emulsion due to the impact of a small deformation on the larger effective interface of the bijel (Figures reprinted from G¨unther et al.15).
fluid domains (note that the simulation volume is kept constant). However, the differ-ence inL(t) for m = 2 and m = 1/2 is small. This can be understood as follows: if a neutrally wetting prolate ellipsoid is adsorbed at a flat interface, it occupies an area AP,F(m > 1) = m1/3Ap,s, whereAp,s is the occupied interface area for a sphere with
the same volume. This corresponds in the case ofm = 2 to the occupied interface being larger by a factor of1.26 as compared to spheres. For an oblate ellipsoid the occupied interface area isAP,F(m < 1) = m−2/3Ap,s which form = 1/2 is by a factor of 1.59
larger than the area occupied by spheres. Since in emulsions the interfaces are generally not flat, these formulae can only provide a qualitative explanation of the behavior ofL(t): If the interface curvature is not neglectable anymore, we lose some of the efficiency of interface stabilization, which is more pronounced form < 1. This explains why the value ofL(t) for m = 1/2 is only slightly smaller than for m = 2.
It seems thatL(t) reaches a steady state after some 105 timesteps for both types of
emulsions and for all three values ofm. However, if one zooms in one can observe that L(t) develops for a longer time period if the particles have a non-spherical shape. As will be demonstrated below, the reason for this phenomenon is the additional rotational degrees of freedom due to the particle anisotropy. Furthermore, the time development ofL(t) for emulsions stabilized by prolate particles requires more time than that for the oblate ones. If a particle changes its orientation as compared to the interface or a neighboring particle this generally changes the interface shape. In this way the domain sizes are influenced, leading to changes ofL(t) – an effect which is not observed for m = 1.
Fig. 8b depicts a zoom-in of the time development ofL(t) for bijels with m = 2 and m = 1/2, respectively. One observes that L(t) decays in both cases. The range of the decay is larger form = 2 than for m = 1/2. The time of reordering is much shorter form = 1/2 as compared to m = 2. These effects can be explained by the presence of
additional rotational degrees of freedom for the anisotropic particles. While oblate particles have only a single additional rotational degree of freedom as compared to spheres, prolate particles show an even more complex behavior due to their second additional rotational degree of freedom.
Th first additional timescale we find can be related to the adsorption of single particles and their further rotation towards the interface during the initial stages of emulsion forma-tions. However, while the decrease ofL(t) in Fig. 8b is certainly of the order of 106LB
timesteps, we find that the single particle dynamics can only be responsible for effects on scales of the order of103timesteps.
By means of time dependent orientational order parameters and pair correlation func-tions we investigate the impact of collective effects of many particles being adsorbed at a flat fluid-fluid interface. For prolate particles, here also the mutual orientation of the particles is important and one has an additional degree of freedom leading to particle ori-entational ordering. The presence of many particles at an interface leads to two additional timescales in the reordering. The first one is the rotation of the particle towards the in-terface. The particle rotates towards its final orientation parallel to the inin-terface. For higher particle concentrations the time needed for a particle to come to its final orienta-tion increases. Hydrodynamic as well as excluded volume effects become more important. Above a critical concentration not every particle reaches its “final” orientation. In sum-mary, time scales involved in many-particle effects are found to be of the order of 104to even several105timesteps.
In emulsions the interfaces are generally not flat. Pickering emulsions usually have (approximately) spherical droplets and a bijel has an even more complicated structure of curved interface. The simplest realization of a curved interface is a single droplet. By simulating the ordering of many prolate ellipsoids at a droplet surface, it was found that the ordering is somewhat faster than in the case of flat interfaces. This is not surprising because capillary interactions between the particles play a substantial role in the case of curved interface.
We can understand one of the additional timescales with the behavior of the ellipsoidal particles at a single droplet. The particles reorder and it can be shown that this leads to a small change of the shape of the droplet which is (almost) exactly spherical in the beginning32. A change of the interface shape caused by reordering of anisotropic particles leads to a change ofL(t). The reordering of particle ensembles at flat as well as spherical interfaces takes of the order of105 timesteps. This reordering takes place in idealized
systems with constant interfaces which do not change their shape considerably. In real emulsions, however, the interface geometry changes substantially during their formation. For example, two droplets of a Pickering emulsion can coalesce. After this unification the particle ordering starts anew. This explains the fact that the additional timescale we find in our emulsions is of the order of several106timesteps.
This reordering is pronounced in the case of curved interfaces, where the movement of the particles leads to interface deformations and capillary interactions. During the for-mation of an emulsion, droplets might coalesce (Pickering emulsions) or domains might merge (bijels). After such an event the particles at the interface have to rearrange in order to adhere to the new interface structure. Due to this, the local reordering is practically being “restarted” leading to an overall increase of the interfacial area on a timescale of at least several106timesteps.
Our findings provide relevant insight in the dynamics of emulsion formation which is generally difficult to investigate experimentally due to the required high temporal res-olution of the measurement method and limited optical transparency of the experimental system. It is well known that in general particle-stabilized emulsions are not thermody-namically stable and therefore the involved fluids will always phase separate – even if this might take several months. Anisotropic particles, however, provide properties which might allow the generation of emulsions that are stable on substantially longer timescales. This is due to the continuous reordering of the particles at liquid interfaces which leads to an increase in interfacial area and as such counteracts the thermodynamically driven reduction of interface area.
4
Conclusion
We introduced our simulation method for fluid interfaces stabilized by colloidal particles, gave some details of the implementation of a highly efficient and massively parallel simu-lation code and reviewed a number of recent applications.
When an interface is stabilized by colloidal particles, it is of fundamental interest to understand how reversible the adsorption of the particles is. For typical particle diameters of the orders of micrometres, the energy required to detach a particle can be of the order of several thousandkBT . For nanoparticles, however, thermal energies might be of the
same order as the gain in free energy due to the particles being adsorped. In Sec. 3.1 and the recent publication by Davies et al.18, we studied this problem in detail and in
particular presented a simplified thermodynamic model to estimate detachment energies of spheroidal particles. The model shows good agreement with our simulations and allows the experimentalist to do ad-hoc estimates in order to for example better understand how stable an emulsion stabilized by colloidal particles of different shapes will be.
Particles and surfactants show some important differences when used for the stabi-lization of fluid interfaces. As shown in Sec. 3.2 and a recent paper by Frijters et al.26,
particles do not directly change the interfacial tension, but reduce the interfacial free en-ergy by removing interface area. Furthermore, particles are massive objects and under extreme circumstances such as high mass density and strong shear forces, the inertia of the particles plays a role. This can, for example, lead to an enhanced deformation of particle armored droplets in shear flow.
The geometrical shape of the particles used to stabilize an emulsion plays an important role for the dynamics and long-term stability of the systems. This was demonstrated in Sec. 3.3 and the recent article by G¨unther at al.15, where we studied the appearance of additional timescales in the dynamics of the formation of emulsions stabilized by oblate or prolate ellipsoids. The dynamics during the adsorption process of a single particle as well as interfacial jamming or capillary interaction between particles play a crucial role. These effects cause a continuous, long-term reordering of already adsorped particles which might be utilized to counteract the inherent thermodynamic instability of these systems by improved kinetic stabilization.
Acknowledgments
We thank all the students and collaborators who have contributed to the results reviewed in this lecture and the corresponding original journal articles: Fernando Bresme, Peter V. Coveney, Gary B. Davies, Stefan Frijters, Florian G¨unther, Florian Janoscheck, Fabian Jansen, Badr Kaoui, Timm Kr¨uger, and Qingguang Xie.
References
1. E. Dickinson, Food emulsions and foams: Stabilization by particles, Current Opinion in Colloid & Interface Science, 15, 40, 2010.
2. D.A. Pink, M. Shajahan, and G. Razul, Computer simulation techniques for food science and engineering: simulating atomic scale and coarse-grained models, Food Structure, 1, 71, 2014.
3. D.E. Tambe and M.M. Sharma, The effect of colloidal particles on fluid-fluid interfa-cial properties and emulsion stability, Advances in Colloid and Interface Science, 52, 1, 1994.
4. S. Melle, M. Lask, and G.G. Fuller, Pickering emulsions with controllable stability, Langmuir, 21, 2158, 2005.
5. E. Kim, K. Stratford, and M.E. Cates, Bijels containing magnetic particles: a simula-tion study, Langmuir, 26, 7928, 2010.
6. B.P. Binks and P.D.I. Fletcher, Particles adsorped at the oil-water interface: A the-oretical comparison between spheres of uniform wettability and “Janus” particles, Langmuir, 17, 4708, 2001.
7. N. Glaser, D.J. Adams, A. B¨oker, and G. Krausch, Janus particles at liquid-liquid interfaces, Langmuir, 22, 5227, 2006.
8. J. Aveyard, Can Janus particles give thermodynamically stable Pickering emulsions?, Soft Matter, 8, 5233, 2012.
9. A. Kumar, B.J. Park, F. Tu, and D. Lee, Amphiphilic Janus particles at fluid interfaces, Soft Matter, 9, 6604, 2013.
10. Q. Xie, G. B. Davies, F. G¨unther, and J. Harting, Tunable dipolar capillary deforma-tions for magnetic Janus particles at fluid-fluid interfaces, Submitted for publication, 2015.
11. M.J.A. Hore and M. Laradji, Prospects of nanorods as an emulsifying agent of immis-cible blends, Journal of Computational Physics, 128, 054901, 2008.
12. W. Zhou, J. Cao, W. Liu, and S. Stoyanov, How rigid rods self-assemble at curved surfaces, Angewandte Chemie - International Edition, 48, 378, 2009.
13. M. Grzelczak, J. Vermant, E.M. Furst, and L.M. Liz-Marz´an, Directed self-assembly of nanoparticles, ACS Nano, 4, 3591, 2010.
14. L. Botto, E.P. Lewandowski, M. Cavallaro, and K.J. Stebe, Capillary interactions between anisotropic particles, Soft Matter, 8, 9957, 2012.
15. F. G¨unther, S. Frijters, and J. Harting, Timescales of emulsion formation caused by anisotropic particles, Soft Matter, 10, 4977, 2014.
16. G.B. Davies, T. Kr¨uger, P.V. Coveney, J. Harting, and F. Bresme, Interface deforma-tions affect the orientation transition of magnetic ellipsoidal particles adsorbed at fluid-fluid interfaces, Soft Matter, 10, 6742, 2014.
17. G.B. Davies, T. Kr¨uger, P.V. Coveney, J. Harting, and F. Bresme, Assembling ellip-soidal particles at fluid interfaces using switchable dipolar capillary interactions, Advanced Materials, 26, 6715, 2014.
18. G.B. Davies, T. Kr¨uger, P.V. Coveney, and J. Harting, Detachment energies of spheroidal particles from fluid-fluid Interfaces, Journal of Chemical Physics, 141, 154902, 2014.
19. B.P. Binks, Particles as surfactants – similarities and differences, Current Opinion in Colloid & Interface Science, 7, 21, 2002.
20. W. Ostwald, Analytische Chemie, Engelmann, 1901.
21. R.D. Vengrenovitch, On the Ostwald ripening theory, Acta Metallurgica, 30, 1079, 1982.
22. Y. Enomoto, K. Kawasaki, and M. Tokuyama, Computer modelling of Ostwald ripen-ing, Acta Metallurgica, 35, 907, 1987.
23. P.W. Voorhees, Ostwald ripening of two-phase mixtures, Annual Review of Materials Science, 22, 197, 1992.
24. P. Taylor, Ostwald ripening in emulsions, Adv. Colloid Interface Sci., 75, 107, 1998. 25. S. Frijters, F. G¨unther, and J. Harting, Domain and droplet sizes in emulsions
stabi-lized by colloidal particles, Physcal Review E, 90, 042307, 2014.
26. S. Frijters, F. G¨unther, and J. Harting, Effects of nanoparticles and surfactant on droplets in shear flow, Soft Matter, 8, 6542, 2012.
27. W. Ramsden, Separation of solids in the surface-layers of solutions and “suspen-sions”, Proceedings of the Royal Society of London, 72, 156, 1903.
28. S.U. Pickering, Emulsions, Journal of the Chemical Society, Transactions, 91, 2001, 1907.
29. K. Stratford, R. Adhikari, I. Pagonabarraga, J.-C. Desplat, and M.E. Cates, Colloidal jamming at interfaces: a route to fluid-bicontinuous gels, Science, 309, 2198, 2005. 30. E.M. Herzig, K.A. White, A.B. Schofield, W.C.K. Poon, and P.S. Clegg, Bicontinuous
emulsions stabilized solely by colloidal particles, Nature Materials, 6, 966, 2007. 31. P.S. Clegg, E.M. Herzig, A.B. Schofield, S.U. Egelhaaf, T.S. Horozov, B.P. Binks,
M.E. Cates, and W.C.K. Poon, Emulsification of partially miscible liquids using col-loidal particles: nonspherical and extended domain structures, Langmuir, 23, 5984, 2007.
32. E. Kim, K. Stratford, R. Adhikari, and M.E. Cates, Arrest of fluid demixing by nanoparticles: a computer simulation study, Langmuir, 24, 6549, 2008.
33. F. Jansen and J. Harting, From bijels to Pickering emulsions: a lattice Boltzmann study, Physcal Review E, 83, 046707, 2011.
34. S. Aland, J. Lowengrub, and A. Voigt, A continuum model of colloid-stabilized inter-faces, Physics of Fluids, 23, 062103, 2011.
35. F. G¨unther, F. Janoschek, S. Frijters, and J. Harting, Lattice Boltzmann simulations of anisotropic particles at liquid interfaces, Computers and Fluids, 80, 184, 2013. 36. D. Groen, O. Henrich, F. Janoschek, P.V. Coveney, and J. Harting, “Lattice-Boltzmann
methods in fluid dynamics: Turbulence and complex colloidal fluids”, in: J¨ulich Blue Gene/P Extreme Scaling Workshop 2011, Wolfgang Frings Bernd Mohr, (Ed.). J¨ulich Supercomputing Centre, 52425 J¨ulich, Germany, apr 2011, FZJ-JSC-IB-2011-02; http://www2.fz-juelich.de/jsc/docs/autoren2011/mohr1/.
37. J. Harting, S. Frijters, F. Janoschek, and F. G¨unther, Coupled lattice Boltzmann and molecular dynamics simulations on massively parallel computers, in: NIC Sympo-sium 2012 – Proceedings, M. Kremer K. Binder, G. M¨unster, (Ed.), vol. 45 of NIC Series, p. 243, Forschungszentrum J¨ulich GmbH Zentralbibliothek, Verlag. 2012, http://hdl.handle.net/2128/4538.
38. T. Kr¨uger, S. Frijters, F. G¨unther, B. Kaoui, and J. Harting, Numerical simulations of complex fluid-fluid interface dynamics, European Physical Journal Special Topics, 222, 177, 2013.
39. R. Benzi, S. Succi, and M. Vergassola, The lattice Boltzmann equation: theory and applications, Physics Reports, 222, 145, 1992.
40. P.L. Bhatnagar, E.P. Gross, and M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Physical Review, 94, 511, 1954.
41. A. Narv´aez, T. Zauner, F. Raischel, R. Hilfer, and J. Harting, Quantitative analysis of numerical estimates for the permeability of porous media from lattice-Boltzmann sim-ulations, Journal of Statistical Mechanics: Theory and Experiment, 2010, P211026, 2010.
42. X. Shan and H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Physical Review E, 47, 1815, 1993.
43. E. Orlandini, M.R. Swift, and J.M. Yeomans, A lattice Boltzmann model of binary-fluid mixtures, Europhysics Letters, 32, 463, 1995.
44. M. R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans, Lattice-Boltzmann simulations of liquid-gas and binary fluid systems, Physical Review E, 54, 5041, 1996. 45. S.V. Lishchuk, C.M. Care, and I. Halliday, Lattice Boltzmann algorithm for surface
tension with greatly reduced microcurrents, Physical Review E, 67, 036701, 2003. 46. T. Lee and P.F. Fischer, Eliminating parasitic currents in the lattice Boltzmann
equa-tion method for nonideal gases, Physical Review E, 74, 046709, 2006.
47. H. Liu, Q. Kang, C. R. Leonardi, B. Jones, S. Schmieschek, A. Narv´aez Salazar, J. R. Williams, A. J. Valocchi, and J. Harting, Multiphase lattice Boltzmann simulations for porous media applications - a review, Computational Geoscience, in press, 2014. 48. C.K. Aidun, Y. Lu, and E.-J. Ding, Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation, Journal of Fluid Mechanics, 373, 287, 1998.
49. C.K. Aidun and J.R. Clausen, Lattice-Boltzmann method for complex flows, Annual Review of Fluid Mechanics, 42, 439, 2010.
50. A.J.C. Ladd, Numerical simulations of particulate duspensions via a discretized Boltzmann equation. Part I. Theoretical foundation, Journal of Fluid Mechanics, 271, 285, 1994.
51. A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltz-mann equation. Part II. Numerical results, Journal of Fluid Mechanics, 271, 311, 1994.
52. A.J.C. Ladd and R. Verberg, Lattice-Boltzmann simulations of particle-fluid suspen-sions, Journal of Statistical Physics, 104, 1191, 2001.
53. B.J. Berne and P. Pechukas, Gaussian model potentials for molecular interactions, Journal of Chemical Physics, 56, 4213, 1972.
54. F. Janoschek, F. Toschi, and J. Harting, Simplified particulate model for coarse-grained hemodynamics simulations, Physcal Review E, 82, 056710, 2010.
55. F. Janoschek, J. Harting, and F. Toschi, Accurate lubrication corrections for spherical and non-spherical particles in discretized fluid simulations, Submitted for publica-tion, 2014.
56. H. Hertz, ¨Uber die Ber¨uhrung fester elastischer K¨orper, Journal f¨ur reine und ange-wandte Mathematik, 92, 156, 1881.
57. N. Satofuka and T. Nishioka, Parallelization of lattice Boltzmann method for incom-pressible flow computations, Computational Mechanics, 23, 164, 1999.
58. J.C. Desplat, I. Pagonabarraga, and P. Bladon, LUDWIG: A parallel lattice-Boltzmann code for complex fluids, Computer Physics Communications, 134, 273, 2001. 59. J. Harting, M. Harvey, J. Chin, M. Venturoli, and P.V. Coveney, Large-scale lattice
Boltzmann simulations of complex fluids: advances through the advent of computa-tional grids, Philosophical Transactions of the Royal Society of London A, 363, 1895, 2005.
60. C. Feichtinger, S. Donath, H. K¨ostler, J. G¨otz, and U. R¨ude, WaLBerla: HPC software design for computational engineering simulations, Journal of Computational Science, 2, 105, 2011.
61. M. Bernaschi, M. Fatica, S. Melchionna, S. Succi, and E. Kaxiras, A flexible high-performance Lattice Boltzmann GPU code for the simulations of fluid flows in com-plex geometries, Concurrency Computat.: Pract. Exper., 22, 1, 2010.
62. E.T. Coon, M.L. Porter, and Q. Kang, Taxila LBM: a parallel, modular lattice Boltz-mann framework for simulating pore-scale flow in porous media, Computational Geo-science, 18, 17, 2014.
63. A. Scheludko, B.V. Toshev, and D.T. Bojadjiev, Attachment of particles to a liquid surface (capillary theory of flotation), Journal of the Chemical Society Faraday Trans-actions I, 72, 2815, 1976.
64. A. W. Neumann and J. K. Spelt, Applied surface thermodynamics, Marcel Dekker, New York, 1996.
65. P.A. Kralchevsky and K. Nagayama, Particles at fluid interfaces and membranes, Elsevier Science, Amsterdam, 2001.
66. A. V Rapacchietta, A. W Neumann, and S. N Omenyi, Force and free-energy analyses of small particles at fluid interfaces: I. Cylinders, Journal of Colloid and Interface Science, 59, 541, 1977.
67. A. V Rapacchietta and A. W Neumann, Force and free-energy analyses of small par-ticles at fluid interfaces: II. Spheres, Journal of Colloid and Interface Science, 59, 555, 1977.
68. P. Singh and D. D. Joseph, Fluid dynamics of floating particles, Journal of Fluid Mechanics, 530, 31–80, 2005.
69. D. D. Joseph, J. Wang, R. Bai, B. H. Yang, and H. H. Hu, Particle motion in a liq-uid film rimming the inside of a partially filled rotating cylinder, Journal of Flliq-uid Mechanics, 496, 139–163, 2003.
70. I.B. Ivanov, P.A. Kralchevsky, and A.D. Nikolov, Film and line tension effects on the attachment of particles to an interface: I. Conditions for mechanical equilibrium of fluid and solid particles at a fluid interface, Journal of Colloid and Interface Science, 112, 97, 1986.
71. R. Aveyard and J.H. Clint, Particle wettability and line tension, Journal of the Chem-ical Society Faraday Transactions I, 92, 85, 1996.
72. J. Faraudo and F. Bresme, Stability of particles adsorbed at liquid/fluid interfaces: Shape effects induced by line tension, Journal of Chemical Physics, 118, 6518, 2003. 73. S. Levine, B.D. Bowen, and S.J. Partridge, Stabilization of emulsions by fine particles I. Partitioning of particles between continuous phase and oil/water interface, Colloids and Surfaces, 38, 325, 1989.
74. T. F. Tadros and B. Vincent, Encyclopedia of emulsion technology, Vol.1, p129, Mar-cel Dekker, New York, 1983.
75. J. Guzowski, M. Tasinkevych, and S. Dietrich, Capillary interactions in Pickering emulsions, Physical Review E, 84, 031401, 2011.
76. B. Binks and T. Horozov, Colloidal Particles at Liquid Interfaces, Cambridge Uni-versity Press, Cambridge, 2006.
77. P. G. de Gennes, Wetting: statics and dynamics, Reviews of Modern Physics, 57, 827, 1985.
78. J. F. Joanny and P. G. de Gennes, A model for contact angle hysteresis, J. Chem. Phys., 81, 552, 1984.
79. S. B. G. O’Brien, The meniscus near a small sphere and its relationship to line pinning of contact lines, Journal of Colloid and Interface Science, 183, 51, 1996.
80. G.I. Taylor, The viscosity of a fluid containing small drops of another fluid, Proceed-ings of the Royal Society A, 138, 41, 1932.
81. G.I. Taylor, The formation of emulsions in definable fields of flow, Proceedings of the Royal Society A, 146, 501, 1934.
