6342
Langmuir 2000, 16, 6342-6350
Mesoscopic Simulation of Drops in Gravitational and Shear Fields Alec T. Clark, Moti Lal,* J. Noel Ruddock, and Patrick B. Warren Unilever Research Port Sunlight, Bebington, Wirral CH63 3JW, England Received December 2, 1999. In Final Form: March 7, 2000 In seeking validation of Dissipative Particle Dynamics (DPD) for the mesoscopic modeling of multiphase fluid-fluid systems in external fields, simulations of a pendant drop and a drop in simple shear flow have been performed. The shape profile of the simulated pendant drop was found to be in perfect agreement with that computed by solving the Laplace equation. At increased values of the gravitational force (g), the drop underwent considerable elongation, developing a “neck” between the solid support and its bulk part. Further increases in g resulted in thinning of the neck, which ruptured as g exceeded a certain value, leading to the detachment of the drop. This picture of the detachment process is consistent with the experimental observations published in the literature. Also, the simulations reproduced the drop volume experiment quantitatively. For the drop in shear flow, the degree of deformation was found to be a linear function of the capillary number (Ca) in the region Ca e 0.11, in good agreement with Taylor’s theory; this is despite the fact that the hydrodynamic regime in the simulations (Re ∼ 1-10) is quite different from that assumed in the theory (Re , 1). At increased shear rates the results showed positive departure from linearity, in agreement with theory and experiment. Further increases in Ca resulted in the drop assuming a dumbbell like shape, the middle part of the “dumbbell” gradually stretching to form a thin neck. The rupture of the neck was occasioned by the instabilities manifested in the form of stochastic oscillations which magnified as the critical point was reached. The time evolution of the shape of the drop as it underwent the breakup process in our simulations bears remarkable similarity to the experimental observations of Torza et al. The critical value of the capillary number obtained in the simulations is in reasonable agreement with the experimental figure.
1. Introduction In the simulation of complex systems, it is often required to cover length and time scales that are orders of magnitude greater than those corresponding to the molecular dimensions and the characteristic times of ordinary molecular processes. Access to these domains, denied in MD simulations due to totally unrealistic computational demands, has been rendered possible by the advent of simulation approaches based on mesoscopic models.1 In a mesoscopic model, the structural unit is a coarse-grained representation of a very large number of molecules. The static and dynamic interactions between these units are assumed to be given by simple functional forms, the dynamics of the system being governed by the forces deriving from these interactions as well as by any external forces acting on it. Such models are neither molecular (all molecular details are subsumed into the structural unit) nor continuum, but occupy a middle (meso) ground between the two extreme descriptions of materials. Consideration of these kinds of models would be appropriate when it is sufficient to define the system with a minimal requirement of information at the molecular level. Familiar examples of complex materials include colloids, emulsions, gels, microporous media, and composites. They are multiphase, multicomponent systems invariably con* Corresponding author. Also at the Centre for Nano-scale Science, The Donnan Laboratories, University of Liverpool, Liverpool L69 7ZD, England. (1) Manneville, P., Boccara, N., Vichiniac, Y. G., Bidaux, R., Eds.; Cellular Automata and Modelling of Complex Systems; Springer: Heidelberg, 1990. (b) Lal, M., Mashelkar, R. A., Kulkarni, B. D., Naik, V. M., Eds.; Structure and Dynamics of Materials in the Mesoscopic Domain; Imperial College Press: 1999. (c) Faraday Discussion 112: Physical Chemistry in the Mesoscopic Regime; The Royal Society of Chemistry: London, 1999. (d) Higgins, J. S., Lal, M., Tildesley, D. J., Eds.; Mesoscale Perspectives of Polymers, Colloids and SurfactantssA special issue of PCCP 1999, 1 (9), 2039-2202.
taining at least one interface (liquid/solid, liquid/liquid, liquid/air, or solid/solid) and often involving external fields such as gravity and flow. Although the mesoscopic simulation does appear to be an effective approach for modeling their behavior, it is imperative that we validate this approach for its ability to treat the various interfaces and the external fields correctly. The need for validation is rendered all the more compelling in view of a certain lack of rigor in the coarse graining of the models and the somewhat empirical nature of the forces involved. Among the various methods in current use, a particle dynamics method, known as Dissipative Particle Dynamics (DPD), promises to be quite apt in modeling systems comprised of immiscible fluids and solid boundaries.2 Our purpose here is to pursue the validation of this method for its application to two-fluid systems in external fields. In the selection of systems for our present purpose, two basic criteria have to be applied. First, the model systems must possess essential features of a real, complex system; that is, they contain at least two phases, an interface and an external field. Second, accurate analytical and/or experimental results are available against which to evaluate the accuracy of the simulation results. Two systems which meet these criteria adequately are (i) a drop hanging from a solid support under a downward (gravitational) force and (ii) a drop in a simple shear field, the subjects of the present study. 2. Dissipative Particle Dynamics (DPD) 2.1. Model. The DPD model, which was originally proposed by Hoogerbrugge and Koelman,3 assumes that (2) Jones, J. L.; Lal, M.; Ruddock, J. N.; Spenley, N. J. Faraday Discuss. 1999, 112, 129. (3) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, 19, 155. (b) Koelman, J. M. V. A.; Hoogerbrugge, P. J. Europhys. Lett. 1993, 21, 363.
10.1021/la991565f CCC: $19.00 © 2000 American Chemical Society Published on Web 04/29/2000
Mesoscopic Simulation of Drops
Langmuir, Vol. 16, No. 15, 2000 6343
the dynamic interactions between the mesoscale units, termed momentum carriers, are composed of two parts, dissipative and stochastic, complementing each other to ensure a constant value for the mean kinetic energy, kT, of the system. The dissipative part, resisting the motion of the momentum carriers, is assumed to be proportional to their relative velocity:
FijD ) -ζij(vi - vj)
(1)
ζij ) γωijD(rij) eijeij
(1a)
FijC ) RωijC(rij)eij
where the coefficient R, a positive quantity, is a measure of the strength of the conservative force and ωijC(rij) is the weight function for the dependence of FijC on rij. The functional forms for ωijR and ωijC prescribed in the model are linear in rij with a truncation distance, r*, which serves as the unit of distance; that is, we set r* ) 1:
ωijR ) ωijC ) (1 - rij)
with
)0
FijD ) - γωijD[eij(vi - vj)]eij
(2)
The dissipative force, since it acts against the particle motion, would tend to reduce the kinetic energy of the system, which is compensated by the random motion produced by the stochastic part, FijR. FijR is assumed to be Brownian-like; that is, it follows a Gaussian distribution with the first moment equal to zero and the second moment proportional to the temperature. So, FijR is expressible as
FijR ∼ σξijeij
FijR ) ωijR(rij)σξijeij
(4)
σ2 ) kT, the kinetic energy of the system 2γ
(5)
(ωijR)2 ) ωijD
(6)
with
With kT as the unit of energy,
for rij g 1
(8)
for rij < 1 for rij g 1
(8a)
It may be noted that the functional form of ωijC corresponds to a repulsive harmonic potential, a kind of soft sphere potential, assumed to result from the coarse graining of the intermolecular interactions between the momentum carriers. No rigorous justification for this assumption exists. However, the coarse-graining simulations of Forrest and Suter,5,6 though less than rigorous, show that the resulting mesoscale potential is monotonically repulsive and nondivergent in the limit rij f 0. A further assumption of the model is that all three forces are pairwise additive, which, in combination with eqs 2 through 8 and 8a, leads to the following expression for the force acting on i due to all other particles (momentum carriers) in the system N
Fi )
N
N
FijD + ∑FijR + ∑FijC ∑ j>i j>i j>i N
- γ(1 - rij)2[eij(vi - vj)]eij + ∑ j>i,r i,r 0.12, the degree of deformation, , deviates appreciably from the linear behavior, the deviation diverging with an increase in Ca, a trend confirmed by the experimental results of Guido and Villone for ν e 1.522 (ν ) 1 for the present system). The theoretical treatments for this regime, which are based on perturbation expansion approaches,24 also predict the positive departure from linearity at ν ) 1, agreeing with the present simulation results. Another possible cause for this deviation is the effect of the Reynolds number. At larger Reynolds numbers, one might expect a tendency for the hydrostatic pressure in the fluid to be reduced in regions where the fluid flows more rapidly (Bernouilli principle). Since this occurs at the tips of the droplet, one might expect the droplet to become more extended. We note that the ratio Re/Ca ) ΓRF/η2 is independent of the strength of the flow field. In a real system, say, water/oil, η ) 10-3 Pa s, F ) 103 kg m-3, and Γ ) 10-2 Nm-1, this ratio is Re/Ca ∼ 107R. Thus, it is only for nanometer size droplets that significant distortion can occur (Ca ∼ 1) while remaining in the Stokes limit (Re , 1). Our present simulations correspond to Re/Ca ≈ 45, equivalent to a 4.5 µm radius droplet. 4.2.2. Drop Breakup. The drop continued to elongate with the increase in shear rate and reached a state when its ellipsoidal shape began to distort. The distortion magnified with further increases in γ˘ , eventually leading to a critical stage at which the drop was unable to sustain its integrity and broke up into two or more fragments. (22) Guido, S.; Villone, M. J. Rheol. 1998, 42, 395. (23) Phillips, W. J.; Graves, R. W.; Flumerfelt, R. W. J. Colloid Interface Sci. 1980, 76, 350. (24) Barthes-Biesel, D.; Acrivos, A. J. Fluid Mech. 1973, 61, 1. (b) Rallison, J. M. Annu. Rev. Fluid Mech. 1984, 16, 45. (c) Stone, H. A. Annu. Rev. Fluid Mech. 1994, 26, 65.
Mesoscopic Simulation of Drops
Langmuir, Vol. 16, No. 15, 2000 6349
Figure 8. Sequence of eight snapshots depicting the drop breakup process in a simulation run performed at Ca ) 0.27 (slightly above Cac). The snapshots are taken successively at 157.50, 161.75, 165.75, 168.50, 170.00, 172.50, 173.50, and 174.25 DPD units of time.
The breakup process revealed by the present simulation, at Ca ) 0.27 (slightly above the critical value), is presented in Figure 8 as a sequence of snapshots covering the last stages of the process. We find that as the drop continues to deform, it begins to develop a “neck”. The neck undergoes thinning with time, and occasionally we see oscillations occurring in this region leading to multiple necks. Further thinning of this region renders the drop highly unstable, resulting in the breakup. The neck disintegrates into fragments, forming satellites between the two main breakaway parts. It is interesting to compare the mechanism of drop breakup brought out by our simulation study with the experimental observations of Torza, Cox, and Mason, who studied the deformation and breakup of silicon oil drops in simple shear fields in various liquid media.25 Using a high-speed movie camera, they were able to capture the breakup process in a series of photographs, shown in Figure 9, for drops in saxtolphthalate, the system with the viscosity ratio, equal to 1.1, closest to that for our model systems. Remarkable similarity between the two sets of snapshots furnishes a clear validation of the DPD approach with respect to its application in the modeling of drop beakup in flow fields. Experimental studies of Grace29 as well as of Torza et al.25 show that the critical capillary number, Cac (the minimum value of Ca required to actuate the breakup process), as a function of viscosity ratio (ν) passes through a minimum occurring at ν ) 1. Their result at this value (25) Torza, S.; Cox, R. G.; Mason, S. G. J. Colloid Interface Sci. 1972, 38, 395.
of ν is 0.36 ( 0.05, which compares reasonably well with the value of 0.26 ( 0.05, found in our simulations. It may be mentioned that the value of the critical Ca at ν ) 1, obtained by Hinch and Acrivos30 from the O(Ca2) theory of Barthes-Biesel and Acrivos,24 is 0.25, very close to the simulation result. 5. Concluding Remarks The satisfactory agreement between the results of our simulations and the available theory and experiment on the pendant drop and the drop in a shear field constitutes adequate validation of the DPD approach for treating multiphase systems in external fields. Detailed mechanisms of the detachment and the breakup processes revealed by the present simulations are found to be entirely consistent with the experimental observations published in the literature. This, together with the previously published work,2 provides a convincing indication of the potential of the method for realistic modeling of processes such as multiphase flow in complex geometries, dropdrop and drop-particle interactions, and spreading of drops on rough solid boundaries under a variety of hydrodynamic and physicochemical conditions. The scope of DPD can be further enhanced by extending the algorithm to include flows other than simple shear. The correspondence between the DPD model and real systems may be established through the dimensionless numbers, Ca, Re, ν, B, and so forth. Depending on the problem, the DPD unit of time may be related to real time through phenomenological quantities such as ηR/Γ, FL2/ η, and L3η/kT, where L is a length characteristic of the
6350
Langmuir, Vol. 16, No. 15, 2000
Clark et al.
Figure 9. Sequence of snapshots of a silicon drop undergoing breakup in a simple shear field, obtained in the experimental study of Torza et al.25 Reprinted with permission from J. Colloid Interface Sci. 1972, 38 (2), 395.
problem. For the pendant drop and the drop in the shear flow, the appropriate quantities are respectively FL2/η and ηR/Γ (with L equal to R), giving values of 2.8 × 10-7and 6.3 × 10-7s for the DPD time units in the two systems for a 4.5 µm radius drop in water (η ) 10-3 Pa s; F ) 103 kg m-3), assuming Γ ) 10-2 N m-1 (oil/water). For equilibrium thermodynamic properties, the equation of state is reckoned to serve as a viable link between the model and the real systems.10 The breakup processes of >0.01 µm size droplets in lowviscosity media, for example, water, occur at Re > 1 (see earlier remarks in section 4.2.1). This is the hydrodynamic regime readily accessible to DPD simulations. The immediate relevance of such processes for the emulsion and colloid technologies would, therefore, make DPD a potentially important tool in the mesoscopic modeling of materials processing. The method does not, however, lend itself to simulations under low-Reynolds-number conditions, which remains its serious limitation. Simulation approaches such as smooth particle hydrodynamics (SPH), which are closely related to DPD,26 would prove more suitable for treating systems in that regime. Alternatively, the Lattice-Boltzmann (L-B) method may be considered, where appreciable progress has been made in the treatment of solid-fluid and fluid-fluid interfaces in moving systems.27
At the phenomenological level, the boundary integral approach has shown exciting promise for modeling complex fluid systems, albeit in the Stokes limit (Re , 1).28 It is hoped that the mesoscopic and phenomenological approaches will play complementary roles in the advancement of computer simulation methodologies for complex materials. LA991565F (26) Espanol, P. Europhys. Lett. 1997, 39, 605. (27) Cates, M. E.; Kendon, V. M.; Bladon, P.; Despalat, J.-C. Faraday Discuss. 1999, 112, 1. (b) Malevanets, A.; Yeomans, J. M. Faraday Discuss. 1999, 112, 237. (28) Zinchenko, A. Z.; Rother, M. A.; Davis, R. H. J. Fluid Mech. 1999, 391, 249. (29) Grace, H. P. Chem. Eng. Commun. 1982, 14, 225. (30) Hinch, E. J.; Acrivos, A. J. Fluid Mech. 1979, 91, 401. (31) Masters, A. J.; Warren, P. B. Europhys. Lett. 1999, 48, 1. (32) Realization of the Stokes regime (Re , 1) in the DPD simulations is inhibited by the dominance of fluctuations at low shear rates. An attempt to reduce the noise by lowering the temperature of the system would result in the occurrence of undesirable density oscillations at the fluid/fluid interface. One might consider increasing the viscosity of the system by increasing the value of γ, but orders of magnitude increases in the viscosity, required for accessing the Stokes regime, cannot be achieved in this way, since, say, a 5-fold increase in γ will lead to only a 2-fold increase in the viscosity.31 Moreover, the increase in γ will necessitate the reduction in the magnitude of the time step (γ δt < 0.3),31 meaning that at large γ the simulation runs would become prohibitively long for the requisite time spans to be covered.
