- Email: [email protected]
Numerical simulation of the dispersion of aggregated Brownian particles under shear flows
Computers & Fluids 86 (2013) 395–404
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Numerical simulation of the dispersion of aggregated Brownian particles under shear flows Takuya Nishiyama a,⇑, Takaji Inamuro b,c, Shugo Yasuda d a
Electronic and Imaging Materials Research Laboratories, Toray Industries Incorporated, Shiga 520-8558, Japan Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto 606-8501, Japan c Advanced Research Institute of Fluid Science and Engineering, Graduate School of Engineering, Kyoto University, Kyoto 606-8501, Japan d Graduate School of Simulation Studies, University of Hyogo, Kobe 650-0047, Japan b
a r t i c l e
i n f o
Article history: Received 24 October 2011 Received in revised form 27 May 2013 Accepted 24 June 2013 Available online 1 July 2013 Keywords: Particle dispersion Shear flow Van der Waals force Brownian motion Lattice Boltzmann method
a b s t r a c t The deformation and breakup processes of a particle-cluster aggregate under shear flows are numerically investigated by the two-phase lattice Boltzmann method. The van der Waals attraction is considered to be the force between particles. Simulations are performed for various fluid forces acting on particles and various inter-particle forces. It is found that the ratio of the fluid force to the maximum inter-particle force, Y, is a key factor in dispersion, and the aggregate of non-Brownian particles is dispersed when Y is over 0.001. The Péclet number, which is the ratio of the diffusion rate due to shear flow to that due to the Brownian motion, is also considered. By comparing the calculated result of the dispersion of Brownian particles with that of non-Brownian particles, it is found that the Brownian motion impedes dispersion and the effect of the Brownian motion is remarkable when the Péclet number is under 105. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The dispersion of micro- and nanoparticles in a liquid is important for making new functional materials such as ceramics, polymers, and electronic products because the characteristics of materials can be controlled by dispersing small particles in a solvent in their production process. However, small particles are easy to be aggregated by an attraction force between the particles, so it is difficult to disperse a large number of particles uniformly in a liquid. Therefore, to develop new functional materials, it is important to investigate particle dispersion in liquids. The breakup mechanism of aggregated particles has been studied both experimentally and theoretically [1–5]. Because the dispersion dynamics of particles is very complicated, it is difficult to investigate the breakup mechanism of aggregated particles only through experimental and theoretical approaches, so several numerical methods have been proposed. However, in general, there are difficulties in treating a moving solid–liquid boundary even in the numerical simulations of particulate flows. Bossis and Brady [6] proposed the Stokesian dynamics, in which a solvent fluid is not treated explicitly, but the hydrodynamic interaction between particles is treated as a resistance of the solvent fluid expressed by relative velocity vectors of the particles. Harada et al. [7] investigated the dispersion of particles by the ⇑ Corresponding author. E-mail address: [email protected] (T. Nishiyama). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.06.023
Stokesian dynamics. However, the Stokesian dynamics remains difficult for applications with complex boundaries. Immersed boundary approaches, in which a moving boundary condition is satisfied through external forces in the Navier–Stokes equation, have been proposed [8–10]. Feng and Michaelides [11], and Melchionna [12] also proposed an immersed boundary approach combined with the lattice Boltzmann method to solve a fluid–particle interaction problem. These models using the immersed boundary method are effective approaches for moving boundary problems, but the immersed boundary method might cause numerical errors while evaluating the hydrodynamic force and torque acting on particles in various situations. Tanaka and Araki [13] proposed the fluid particle dynamics method, in which a particle is treated as a fluid with large viscosity. Although this method can eliminate the difficulty originating from the solid-liquid boundary condition, it is difficult to retain the spherical shape for a long time. Ladd [14] proposed a simulation model for moving particles using the lattice Boltzmann method in which the hydrodynamic force exerted on a particle is calculated on the basis of momentum exchange. Recently, the two-phase lattice Boltzmann method with the same density [15] has been applied to the simulations of the dispersion of aggregated particles under shear flows by the authors [16]. In this method, the particle is modeled by a hard droplet with large viscosity and strong surface tension, and consequently, we do not need to explicitly track the moving solid–liquid boundary. Strong surface tension is applied to maintain a spherical droplet without any other artificial treatments. Some models of the lattice
396
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
Boltzmann method that involve droplets and tunable surface tension are also available [17–20]. Here the free energy functional scheme proposed by Swift et al. [17] is used to apply strong surface tension. The advantage of the two-phase lattice Boltzmann method is its simplicity in treating the solid-liquid boundary. Zaleski et al. [21], e.g., proposed the volume-of-fluid method, which involves droplets and tunable surface tension. However, the evaluation of the boundary gradient is complex. In the two-phase lattice Boltzmann method, the boundary surface is autonomously determined to minimize the free energy of the system. In addition, a colored order parameter is assigned to each droplet for preventing the droplets from merging into bigger droplets. In this study, we improve the above-mentioned method to treat the Brownian motion of nanoparticles and simulate a large number of particles with a smaller number of colored order parameters. Next, we use the improved method to investigate the behavior of particles under shear flows for various conditions. In general, particles diffuse by the Brownian motion in fluid, so, the Brownian motion may promote dispersion. However, it is also known that particles, which are initially separated, aggregate by the Brownian motion (see e.g., Ref. [16]). Therefore, the effect of the Brownian motion on the dispersion of initially aggregated particles is not yet clear. To classify the calculated results, we introduce two important dimensionless parameters: the ratio of the fluid force to the maximum inter-particle force and the Péclet number, which is the ratio of the diffusion rate due to shear flow to that due to the Brownian motion. The paper is organized as follows: In Section 2, we explain the numerical method. Simulation results are presented in Section 3, in which we investigate the dispersion of the aggregations of 36 and six particles. We conclude the study in Section 4. 2. Numerical method 2.1. Two-phase lattice Boltzmann method for immiscible fluids with the same density Non-dimensional variables are used as in Ref. [22]. The lattice kinetic scheme (LKS) [23], which is an extension method of the lattice Boltzmann methods, is used to formulate the method. In the LKS, macroscopic variables are calculated without particle velocity distribution functions, and thus the scheme can save computer memory, since there is no need to store the particle velocity distribution functions. In addition, in order to represent many hard droplets, which cannot merge into bigger droplets, we introduce colored order parameters to make different colored droplets. Differently colored droplets do not merge if they collide because the boundary surface of each colored droplet is independently and autonomously determined. Note that the color is physically meaningless and is used only for distinguishing each droplet from the others. In the present study, the Stokes flow is assumed because the diameter of a particle is very small (e.g., 1 lm). The 15-velocity model with particle velocities ci(i = 1, 2, . . . , 15) is used in this study. The velocity vectors of this model are given by ½c 1 ; c 2 ; c 3 ; c 4 ; c 5 ; c 6 ; c 7 ; c 8 ; c 9 ; c 10 ; c 11 ; c 12 ; c 13 ; c 14 ; c 15 2 3 0 1 0 0 1 0 0 1 1 1 1 1 1 1 1 6 7 ¼ 4 0 0 1 0 0 1 0 1 1 1 1 1 1 1 1 5: 0 0 0 1 0 0 1 1 1 1 1 1 1 1 1
ð1Þ
i¼1
uðx; t þ DtÞ ¼
ð2Þ
15 1X g eq ðx ci Dx; tÞ; 3 i¼1 i
ð3Þ
15 X c i g eq i ðx c i Dx; tÞ:
ð4Þ
i¼1
where flieq and g eq are equilibrium distribution functions, Dx is a i spacing of the cubic lattice, and Dt is a time step during which the particles travel the lattice spacing. The equilibrium distribution functions in Eqs. (2)–(4) are given by
h i jf flieq ¼ Hi /l þ F i p0 ð/l Þ jf /l r2 /l jr/l j2 þ 3Ei /l cia ua 6 þ Ei jf Gab ð/l Þcia cib ;
ð5Þ
@ub @ua cia cib g eq þ i ¼ Ei 3p þ 3c ia ua þ ADx @xa @xb þ Ei jg Gab ð/l Þcia cib þ 3Ei cia Dx
N X
F la Ul ;
ð6Þ
l¼1
where
E1 ¼ 2=9;
E2 ¼ E3 ¼ E4 ¼ ¼ E7 ¼ 1=9;
E8 ¼ E9 ¼ E10 ¼ ¼ E15 ¼ 1=72; H1 ¼ 1;
H2 ¼ H3 ¼ H4 ¼ ¼ H15 ¼ 0;
F 1 ¼ 7=3;
F i ¼ 3Ei ði ¼ 2; 3; 4; . . . ; 15Þ;
ð7Þ
and
Gab ð/Þ ¼
9 @/ @/ 3 @/ @/ dab ; 2 @xa @xb 2 @xc @xc
ð8Þ
with a, b, c = x, y, z (subscripts a, b, and c represent Cartesian coordinates, and the summation convention is used). In the above equations, dab is the Kronecker delta, jf is a constant parameter determining the width of the interface, jg is a constant parameter determining the strength of the surface tension, and A is a constant parameter related to fluid viscosity. In Eq. (4), p0(/) is given by
p0 ð/Þ ¼ /T /
1 a/2 ; 1 b/
ð9Þ
where a, b, and T/ are free parameters determining the maximum and minimum values of the order parameter, /max and /min. It is noted that flieq is the same as that for the model of Swift et al. [17], except that flieq in Eq. (5) has no second terms of ua because we assume the Stokes flow. The last term of Eq. (6) is the force Fl acting on the lth particle per unit mass:
Fl ¼
N X
F v lm þ F Br l ;
ð10Þ
m¼1
where Fv lm is the attractive force per unit mass from the mth particle to the lth particle, and FBr l is the Brownian force per unit mass acting on the lth particle. The detailed formulations of forces are given in the following subsections. In addition, Ul in Eq. (6) represents the region of a droplet; i.e., it is unity inside the droplet and zero outside the droplet:
The physical space is divided into a cubic lattice, and the colored order parameter /l(x, t) (l = 1, 2, . . . , N), where N is the number of colors, the pressure p(x, t) and the velocity u(x, t) of whole fluid at the lattice point x and at time t are computed as follows: 15 X /l ðx; t þ DtÞ ¼ flieq ðx c i Dx; tÞ;
pðx; t þ DtÞ ¼
Ul ¼
1; /l P /B ; 0; /l < /B :
ð11Þ
/B is the threshold value of the boundary and is defined as follows:
/B ¼
f/min þ /max ; fþ1
ð12Þ
where f is a weight parameter and is described in detail in Appendix A.
397
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
The following finite-difference approximations are used to calculate the derivatives in Eqs. (5), (6), and (8): 15 @w 1 X cia wðx þ ci DxÞ; @xa 10Dx i¼1 " # 15 X 1 2 r w wðx þ c i DxÞ 14wðxÞ : 5ðDxÞ2 i¼2
ð13Þ ð14Þ
2.2. Inter-particle force
Applying the asymptotic theory [24] to Eqs. (2)–(9), it is found that the macroscopic variables (pressure p and velocity u) satisfy the incompressible Navier–Stokes equation with relative errors of O[(Dx)2] [22]. The kinematic viscosity, m, of the fluid is given by
1 2 m ¼ A Dx: 6 9
ð15Þ
Thus, kinematic viscosity can be changed the by A. In the interface, kinematic viscosity is changed as follows:
m ¼ mc þ
mp mc 2
/0 /B p þ1 ; sin /max /min
ð16Þ
where mp is the kinematic viscosity of the droplet, mc is that of the solvent, and /0 ¼ maxl2f1;2;...;Ng /l at the lattice point. We set the viscosity ratio of the droplet to the surrounding liquid to be g = md/ mc = 10, and the droplet and the liquid have equal density. According to the Hadamard–Rybczynski drag law [25], the ratio of the drag force for a spherical fluid particle to that for a rigid spherical particle is 0.970. Thus, the Stokes drag law is recovered within 3% error. The interfacial tension, rl, for the lth colored phase is given by
rl ¼ jg
2 @/l dn; @n 1
Z
physically meaningful current is set to be sufficiently large than the spurious current by tuning the parameters. The parameters we choose are described in Sections 3.1 and 3.2. By using the above method, we can form many hard droplets with large viscosity and strong surface tension to model solid particles. The validity of the method is examined in Appendix A.
The van der Waals attractive force is applied as the inter-particle force. The attractive force per unit mass from the mth particle to the lth particle Fvlm is given by [26]
Fv
lm
¼
8 > Ah > > < pD3 q
D2 r lm 2 ðr 2lm D2 Þ
p
> > > :
2
2 lm þ rD3 r22rD 2 þ r lm
lm
for D þ LD 6 r lm 6 2D; 0;
lm
r lm r lm
; ð18Þ
otherwise;
where D is the diameter of particle, qp is the density of the particle, rlm is the vector from the center of the lth particle to that of the mth particle, rlm = jrlmj, Ah is the Hamaker constant, and LD corresponds to the Debye length, which gives the thickness from the surface of the particle. The attractive force reduces to zero, if the distance between particle surfaces is smaller than LD. In addition, to avoid allto-all computations, the inter-particle force is cut off when the distance between the centers of the two particles is greater than 2D, because the van der Waals force is very small when two particles are far apart. The correction of the lubrication force between particles is not considered in the present study because its effect is negligible in the dispersion process of aggregated particles.
1
ð17Þ
with n being the coordinate normal to the interface. It is noted that a small spurious current is generated around the droplet when strong surface tension is applied. However, in the simulation,
2.3. Brownian motion We apply the Brownian motion to particles using the Langevin approach. The advantages of the Langevin approach are its simplicity and computational efficiency.
φ1 φ1 φ1 φ1
(a)
(b) Fig. 1. Color reassignment. (a) The same color is assigned to particles when they are far apart from each other. (b) Exchange of order parameters in certain areas.
398
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
The Brownian force acting on the lth particle per unit mass at every time step is formularized as follows:
F Br la
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6pqc mc DkB T Xa; ¼ Dt pD3 qp 6
ð19Þ
where kB is the Boltzmann constant, qc is the density of the solvent, T is the temperature of the solvent, and Xa is the normal random number whose average is 0 and variance is 1. Note that Xa is independently chosen for a = x, y, and z. We also note that T is not the temperature of the solvent in the experimental system but the temperature for a single Brownian particle obeying the Langevin equation in an infinite domain. The temperature, T, is used as a reference value of the solvent temperature in the present simulation. The difference between the temperatures T and Tc estimated from a numerical simulation of particle diffusion is given in Appendix B. Fig. 3. Computational domain and six aggregated particles.
2.4. Color reassignment As described in Section 2.1, we assign a different colored order parameter to each particle to prevent the particles from merging into a larger particle. However, if the number of particles increases, many colors are necessary, which increases computational time. To save computational time, the same color is assigned to particles that are far apart from each other, as shown in Fig. 1(a). Next, if particles of the same color approach each other, one particle’s color is changed. To change the color of the particle /1 to /2, the order parameter /1 and /2 are exchanged in the region within a specified distance (e.g., 1.6D/2) from the center of mass of the particle, as shown in Fig. 1(b). The accuracy of this scheme is examined in Appendix C. 3. Results and discussion In our previous paper [16], four cases with two, six, 18, and 36 non-Brownian particles are calculated. In that study, the aggregates with 18 and 36 particles are loose structures. In the present study, an aggregate composed of particles with hexagonal closepacked structure is used. In the following subsection, we compare the present results to the previous results from the view point of the effect of the initial arrangement of particles. Aggregated particles composed of same particles with diameter D are placed in a liquid inside a rectangular domain. To evaluate the cluster size effect, two different aggregations are considered: the aggregation of 36 particles (Fig. 2) and the aggregation of six particles (Fig. 3). Initially, the particles are arranged in a hexagonal close-packed structure for both aggregations. In addition, for the
aggregation of six particles, two domain sizes are used to investigate the domain size effect. At t = 0, the top and bottom walls begin to move in the x-direction with velocities uw and uw, respectively. Periodic boundary conditions are used on the other sides of the domain. The solid–liquid boundary is determined by the threshold value (/max + /min)/2, i.e., f = 1 in Eq. (12), so the diameter is determined by D1/2 (the definition of D1/2 is given in Appendix A). Although D1/4 is more suitable, as observed in Appendix A, we use D1/2 throughout this study for simplicity. The other parameters are set at the same values as those mentioned in Ref. [16]. In order to classify the calculated results, we introduce the following non-dimensional parameters: (a) The ratio of the fluid force [/qcmc(uwD/Lz)D] to the inter-particle force [/Ah/D]:
I¼
qc mc uw D3 Ah Lz
ð20Þ
:
(b) The ratio of LD to the diameter of the particle:
s ¼ LD =D:
ð21Þ
(c) The ratio of the fluid force [/qcmc(uwD/Lz)D] to the maximum inter-particle force [/Ah/s2D for s 1]:
Y¼
qc mc uw D3 s2 Ah Lz
¼ Is2 :
ð22Þ
Note that the inter-particle force is maximum when rlm = D + LD in Eq. (18). Table 1 Computational conditions.
Fig. 2. Computational domain and 36 aggregated particles.
Case Case Case Case Case Case Case Case Case Case Case Case Case Case Case
A B C D E F G H I J K L M N O
I
s
Y
Pe
1000 100 10 10 10 1000 100 100 0.625 10 2.5 1000 100 10 0.625
0.005 0.005 0.005 0.016 0.05 0.016 0.05 0.016 0.2 0.01 0.01 0.005 0.016 0.05 0.2
0.025 0.0025 0.00025 0.0025 0.025 0.25 0.25 0.025 0.025 0.001 0.00025 0.025 0.025 0.025 0.025
1 1 1 1 1 1 1 1 1 1 1 350,000 35,000 3500 220
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
(d) The Péclet number, which is the ratio of the rate of diffusion due to shear flow[/2uwD/Lz] to the rate of diffusion due to the Brownian motion [/kBT/3pqcmcD2]:
Pe ¼
6pqc mc uw D3 6pAh I: ¼ kB TLz kB T
ð23Þ
The computational conditions are given in Table 1. All conditions are implemented for the aggregation of 36 particles, whereas
399
Cases A–K (non-Brownian cases) are implemented for the aggregation of six particles. In Cases A–K, the Brownian motion is not considered. In Cases A–C, the effect of the intensity of the shear flow is investigated. In Cases C–E, the effect of the value of LD is investigated. In Cases F–K, the intensity of the shear flow and LD are varied. In Cases L–O, the Brownian motion is applied to investigate the dispersion of nanoparticles. In reality, Case A, for example, corresponds to the following situation, in which the parameters with a tilde are dimensional parameters in the real system. The particle is made of titanium oxide with a
Fig. 4. Time evolution of the aggregate under shear flows for (a) I = 1000, s = 0.005, Y = 0.025, and Pe = 1 (Case A); and (b) I = 10, s = 0.005, Y = 0.00025, and Pe = 1 (Case C). t⁄ = 2uwt/Lz.
400
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
e h ¼ 8:0 1020 J, diameter D e ¼ 1 lm, and Hamaker constant of A e L D ¼ 5 nm. The solvent is water with a kinematic viscosity of m~c ¼ 1:0 106 m2 s1 . The top and bottom walls at a distance of 10 lm move in opposite directions with a velocity of ~ w ¼ 0:8 m s1 . Moreover, Case N corresponds to the following situau tion. The particle is titanium oxide with a Hamaker constant of e h ¼ 8:0 1020 J, diameter D e ¼ 100 nm, and e A L D ¼ 5 nm. The sol~c ¼ 1:0 106 m2 s1 vent is water with a kinematic viscosity of m and a temperature of 308 K. The top and bottom walls at a distance ~ w ¼ 0:8 m s1 . of 1 lm move in opposite directions with a velocity of u
3.1. Results for 36 particles The domain is divided into a 240 100 200 cubic lattice. The diameter of the particle is D = 20Dx. The 18–20 colored order parameters are used to save computation time. We chose a = 9/ 49, b = 2/21, and T/ = 0.55 in Eq. (9). The other parameters are fixed at jf = 0.01Dx2, jg = 0.01Dx2, and uw = 0.1. The calculated results of Case A, C, E, and N are shown in Figs. 4 and 5. In these figures, t⁄ = 2uwt/Lz is dimensionless time. In these calculations, color reassignments were implemented 72, 39, 73, and 64 times for Cases A,
Fig. 5. Time evolution of the aggregate under shear flows for (a) I = 10, s = 0.05, Y = 0.025, and Pe = 1 (Case E) and (b) I = 10, s = 0.05, Y = 0.025, and Pe = 3500 (Case N). t⁄ = 2uwt/Lz.
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
Fig. 6. Time evolution of the standard deviation, r, for the aggregation of 36 particles (I = 1000, Y = 0.025, and Pe = 1 for Case A; I = 10, Y = 0.00025, and Pe = 1 for Case C; I = 10, Y = 0.025, and Pe = 1 for Case E; and I = 10, Y = 0.025, and Pe = 3500 for Case N). t⁄ = 2uwt/Lz.
401
Fig. 9. Time evolution of the standard deviation, r, for the aggregation of six particles in domain A. r0 is the initial standard deviation, and t⁄ = 2uwt/Lz.
Fig. 10. Time evolution of the standard deviation, r, for the aggregation of six particles in domain B. r0 is the initial standard deviation, and t⁄ = 2uwt/Lz. , versus Y for the aggregation of 36 Fig. 7. Time-averaged standard deviation, r is the time-averaged value for 140 6 t⁄ 6 164, and r0 is the particles for Pe = 1. r initial standard deviation.
1.2
1
Case L / Case A Case M / Case H Case N / Case E Case O / Case I
0.8
0.6
0.4 102
103
104
105
106
Pe , to Fig. 8. Ratio of the time-averaged standard deviation of the Brownian case, r ðPe¼1Þ , at Y = 0.025 for the aggregation of 36 that of the non-Brownian case, r ⁄ and r ðPe¼1Þ are time-averaged values for 140 6 t 6 164. particles. r
C, E, and N, respectively. From the results of Figs. 4 and 5, we find that the aggregated particles are dispersed by shear flow for Case A, whereas they are not dispersed for Case C because the shear flow is
weak (the value of I in Case C is one-hundredth of that in Case A). However, the aggregated particles are dispersed for Case E although I has the same value as in Case C because the value of s in Case E is ten times larger than that in Case C. In Case N, the Brownian motion is considered, but other conditions are the same as those in Case E. Comparing the results of Case E with those of Case N, it is observed that the Brownian motion impedes particle dispersion. To measure the dispersion state of the particles, the standard deviation r of the lengths between the center of the aggregation and those of the particles is calculated during the deformation of the aggregates. The time evolutions of r for Cases A, C, E, and N are shown in Fig. 6. In this figure, r0 is the standard deviation for the initial structure. The results of Fig. 6 show that the standard deviations for Cases A and E are close each other, which indicates that the parameter Y, which is a combination of I and s, should be a key factor determining particle dispersion. I is the parameter relating two particles far apart from each other, and Y is the parameter relating two particles that are almost in contact. On the other hand, it is clearly found from Fig. 6 that the standard deviation of 36 particles for Case N increases more slowly than that for Case E, which means that the Brownian motion impedes particle dispersion and that the dispersion of nanoparticles is more difficult. , for Fig. 7 shows the time-averaged standard deviation, r 140 6 t⁄ 6 164 for all cases with 36 non-Brownian particles from Table 1. The results in Fig. 7 show that dispersion begins when Y
402
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
is over 0.001. In addition, our previous results [16] indicate that aggregates with a loose structure begin to disperse when Y is over 0.0025, which is of the same order as 0.001, although only four cases for s = 0.005 are available. Thus, we consider that for dispersion, the initial particle arrangement does not considerably change the threshold value of Y. Fig. 8 shows the relationship between the Péclet number and time-averaged standard deviation for 140 6 t⁄ 6 164 of Brownian particles normalized by that of non-Brownian particles at the same value of Y = 0.025. It is found that the dispersion-retarding effect of the Brownian motion is remarkable when the Péclet number is under 105. 3.2. Results for 6 particles In order to examine the cluster- and domain-size effects, nonBrownian cases (Cases A–K in Table 1) are simulated for the aggregation of 6 particles with two different domain sizes. One domain is divided into a 132 80 110 cubic lattice, and the other is divided into a 198 120 165 cubic lattice. We call the former domain ‘‘domain A’’ and the latter ‘‘domain B.’’ The diameter of the particle is D = 20Dx. We use the same values for the parameters a, b, T/, jf, and jg as those mentioned in Section 3.1 and uw = 0.055 and 0.0825 for domains A and B, respectively. The time evolution of the standard deviation, r, for domain A (B) is shown in Fig. 9 (10). We can observe that the particles disperse in all cases except Cases C and K in which Y is over 0.001 for both domains A and B. Thus, the threshold value of Y for dispersion is 0.001 for six particles, which is the same as that for 36 particles, and the cluster- and domain-size effects are not clearly observed. It is noted that the dispersion of particles is slower in domain B, which is larger than domain A, because the development of shear flows depends on the height Lz of the domain.
Fig. A.1. Computational domain for the measurement of Stokes drag on a particle.
D) is placed at the center of a cubic domain of size L, as shown in Fig. A.1. Because the interface of the droplet is approximately 3Dx in thickness, we define the diameter of the droplet by using the threshold value, /B, given in Eq. (12). We denote the diameter by D1/1+f, which corresponds to the weight parameter, f. The usual periodic boundary conditions are used on all the sides of the domain, and periodic boundary conditions with pressure difference Dp are used at x = 0 and x = L. Because of the flow induced by the pressure difference, the hard droplet moves in the x-direction. We apply a uniform body force to the internal region of the droplet to keep the droplet stationary. The body force is applied to the region / > (9/max + /min)/10. After a steady state is obtained, we calculate the total body force, F, and the average velocity, U, at x = 0. Computational conditions and calculated results are given in Table A.1, in which the deformation factor, k1/1+f, is also shown. The droplet is slightly deformed by the flow, and we measure the deformation factor defined by
4. Concluding remarks To numerically investigate the dispersion of aggregated nanoparticles, we applied the lattice Boltzmann method for the twophase fluid flow with uniform density for simulating the dispersion of aggregated particles under shear flow. The ratio of fluid force to the maximum inter-particle force, Y, is introduced, and it is found that Y is the key factor for dispersion and the aggregates of nonBrownian particles are dispersed when Y is over 0.001. By comparing the calculated results of the dispersion of Brownian particles to those of non-Brownian particles, it is found that the Brownian motion impedes dispersion and the effect of the Brownian motion is remarkable when the Péclet number is under 105.
k1=1þf
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 !2 u I u D DII1=1þf 1=1þf t ¼ 1 þ 1 D1=1þf D1=1þf
ðA:1Þ
where DI1=1þf and DII1=1þf are diameters in the x- and the z-directions, respectively. It is found from Table A.1 that k1/1+f is greater than 5% in the smallest domain (Case I), but it is relatively small for the other cases. Fig. A.2 shows the calculated results of the relationship between the dimensionless drag force, Q, and the filling factor, Ct, which are defined as follows:
Appendix A. Validation of present method: drag force acting on a sphere
Q¼
3plD1=1þf U F
ðA:2Þ
and
We calculated a Stokes flow through the periodic arrays of spheres and compared the drag force acting on the sphere with the results of Zick and Homsy [27]. A hard droplet with an initial diameter D (initially, we assign /max inside D, and /min outside
Ct ¼
pD31=1þf =6 L3
ðA:3Þ
:
Table A.1 Computational conditions and calculated results for the measurement of Stokes drag. L Case Case Case Case Case Case Case
I II III IV V VI VII
31 35 41 45 60 70 90
Dp 3
1.3 10 7.5 104 4.0 104 2.5 104 7.0 104 3.5 104 1.5 104
U
F
D1/10
D1/4
D1/2
k1/10 (%)
k1/4 (%)
k1/2 (%)
0.067 0.066 0.066 0.069 0.070 0.067 0.066
1.25 0.919 0.671 0.607 2.52 1.71 1.19
22.6 22.5 22.4 23.4 43.2 43.4 44.7
21.2 21.0 21.0 20.8 41.8 42.3 43.4
19.9 19.7 19.6 19.5 40.5 40.9 42.1
5.6 2.6 1.7 1.3 2.5 1.3 0.32
5.4 3.1 1.9 1.2 3.2 0.92 0.66
5.3 2.4 1.1 1.1 2.6 1.0 0.34
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
403
Fig. B.1. Time evolution of mean-square displacement hjxBrj2i; t is the time step.
Fig. A.2. Comparison between the results of Zick and Homsy and those of the present method for various /B. Q is the dimensionless drag force, and Ct is the filling factor.
It is found from Fig. A.2 that the results estimated using D1/4 agree well with those of Zick and Homsy regardless of the droplet size. Appendix B. Validation of the present method: free diffusion of particle We calculated the Brownian motion of a single spherical particle in a three-dimensional cubic domain and compared the estimated temperature, Tc, with the given temperature, T, from the fluctuation–dissipation relation. A hard droplet with diameter D = 21Dx (D1/4) is placed at the center of a cubic domain of size
L = 160Dx. The usual periodic boundary conditions are used on all sides of the domain. The same random force per unit mass as Eq. (19) (except for Q) is applied to the particle at each time step:
F Br la
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6pqc mc DkB T Xa: ¼ Q Dt pD3 6
ðB:1Þ
We use kBT = 4.4 103, qcmc = 0.0167, and Q = 0.815. Here Q is defined by Eq. (A.2) as the ratio of the drag force acting on the particle in the periodic domain to that in the infinite domain. In addition, the additional Q in Eq. (B.1) is introduced in the same manner as that in Ref. [28]. However, we assume Q = 1 in the simulation of the dispersion of particles, as mentioned in Eq. (19). Fig. B.1 shows the time evolution of the mean-square displacement. It is found that the mean-square displacement is proportional to the time step, t, when t is over 3 104. The calculated
Fig. C.1. Time evolution of the aggregation of six particles under shear flows for Case A using (a) six colors and (b) seven colors with color reassignment; t⁄ = 2uwt/Lz. (c) Time evolution of total mass (enlarged view). (d) Time evolution of total mass (overall view).
404
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
diffusion coefficient is Dc0 ¼ 1:2 103 . The temperature of the system, Tc, estimated from the Stokes–Einstein relation is c kT ¼ 3pqc mc DDc0 =Q ¼ 4:9 103 . Therefore, the fluctuation–dissipation relation is almost satisfied; i.e. the estimated temperature, Tc, is 12% larger than the given temperature, T. Appendix C. Validation of the present method: color scheme We calculated the dispersion of six particles by using seven colors with the color reassignment scheme and compared the results with those obtained when the color reassignment scheme is not used. The domain is the same as domain A (132 80 110 cubic lattice), and the computational condition is the same as Case A from Section 3.2. Initially the particle color is the same as the particle number. The color of the sixth particle alternates between the sixth or seventh color for every 100 steps. The colors of other particles do not change. When the color of the sixth particle changes, the order parameters /6 and /7 are exchanged in the region within 16Dx from its center of mass. Figs. C.1(a) and (b) show the time evolution of the aggregate when the color reassignment scheme is not used and when seven colors are used with the color reassignment scheme, respectively. In the figure, t is the time step, and t⁄ = 2uwt/Lz is the dimensionless time. The two results are almost the same. In addition, Fig. C.1(c) and (d) show the time evolution of the total mass of the order parameters /6 and /7 with seven colors and the color reassignment scheme, respectively. It is evident that the total mass is conserved. Note that Fig. C.1(c) shows that the total mass varies over 100 steps in correspondence with the color reassignment because the mass of the droplet phase is larger than that of the liquid phase. In addition seven colors are required to change the color of a particle because six particles are close to each other. The spherical region in which the order parameter is exchanged must not contain other droplets of the same color, or else the droplet will wane when the color is changed. The example of this section is only for confirmation although the number of colors is larger than the number of particles. If the number of particles is large, the required number of colors is less than the number of particles. References [1] Batchelor GK, Green JT. The hydrodynamic interaction of two small freelymoving spheres in a linear flow field. J Fluid Mech 1972;56:375–400. [2] Bagster DF, Tomi D. The stresses within a sphere in simple flow fields. Chem Eng Sci 1974;29:1773–83.
[3] Adler PM, Mills PM. Motion and rupture of a porous sphere in a linear flow field. J Rheol 1979;23:25–37. [4] Sonntag RC, Russel WB. Structure and breakup of flocs subjected to fluid stresses. I. Shear experiments. J Colloid Interface Sci 1986;113:399–413. [5] Sonntag RC, Russel WB. Structure and breakup of flocs subjected to fluid stresses. II. Theory. J Colloid Interface Sci 1987;115:378–89. [6] Bossis G, Brady JF. Dynamic simulation of sheared suspensions. I. General method. J Chem Phys 1987;80:5141–54. [7] Harada S, Tanaka R, Nogami H, Sawada M. Dependence of fragmentation behavior of colloidal aggregates on their fractal structure. J Colloid Interface Sci 2006;301:123–9. [8] Fujita M, Yamaguchi Y. Multiscale simulation method for self-organization of nanoparticles in dense suspension. J Comput Phys 2007;223:108–20. [9] Nakayama Y, Yamamoto R. Simulation method to resolve hydrodynamic interactions in colloidal dispersions. Phys Rev E 2005;71:036707. [10] Takeuchi S, Morita I, Kajishima T. Motion of particle agglomerate involving interparticle force in dilute suspension. Powder Technol 2008;184:232–40. [11] Feng ZG, Michaelides EE. The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems. J Comput Phys 2004;195:602–28. [12] Melchionna S. Incorporation of smooth spherical bodies in the lattice Boltzmann method. J Comput Phys 2011;230:3966–76. [13] Tanaka H, Araki T. Simulation method of colloidal suspensions with hydrodynamic interactions: fluid particle dynamics. Phys Rev Lett 2000;85:1338–41. [14] Ladd AJC. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation. J Fluid Mech 1994;271:285–309. [15] Inamuro T, Tomita R, Ogino F. Lattice Boltzmann simulations of drop deformation and breakup in shear flows. Int J Modern Phys B 1997;17:21–6. [16] Inamuro T, Ii T. Lattice Boltzmann simulation of the dispersion of aggregated particles under shear flows. Math Comput Simul 2006;72:141–6. [17] Swift MR, Osborn WR, Yeomans JM. Lattice Boltzmann simulation of nonideal fluids. Phys Rev Lett 1995;75:830–4. [18] Falcucci G, Bella G, Chiatti G, Chibbaro S, Sbragaglia M, Succi S. Lattice Boltzmann models with mid-range interactions. Commun Comput Phys 2007;2:1071–84. [19] Falcucci G, Ubertini S, Succi S. Lattice Boltzmann simulations of phaseseparating flows at large density ratios: the case of doubly-attractive pseudopotentials. Soft Matter 2010;6:4357–65. [20] Sbragaglia M, Benzi R, Biferale L, Succi S, Sugiyama K, Toschi F. Generalized lattice Boltzmann method with multirange pseudopotential. Phys Rev E 2007;75:026702. [21] Zaleski S, Li J, Succi S. Two-dimensional navier-stokes simulation of deformation and breakup of liquid patches. Phys Rev Lett 1995;75:244–7. [22] Inamuro T, Ogata T, Tajima S, Konishi N. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J Comput Phys 2004;198:628–44. [23] Inamuro T. A lattice kinetic scheme for incompressible viscous flows with heat transfer. Philos Trans Roy Soc Lond, Ser A 2002;360:477–84. [24] Sone Y. Molecular gas dynamics. Birkhäuser Boston; 2007. p. 73–167. [25] Leal LG. Advanced transport phenomena: fluid mechanics and convective transport processes. New York: Cambridge University Press; 2007. p. 482–3. [26] Mahanty J, Ninham BW. Dispersion forces. US: Academic Press Inc.; 1976. p. 10. [27] Zick AA, Homsy GM. Stokes flow through periodic arrays of spheres. J Fluid Mech 1982;115:13–26. [28] Iwashita T, Nakayama Y, Yamamoto R. A numerical model for Brownian particles fluctuating in incompressible fluids. J Phys Soc Jpn 2008;77:074007.
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Numerical simulation of the dispersion of aggregated Brownian particles under shear flows Takuya Nishiyama a,⇑, Takaji Inamuro b,c, Shugo Yasuda d a
Electronic and Imaging Materials Research Laboratories, Toray Industries Incorporated, Shiga 520-8558, Japan Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto 606-8501, Japan c Advanced Research Institute of Fluid Science and Engineering, Graduate School of Engineering, Kyoto University, Kyoto 606-8501, Japan d Graduate School of Simulation Studies, University of Hyogo, Kobe 650-0047, Japan b
a r t i c l e
i n f o
Article history: Received 24 October 2011 Received in revised form 27 May 2013 Accepted 24 June 2013 Available online 1 July 2013 Keywords: Particle dispersion Shear flow Van der Waals force Brownian motion Lattice Boltzmann method
a b s t r a c t The deformation and breakup processes of a particle-cluster aggregate under shear flows are numerically investigated by the two-phase lattice Boltzmann method. The van der Waals attraction is considered to be the force between particles. Simulations are performed for various fluid forces acting on particles and various inter-particle forces. It is found that the ratio of the fluid force to the maximum inter-particle force, Y, is a key factor in dispersion, and the aggregate of non-Brownian particles is dispersed when Y is over 0.001. The Péclet number, which is the ratio of the diffusion rate due to shear flow to that due to the Brownian motion, is also considered. By comparing the calculated result of the dispersion of Brownian particles with that of non-Brownian particles, it is found that the Brownian motion impedes dispersion and the effect of the Brownian motion is remarkable when the Péclet number is under 105. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The dispersion of micro- and nanoparticles in a liquid is important for making new functional materials such as ceramics, polymers, and electronic products because the characteristics of materials can be controlled by dispersing small particles in a solvent in their production process. However, small particles are easy to be aggregated by an attraction force between the particles, so it is difficult to disperse a large number of particles uniformly in a liquid. Therefore, to develop new functional materials, it is important to investigate particle dispersion in liquids. The breakup mechanism of aggregated particles has been studied both experimentally and theoretically [1–5]. Because the dispersion dynamics of particles is very complicated, it is difficult to investigate the breakup mechanism of aggregated particles only through experimental and theoretical approaches, so several numerical methods have been proposed. However, in general, there are difficulties in treating a moving solid–liquid boundary even in the numerical simulations of particulate flows. Bossis and Brady [6] proposed the Stokesian dynamics, in which a solvent fluid is not treated explicitly, but the hydrodynamic interaction between particles is treated as a resistance of the solvent fluid expressed by relative velocity vectors of the particles. Harada et al. [7] investigated the dispersion of particles by the ⇑ Corresponding author. E-mail address: [email protected] (T. Nishiyama). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.06.023
Stokesian dynamics. However, the Stokesian dynamics remains difficult for applications with complex boundaries. Immersed boundary approaches, in which a moving boundary condition is satisfied through external forces in the Navier–Stokes equation, have been proposed [8–10]. Feng and Michaelides [11], and Melchionna [12] also proposed an immersed boundary approach combined with the lattice Boltzmann method to solve a fluid–particle interaction problem. These models using the immersed boundary method are effective approaches for moving boundary problems, but the immersed boundary method might cause numerical errors while evaluating the hydrodynamic force and torque acting on particles in various situations. Tanaka and Araki [13] proposed the fluid particle dynamics method, in which a particle is treated as a fluid with large viscosity. Although this method can eliminate the difficulty originating from the solid-liquid boundary condition, it is difficult to retain the spherical shape for a long time. Ladd [14] proposed a simulation model for moving particles using the lattice Boltzmann method in which the hydrodynamic force exerted on a particle is calculated on the basis of momentum exchange. Recently, the two-phase lattice Boltzmann method with the same density [15] has been applied to the simulations of the dispersion of aggregated particles under shear flows by the authors [16]. In this method, the particle is modeled by a hard droplet with large viscosity and strong surface tension, and consequently, we do not need to explicitly track the moving solid–liquid boundary. Strong surface tension is applied to maintain a spherical droplet without any other artificial treatments. Some models of the lattice
396
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
Boltzmann method that involve droplets and tunable surface tension are also available [17–20]. Here the free energy functional scheme proposed by Swift et al. [17] is used to apply strong surface tension. The advantage of the two-phase lattice Boltzmann method is its simplicity in treating the solid-liquid boundary. Zaleski et al. [21], e.g., proposed the volume-of-fluid method, which involves droplets and tunable surface tension. However, the evaluation of the boundary gradient is complex. In the two-phase lattice Boltzmann method, the boundary surface is autonomously determined to minimize the free energy of the system. In addition, a colored order parameter is assigned to each droplet for preventing the droplets from merging into bigger droplets. In this study, we improve the above-mentioned method to treat the Brownian motion of nanoparticles and simulate a large number of particles with a smaller number of colored order parameters. Next, we use the improved method to investigate the behavior of particles under shear flows for various conditions. In general, particles diffuse by the Brownian motion in fluid, so, the Brownian motion may promote dispersion. However, it is also known that particles, which are initially separated, aggregate by the Brownian motion (see e.g., Ref. [16]). Therefore, the effect of the Brownian motion on the dispersion of initially aggregated particles is not yet clear. To classify the calculated results, we introduce two important dimensionless parameters: the ratio of the fluid force to the maximum inter-particle force and the Péclet number, which is the ratio of the diffusion rate due to shear flow to that due to the Brownian motion. The paper is organized as follows: In Section 2, we explain the numerical method. Simulation results are presented in Section 3, in which we investigate the dispersion of the aggregations of 36 and six particles. We conclude the study in Section 4. 2. Numerical method 2.1. Two-phase lattice Boltzmann method for immiscible fluids with the same density Non-dimensional variables are used as in Ref. [22]. The lattice kinetic scheme (LKS) [23], which is an extension method of the lattice Boltzmann methods, is used to formulate the method. In the LKS, macroscopic variables are calculated without particle velocity distribution functions, and thus the scheme can save computer memory, since there is no need to store the particle velocity distribution functions. In addition, in order to represent many hard droplets, which cannot merge into bigger droplets, we introduce colored order parameters to make different colored droplets. Differently colored droplets do not merge if they collide because the boundary surface of each colored droplet is independently and autonomously determined. Note that the color is physically meaningless and is used only for distinguishing each droplet from the others. In the present study, the Stokes flow is assumed because the diameter of a particle is very small (e.g., 1 lm). The 15-velocity model with particle velocities ci(i = 1, 2, . . . , 15) is used in this study. The velocity vectors of this model are given by ½c 1 ; c 2 ; c 3 ; c 4 ; c 5 ; c 6 ; c 7 ; c 8 ; c 9 ; c 10 ; c 11 ; c 12 ; c 13 ; c 14 ; c 15 2 3 0 1 0 0 1 0 0 1 1 1 1 1 1 1 1 6 7 ¼ 4 0 0 1 0 0 1 0 1 1 1 1 1 1 1 1 5: 0 0 0 1 0 0 1 1 1 1 1 1 1 1 1
ð1Þ
i¼1
uðx; t þ DtÞ ¼
ð2Þ
15 1X g eq ðx ci Dx; tÞ; 3 i¼1 i
ð3Þ
15 X c i g eq i ðx c i Dx; tÞ:
ð4Þ
i¼1
where flieq and g eq are equilibrium distribution functions, Dx is a i spacing of the cubic lattice, and Dt is a time step during which the particles travel the lattice spacing. The equilibrium distribution functions in Eqs. (2)–(4) are given by
h i jf flieq ¼ Hi /l þ F i p0 ð/l Þ jf /l r2 /l jr/l j2 þ 3Ei /l cia ua 6 þ Ei jf Gab ð/l Þcia cib ;
ð5Þ
@ub @ua cia cib g eq þ i ¼ Ei 3p þ 3c ia ua þ ADx @xa @xb þ Ei jg Gab ð/l Þcia cib þ 3Ei cia Dx
N X
F la Ul ;
ð6Þ
l¼1
where
E1 ¼ 2=9;
E2 ¼ E3 ¼ E4 ¼ ¼ E7 ¼ 1=9;
E8 ¼ E9 ¼ E10 ¼ ¼ E15 ¼ 1=72; H1 ¼ 1;
H2 ¼ H3 ¼ H4 ¼ ¼ H15 ¼ 0;
F 1 ¼ 7=3;
F i ¼ 3Ei ði ¼ 2; 3; 4; . . . ; 15Þ;
ð7Þ
and
Gab ð/Þ ¼
9 @/ @/ 3 @/ @/ dab ; 2 @xa @xb 2 @xc @xc
ð8Þ
with a, b, c = x, y, z (subscripts a, b, and c represent Cartesian coordinates, and the summation convention is used). In the above equations, dab is the Kronecker delta, jf is a constant parameter determining the width of the interface, jg is a constant parameter determining the strength of the surface tension, and A is a constant parameter related to fluid viscosity. In Eq. (4), p0(/) is given by
p0 ð/Þ ¼ /T /
1 a/2 ; 1 b/
ð9Þ
where a, b, and T/ are free parameters determining the maximum and minimum values of the order parameter, /max and /min. It is noted that flieq is the same as that for the model of Swift et al. [17], except that flieq in Eq. (5) has no second terms of ua because we assume the Stokes flow. The last term of Eq. (6) is the force Fl acting on the lth particle per unit mass:
Fl ¼
N X
F v lm þ F Br l ;
ð10Þ
m¼1
where Fv lm is the attractive force per unit mass from the mth particle to the lth particle, and FBr l is the Brownian force per unit mass acting on the lth particle. The detailed formulations of forces are given in the following subsections. In addition, Ul in Eq. (6) represents the region of a droplet; i.e., it is unity inside the droplet and zero outside the droplet:
The physical space is divided into a cubic lattice, and the colored order parameter /l(x, t) (l = 1, 2, . . . , N), where N is the number of colors, the pressure p(x, t) and the velocity u(x, t) of whole fluid at the lattice point x and at time t are computed as follows: 15 X /l ðx; t þ DtÞ ¼ flieq ðx c i Dx; tÞ;
pðx; t þ DtÞ ¼
Ul ¼
1; /l P /B ; 0; /l < /B :
ð11Þ
/B is the threshold value of the boundary and is defined as follows:
/B ¼
f/min þ /max ; fþ1
ð12Þ
where f is a weight parameter and is described in detail in Appendix A.
397
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
The following finite-difference approximations are used to calculate the derivatives in Eqs. (5), (6), and (8): 15 @w 1 X cia wðx þ ci DxÞ; @xa 10Dx i¼1 " # 15 X 1 2 r w wðx þ c i DxÞ 14wðxÞ : 5ðDxÞ2 i¼2
ð13Þ ð14Þ
2.2. Inter-particle force
Applying the asymptotic theory [24] to Eqs. (2)–(9), it is found that the macroscopic variables (pressure p and velocity u) satisfy the incompressible Navier–Stokes equation with relative errors of O[(Dx)2] [22]. The kinematic viscosity, m, of the fluid is given by
1 2 m ¼ A Dx: 6 9
ð15Þ
Thus, kinematic viscosity can be changed the by A. In the interface, kinematic viscosity is changed as follows:
m ¼ mc þ
mp mc 2
/0 /B p þ1 ; sin /max /min
ð16Þ
where mp is the kinematic viscosity of the droplet, mc is that of the solvent, and /0 ¼ maxl2f1;2;...;Ng /l at the lattice point. We set the viscosity ratio of the droplet to the surrounding liquid to be g = md/ mc = 10, and the droplet and the liquid have equal density. According to the Hadamard–Rybczynski drag law [25], the ratio of the drag force for a spherical fluid particle to that for a rigid spherical particle is 0.970. Thus, the Stokes drag law is recovered within 3% error. The interfacial tension, rl, for the lth colored phase is given by
rl ¼ jg
2 @/l dn; @n 1
Z
physically meaningful current is set to be sufficiently large than the spurious current by tuning the parameters. The parameters we choose are described in Sections 3.1 and 3.2. By using the above method, we can form many hard droplets with large viscosity and strong surface tension to model solid particles. The validity of the method is examined in Appendix A.
The van der Waals attractive force is applied as the inter-particle force. The attractive force per unit mass from the mth particle to the lth particle Fvlm is given by [26]
Fv
lm
¼
8 > Ah > > < pD3 q
D2 r lm 2 ðr 2lm D2 Þ
p
> > > :
2
2 lm þ rD3 r22rD 2 þ r lm
lm
for D þ LD 6 r lm 6 2D; 0;
lm
r lm r lm
; ð18Þ
otherwise;
where D is the diameter of particle, qp is the density of the particle, rlm is the vector from the center of the lth particle to that of the mth particle, rlm = jrlmj, Ah is the Hamaker constant, and LD corresponds to the Debye length, which gives the thickness from the surface of the particle. The attractive force reduces to zero, if the distance between particle surfaces is smaller than LD. In addition, to avoid allto-all computations, the inter-particle force is cut off when the distance between the centers of the two particles is greater than 2D, because the van der Waals force is very small when two particles are far apart. The correction of the lubrication force between particles is not considered in the present study because its effect is negligible in the dispersion process of aggregated particles.
1
ð17Þ
with n being the coordinate normal to the interface. It is noted that a small spurious current is generated around the droplet when strong surface tension is applied. However, in the simulation,
2.3. Brownian motion We apply the Brownian motion to particles using the Langevin approach. The advantages of the Langevin approach are its simplicity and computational efficiency.
φ1 φ1 φ1 φ1
(a)
(b) Fig. 1. Color reassignment. (a) The same color is assigned to particles when they are far apart from each other. (b) Exchange of order parameters in certain areas.
398
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
The Brownian force acting on the lth particle per unit mass at every time step is formularized as follows:
F Br la
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6pqc mc DkB T Xa; ¼ Dt pD3 qp 6
ð19Þ
where kB is the Boltzmann constant, qc is the density of the solvent, T is the temperature of the solvent, and Xa is the normal random number whose average is 0 and variance is 1. Note that Xa is independently chosen for a = x, y, and z. We also note that T is not the temperature of the solvent in the experimental system but the temperature for a single Brownian particle obeying the Langevin equation in an infinite domain. The temperature, T, is used as a reference value of the solvent temperature in the present simulation. The difference between the temperatures T and Tc estimated from a numerical simulation of particle diffusion is given in Appendix B. Fig. 3. Computational domain and six aggregated particles.
2.4. Color reassignment As described in Section 2.1, we assign a different colored order parameter to each particle to prevent the particles from merging into a larger particle. However, if the number of particles increases, many colors are necessary, which increases computational time. To save computational time, the same color is assigned to particles that are far apart from each other, as shown in Fig. 1(a). Next, if particles of the same color approach each other, one particle’s color is changed. To change the color of the particle /1 to /2, the order parameter /1 and /2 are exchanged in the region within a specified distance (e.g., 1.6D/2) from the center of mass of the particle, as shown in Fig. 1(b). The accuracy of this scheme is examined in Appendix C. 3. Results and discussion In our previous paper [16], four cases with two, six, 18, and 36 non-Brownian particles are calculated. In that study, the aggregates with 18 and 36 particles are loose structures. In the present study, an aggregate composed of particles with hexagonal closepacked structure is used. In the following subsection, we compare the present results to the previous results from the view point of the effect of the initial arrangement of particles. Aggregated particles composed of same particles with diameter D are placed in a liquid inside a rectangular domain. To evaluate the cluster size effect, two different aggregations are considered: the aggregation of 36 particles (Fig. 2) and the aggregation of six particles (Fig. 3). Initially, the particles are arranged in a hexagonal close-packed structure for both aggregations. In addition, for the
aggregation of six particles, two domain sizes are used to investigate the domain size effect. At t = 0, the top and bottom walls begin to move in the x-direction with velocities uw and uw, respectively. Periodic boundary conditions are used on the other sides of the domain. The solid–liquid boundary is determined by the threshold value (/max + /min)/2, i.e., f = 1 in Eq. (12), so the diameter is determined by D1/2 (the definition of D1/2 is given in Appendix A). Although D1/4 is more suitable, as observed in Appendix A, we use D1/2 throughout this study for simplicity. The other parameters are set at the same values as those mentioned in Ref. [16]. In order to classify the calculated results, we introduce the following non-dimensional parameters: (a) The ratio of the fluid force [/qcmc(uwD/Lz)D] to the inter-particle force [/Ah/D]:
I¼
qc mc uw D3 Ah Lz
ð20Þ
:
(b) The ratio of LD to the diameter of the particle:
s ¼ LD =D:
ð21Þ
(c) The ratio of the fluid force [/qcmc(uwD/Lz)D] to the maximum inter-particle force [/Ah/s2D for s 1]:
Y¼
qc mc uw D3 s2 Ah Lz
¼ Is2 :
ð22Þ
Note that the inter-particle force is maximum when rlm = D + LD in Eq. (18). Table 1 Computational conditions.
Fig. 2. Computational domain and 36 aggregated particles.
Case Case Case Case Case Case Case Case Case Case Case Case Case Case Case
A B C D E F G H I J K L M N O
I
s
Y
Pe
1000 100 10 10 10 1000 100 100 0.625 10 2.5 1000 100 10 0.625
0.005 0.005 0.005 0.016 0.05 0.016 0.05 0.016 0.2 0.01 0.01 0.005 0.016 0.05 0.2
0.025 0.0025 0.00025 0.0025 0.025 0.25 0.25 0.025 0.025 0.001 0.00025 0.025 0.025 0.025 0.025
1 1 1 1 1 1 1 1 1 1 1 350,000 35,000 3500 220
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
(d) The Péclet number, which is the ratio of the rate of diffusion due to shear flow[/2uwD/Lz] to the rate of diffusion due to the Brownian motion [/kBT/3pqcmcD2]:
Pe ¼
6pqc mc uw D3 6pAh I: ¼ kB TLz kB T
ð23Þ
The computational conditions are given in Table 1. All conditions are implemented for the aggregation of 36 particles, whereas
399
Cases A–K (non-Brownian cases) are implemented for the aggregation of six particles. In Cases A–K, the Brownian motion is not considered. In Cases A–C, the effect of the intensity of the shear flow is investigated. In Cases C–E, the effect of the value of LD is investigated. In Cases F–K, the intensity of the shear flow and LD are varied. In Cases L–O, the Brownian motion is applied to investigate the dispersion of nanoparticles. In reality, Case A, for example, corresponds to the following situation, in which the parameters with a tilde are dimensional parameters in the real system. The particle is made of titanium oxide with a
Fig. 4. Time evolution of the aggregate under shear flows for (a) I = 1000, s = 0.005, Y = 0.025, and Pe = 1 (Case A); and (b) I = 10, s = 0.005, Y = 0.00025, and Pe = 1 (Case C). t⁄ = 2uwt/Lz.
400
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
e h ¼ 8:0 1020 J, diameter D e ¼ 1 lm, and Hamaker constant of A e L D ¼ 5 nm. The solvent is water with a kinematic viscosity of m~c ¼ 1:0 106 m2 s1 . The top and bottom walls at a distance of 10 lm move in opposite directions with a velocity of ~ w ¼ 0:8 m s1 . Moreover, Case N corresponds to the following situau tion. The particle is titanium oxide with a Hamaker constant of e h ¼ 8:0 1020 J, diameter D e ¼ 100 nm, and e A L D ¼ 5 nm. The sol~c ¼ 1:0 106 m2 s1 vent is water with a kinematic viscosity of m and a temperature of 308 K. The top and bottom walls at a distance ~ w ¼ 0:8 m s1 . of 1 lm move in opposite directions with a velocity of u
3.1. Results for 36 particles The domain is divided into a 240 100 200 cubic lattice. The diameter of the particle is D = 20Dx. The 18–20 colored order parameters are used to save computation time. We chose a = 9/ 49, b = 2/21, and T/ = 0.55 in Eq. (9). The other parameters are fixed at jf = 0.01Dx2, jg = 0.01Dx2, and uw = 0.1. The calculated results of Case A, C, E, and N are shown in Figs. 4 and 5. In these figures, t⁄ = 2uwt/Lz is dimensionless time. In these calculations, color reassignments were implemented 72, 39, 73, and 64 times for Cases A,
Fig. 5. Time evolution of the aggregate under shear flows for (a) I = 10, s = 0.05, Y = 0.025, and Pe = 1 (Case E) and (b) I = 10, s = 0.05, Y = 0.025, and Pe = 3500 (Case N). t⁄ = 2uwt/Lz.
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
Fig. 6. Time evolution of the standard deviation, r, for the aggregation of 36 particles (I = 1000, Y = 0.025, and Pe = 1 for Case A; I = 10, Y = 0.00025, and Pe = 1 for Case C; I = 10, Y = 0.025, and Pe = 1 for Case E; and I = 10, Y = 0.025, and Pe = 3500 for Case N). t⁄ = 2uwt/Lz.
401
Fig. 9. Time evolution of the standard deviation, r, for the aggregation of six particles in domain A. r0 is the initial standard deviation, and t⁄ = 2uwt/Lz.
Fig. 10. Time evolution of the standard deviation, r, for the aggregation of six particles in domain B. r0 is the initial standard deviation, and t⁄ = 2uwt/Lz. , versus Y for the aggregation of 36 Fig. 7. Time-averaged standard deviation, r is the time-averaged value for 140 6 t⁄ 6 164, and r0 is the particles for Pe = 1. r initial standard deviation.
1.2
1
Case L / Case A Case M / Case H Case N / Case E Case O / Case I
0.8
0.6
0.4 102
103
104
105
106
Pe , to Fig. 8. Ratio of the time-averaged standard deviation of the Brownian case, r ðPe¼1Þ , at Y = 0.025 for the aggregation of 36 that of the non-Brownian case, r ⁄ and r ðPe¼1Þ are time-averaged values for 140 6 t 6 164. particles. r
C, E, and N, respectively. From the results of Figs. 4 and 5, we find that the aggregated particles are dispersed by shear flow for Case A, whereas they are not dispersed for Case C because the shear flow is
weak (the value of I in Case C is one-hundredth of that in Case A). However, the aggregated particles are dispersed for Case E although I has the same value as in Case C because the value of s in Case E is ten times larger than that in Case C. In Case N, the Brownian motion is considered, but other conditions are the same as those in Case E. Comparing the results of Case E with those of Case N, it is observed that the Brownian motion impedes particle dispersion. To measure the dispersion state of the particles, the standard deviation r of the lengths between the center of the aggregation and those of the particles is calculated during the deformation of the aggregates. The time evolutions of r for Cases A, C, E, and N are shown in Fig. 6. In this figure, r0 is the standard deviation for the initial structure. The results of Fig. 6 show that the standard deviations for Cases A and E are close each other, which indicates that the parameter Y, which is a combination of I and s, should be a key factor determining particle dispersion. I is the parameter relating two particles far apart from each other, and Y is the parameter relating two particles that are almost in contact. On the other hand, it is clearly found from Fig. 6 that the standard deviation of 36 particles for Case N increases more slowly than that for Case E, which means that the Brownian motion impedes particle dispersion and that the dispersion of nanoparticles is more difficult. , for Fig. 7 shows the time-averaged standard deviation, r 140 6 t⁄ 6 164 for all cases with 36 non-Brownian particles from Table 1. The results in Fig. 7 show that dispersion begins when Y
402
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
is over 0.001. In addition, our previous results [16] indicate that aggregates with a loose structure begin to disperse when Y is over 0.0025, which is of the same order as 0.001, although only four cases for s = 0.005 are available. Thus, we consider that for dispersion, the initial particle arrangement does not considerably change the threshold value of Y. Fig. 8 shows the relationship between the Péclet number and time-averaged standard deviation for 140 6 t⁄ 6 164 of Brownian particles normalized by that of non-Brownian particles at the same value of Y = 0.025. It is found that the dispersion-retarding effect of the Brownian motion is remarkable when the Péclet number is under 105. 3.2. Results for 6 particles In order to examine the cluster- and domain-size effects, nonBrownian cases (Cases A–K in Table 1) are simulated for the aggregation of 6 particles with two different domain sizes. One domain is divided into a 132 80 110 cubic lattice, and the other is divided into a 198 120 165 cubic lattice. We call the former domain ‘‘domain A’’ and the latter ‘‘domain B.’’ The diameter of the particle is D = 20Dx. We use the same values for the parameters a, b, T/, jf, and jg as those mentioned in Section 3.1 and uw = 0.055 and 0.0825 for domains A and B, respectively. The time evolution of the standard deviation, r, for domain A (B) is shown in Fig. 9 (10). We can observe that the particles disperse in all cases except Cases C and K in which Y is over 0.001 for both domains A and B. Thus, the threshold value of Y for dispersion is 0.001 for six particles, which is the same as that for 36 particles, and the cluster- and domain-size effects are not clearly observed. It is noted that the dispersion of particles is slower in domain B, which is larger than domain A, because the development of shear flows depends on the height Lz of the domain.
Fig. A.1. Computational domain for the measurement of Stokes drag on a particle.
D) is placed at the center of a cubic domain of size L, as shown in Fig. A.1. Because the interface of the droplet is approximately 3Dx in thickness, we define the diameter of the droplet by using the threshold value, /B, given in Eq. (12). We denote the diameter by D1/1+f, which corresponds to the weight parameter, f. The usual periodic boundary conditions are used on all the sides of the domain, and periodic boundary conditions with pressure difference Dp are used at x = 0 and x = L. Because of the flow induced by the pressure difference, the hard droplet moves in the x-direction. We apply a uniform body force to the internal region of the droplet to keep the droplet stationary. The body force is applied to the region / > (9/max + /min)/10. After a steady state is obtained, we calculate the total body force, F, and the average velocity, U, at x = 0. Computational conditions and calculated results are given in Table A.1, in which the deformation factor, k1/1+f, is also shown. The droplet is slightly deformed by the flow, and we measure the deformation factor defined by
4. Concluding remarks To numerically investigate the dispersion of aggregated nanoparticles, we applied the lattice Boltzmann method for the twophase fluid flow with uniform density for simulating the dispersion of aggregated particles under shear flow. The ratio of fluid force to the maximum inter-particle force, Y, is introduced, and it is found that Y is the key factor for dispersion and the aggregates of nonBrownian particles are dispersed when Y is over 0.001. By comparing the calculated results of the dispersion of Brownian particles to those of non-Brownian particles, it is found that the Brownian motion impedes dispersion and the effect of the Brownian motion is remarkable when the Péclet number is under 105.
k1=1þf
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 !2 u I u D DII1=1þf 1=1þf t ¼ 1 þ 1 D1=1þf D1=1þf
ðA:1Þ
where DI1=1þf and DII1=1þf are diameters in the x- and the z-directions, respectively. It is found from Table A.1 that k1/1+f is greater than 5% in the smallest domain (Case I), but it is relatively small for the other cases. Fig. A.2 shows the calculated results of the relationship between the dimensionless drag force, Q, and the filling factor, Ct, which are defined as follows:
Appendix A. Validation of present method: drag force acting on a sphere
Q¼
3plD1=1þf U F
ðA:2Þ
and
We calculated a Stokes flow through the periodic arrays of spheres and compared the drag force acting on the sphere with the results of Zick and Homsy [27]. A hard droplet with an initial diameter D (initially, we assign /max inside D, and /min outside
Ct ¼
pD31=1þf =6 L3
ðA:3Þ
:
Table A.1 Computational conditions and calculated results for the measurement of Stokes drag. L Case Case Case Case Case Case Case
I II III IV V VI VII
31 35 41 45 60 70 90
Dp 3
1.3 10 7.5 104 4.0 104 2.5 104 7.0 104 3.5 104 1.5 104
U
F
D1/10
D1/4
D1/2
k1/10 (%)
k1/4 (%)
k1/2 (%)
0.067 0.066 0.066 0.069 0.070 0.067 0.066
1.25 0.919 0.671 0.607 2.52 1.71 1.19
22.6 22.5 22.4 23.4 43.2 43.4 44.7
21.2 21.0 21.0 20.8 41.8 42.3 43.4
19.9 19.7 19.6 19.5 40.5 40.9 42.1
5.6 2.6 1.7 1.3 2.5 1.3 0.32
5.4 3.1 1.9 1.2 3.2 0.92 0.66
5.3 2.4 1.1 1.1 2.6 1.0 0.34
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
403
Fig. B.1. Time evolution of mean-square displacement hjxBrj2i; t is the time step.
Fig. A.2. Comparison between the results of Zick and Homsy and those of the present method for various /B. Q is the dimensionless drag force, and Ct is the filling factor.
It is found from Fig. A.2 that the results estimated using D1/4 agree well with those of Zick and Homsy regardless of the droplet size. Appendix B. Validation of the present method: free diffusion of particle We calculated the Brownian motion of a single spherical particle in a three-dimensional cubic domain and compared the estimated temperature, Tc, with the given temperature, T, from the fluctuation–dissipation relation. A hard droplet with diameter D = 21Dx (D1/4) is placed at the center of a cubic domain of size
L = 160Dx. The usual periodic boundary conditions are used on all sides of the domain. The same random force per unit mass as Eq. (19) (except for Q) is applied to the particle at each time step:
F Br la
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6pqc mc DkB T Xa: ¼ Q Dt pD3 6
ðB:1Þ
We use kBT = 4.4 103, qcmc = 0.0167, and Q = 0.815. Here Q is defined by Eq. (A.2) as the ratio of the drag force acting on the particle in the periodic domain to that in the infinite domain. In addition, the additional Q in Eq. (B.1) is introduced in the same manner as that in Ref. [28]. However, we assume Q = 1 in the simulation of the dispersion of particles, as mentioned in Eq. (19). Fig. B.1 shows the time evolution of the mean-square displacement. It is found that the mean-square displacement is proportional to the time step, t, when t is over 3 104. The calculated
Fig. C.1. Time evolution of the aggregation of six particles under shear flows for Case A using (a) six colors and (b) seven colors with color reassignment; t⁄ = 2uwt/Lz. (c) Time evolution of total mass (enlarged view). (d) Time evolution of total mass (overall view).
404
T. Nishiyama et al. / Computers & Fluids 86 (2013) 395–404
diffusion coefficient is Dc0 ¼ 1:2 103 . The temperature of the system, Tc, estimated from the Stokes–Einstein relation is c kT ¼ 3pqc mc DDc0 =Q ¼ 4:9 103 . Therefore, the fluctuation–dissipation relation is almost satisfied; i.e. the estimated temperature, Tc, is 12% larger than the given temperature, T. Appendix C. Validation of the present method: color scheme We calculated the dispersion of six particles by using seven colors with the color reassignment scheme and compared the results with those obtained when the color reassignment scheme is not used. The domain is the same as domain A (132 80 110 cubic lattice), and the computational condition is the same as Case A from Section 3.2. Initially the particle color is the same as the particle number. The color of the sixth particle alternates between the sixth or seventh color for every 100 steps. The colors of other particles do not change. When the color of the sixth particle changes, the order parameters /6 and /7 are exchanged in the region within 16Dx from its center of mass. Figs. C.1(a) and (b) show the time evolution of the aggregate when the color reassignment scheme is not used and when seven colors are used with the color reassignment scheme, respectively. In the figure, t is the time step, and t⁄ = 2uwt/Lz is the dimensionless time. The two results are almost the same. In addition, Fig. C.1(c) and (d) show the time evolution of the total mass of the order parameters /6 and /7 with seven colors and the color reassignment scheme, respectively. It is evident that the total mass is conserved. Note that Fig. C.1(c) shows that the total mass varies over 100 steps in correspondence with the color reassignment because the mass of the droplet phase is larger than that of the liquid phase. In addition seven colors are required to change the color of a particle because six particles are close to each other. The spherical region in which the order parameter is exchanged must not contain other droplets of the same color, or else the droplet will wane when the color is changed. The example of this section is only for confirmation although the number of colors is larger than the number of particles. If the number of particles is large, the required number of colors is less than the number of particles. References [1] Batchelor GK, Green JT. The hydrodynamic interaction of two small freelymoving spheres in a linear flow field. J Fluid Mech 1972;56:375–400. [2] Bagster DF, Tomi D. The stresses within a sphere in simple flow fields. Chem Eng Sci 1974;29:1773–83.
[3] Adler PM, Mills PM. Motion and rupture of a porous sphere in a linear flow field. J Rheol 1979;23:25–37. [4] Sonntag RC, Russel WB. Structure and breakup of flocs subjected to fluid stresses. I. Shear experiments. J Colloid Interface Sci 1986;113:399–413. [5] Sonntag RC, Russel WB. Structure and breakup of flocs subjected to fluid stresses. II. Theory. J Colloid Interface Sci 1987;115:378–89. [6] Bossis G, Brady JF. Dynamic simulation of sheared suspensions. I. General method. J Chem Phys 1987;80:5141–54. [7] Harada S, Tanaka R, Nogami H, Sawada M. Dependence of fragmentation behavior of colloidal aggregates on their fractal structure. J Colloid Interface Sci 2006;301:123–9. [8] Fujita M, Yamaguchi Y. Multiscale simulation method for self-organization of nanoparticles in dense suspension. J Comput Phys 2007;223:108–20. [9] Nakayama Y, Yamamoto R. Simulation method to resolve hydrodynamic interactions in colloidal dispersions. Phys Rev E 2005;71:036707. [10] Takeuchi S, Morita I, Kajishima T. Motion of particle agglomerate involving interparticle force in dilute suspension. Powder Technol 2008;184:232–40. [11] Feng ZG, Michaelides EE. The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems. J Comput Phys 2004;195:602–28. [12] Melchionna S. Incorporation of smooth spherical bodies in the lattice Boltzmann method. J Comput Phys 2011;230:3966–76. [13] Tanaka H, Araki T. Simulation method of colloidal suspensions with hydrodynamic interactions: fluid particle dynamics. Phys Rev Lett 2000;85:1338–41. [14] Ladd AJC. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation. J Fluid Mech 1994;271:285–309. [15] Inamuro T, Tomita R, Ogino F. Lattice Boltzmann simulations of drop deformation and breakup in shear flows. Int J Modern Phys B 1997;17:21–6. [16] Inamuro T, Ii T. Lattice Boltzmann simulation of the dispersion of aggregated particles under shear flows. Math Comput Simul 2006;72:141–6. [17] Swift MR, Osborn WR, Yeomans JM. Lattice Boltzmann simulation of nonideal fluids. Phys Rev Lett 1995;75:830–4. [18] Falcucci G, Bella G, Chiatti G, Chibbaro S, Sbragaglia M, Succi S. Lattice Boltzmann models with mid-range interactions. Commun Comput Phys 2007;2:1071–84. [19] Falcucci G, Ubertini S, Succi S. Lattice Boltzmann simulations of phaseseparating flows at large density ratios: the case of doubly-attractive pseudopotentials. Soft Matter 2010;6:4357–65. [20] Sbragaglia M, Benzi R, Biferale L, Succi S, Sugiyama K, Toschi F. Generalized lattice Boltzmann method with multirange pseudopotential. Phys Rev E 2007;75:026702. [21] Zaleski S, Li J, Succi S. Two-dimensional navier-stokes simulation of deformation and breakup of liquid patches. Phys Rev Lett 1995;75:244–7. [22] Inamuro T, Ogata T, Tajima S, Konishi N. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J Comput Phys 2004;198:628–44. [23] Inamuro T. A lattice kinetic scheme for incompressible viscous flows with heat transfer. Philos Trans Roy Soc Lond, Ser A 2002;360:477–84. [24] Sone Y. Molecular gas dynamics. Birkhäuser Boston; 2007. p. 73–167. [25] Leal LG. Advanced transport phenomena: fluid mechanics and convective transport processes. New York: Cambridge University Press; 2007. p. 482–3. [26] Mahanty J, Ninham BW. Dispersion forces. US: Academic Press Inc.; 1976. p. 10. [27] Zick AA, Homsy GM. Stokes flow through periodic arrays of spheres. J Fluid Mech 1982;115:13–26. [28] Iwashita T, Nakayama Y, Yamamoto R. A numerical model for Brownian particles fluctuating in incompressible fluids. J Phys Soc Jpn 2008;77:074007.