On Simulating Colloids by Dissipative Particle Dynamics: Issues and Complications

Journal of Colloid and Interface Science 242, 106–109 (2001) doi:10.1006/jcis.2001.7759, available online at http://www.idealibrary.com on

NOTE On Simulating Colloids by Dissipative Particle Dynamics: Issues and Complications We discuss the application of dissipative particle dynamics (DPD) to the study of colloidal systems composed of spherical particles. We develop a new model based on DPD particles with a central core and consider the effects of various parameters. The model is used to simulate (a) colloidal particle scattering as a means of probing the depletion interaction and (b) bulk dispersion viscosities for comparison with related models and to find evidence for aggregation in these systems. These preliminary simulations form a testing ground for assessing the feasibility of using DPD for studies of the rheology of particle gel structures. °C 2001 Academic Press Key Words: dissipative particle dynamics; mesoscopic methods; colloidal interaction; hydrodynamic interaction; computer simulation.

INTRODUCTION Brownian and Stokesian dynamic simulations have been used extensively to study colloidal systems modeled by spherical particles suspended in a hydrodynamic continuum (1, 2). Because of the long-range nature of the force, the hydrodynamic interaction between N such particles is a many-body problem with computation times scaling as N 3 and its strict inclusion for more than ∼100 particles is prohibitive. Many simulations of colloidal systems require larger scale models with N ∼ 103 to obtain satisfactory results. Our own experience with particle gels is typical (2, 3). The fractal structure of these gels is not apparent over short length scales, and so a meaningful study necessitates a system with at least ∼10 particles spanning in each direction. For this reason many studies, including our own, have neglected the full inclusion of the problematic hydrodynamic interaction. For this type of simulation, which we call simple Brownian Dynamics, only the first-order hydrodynamic interaction of the spheres with the continuum is considered in the so-called freedraining limit. Although the aggregation process and final structure of particle gels are not thought to be significantly influenced by hydrodynamic interactions, the rheology, particularly at high shear rates (or dimensionless Peclet number), almost certainly is. For this reason our own studies of rheology have been limited to low shear rates (4, 5). Even these simulations are limited by the method of including the shear field. For simple Brownian dynamics the homogeneous shear technique must be used, which imposes a velocity profile on the particles. This procedure takes no account of the heterogeneity of the material, which, in a real system, would set up corresponding inhomogeneities in the stress and strain fields. This could be overcome by applying the shear field using moving boundaries and allowing the internal stress and strain to respond naturally. Attempts to use such Lees–Edwards boundary conditions fail, however, for simple Brownian dynamics because the algorithm does not include a term to accelerate the suspending fluid in response to particle motion (5). Nevertheless, the rheology of particle gels, and colloidal systems in general, at high shear rates is of 0021-9797/01 $35.00

C 2001 by Academic Press Copyright ° All rights of reproduction in any form reserved.

great interest. A related problem, the study of fluid flow through a structured porous matrix of particles, perhaps a fractal gel, also requires proper treatment of hydrodynamics. Given that the exact inclusion of hydrodynamic interactions remains beyond our means at this time, one possible route for at least the approximate inclusion of hydrodynamic interactions in this type of simulation is through dissipative particle dynamics (DPD). In this method, rather than the fluid phase being treated as a continuum, it is simulated by particles that act as momentum carriers. The colloidal particles have been simulated either by rigid assemblies of these particles (6, 7) or by larger spheres (8, 9). Increasing the size of these “colloidal” particles compared with the fluid particles is expected to improve the quality of the hydrodynamic representation, especially at short range. However, this also increases the total number of particles in the simulation for a given number of “colloidal” particles and is limited by computer power. An additional complication arises from the depletion effects caused by structuring of the small particles between adjacent larger particles. The presence of this unavoidable effect of coarse-graining the hydrodynamic medium has been inferred in simulations using fused DPD particles to represent the colloidal particles (7). In this Note we examine some possible extensions of this DPD approach to colloidal systems and consider issues concerning the viability of the method for future work in this area. One purpose of the calculations is to explore to what extent depletion effects are manifest in the colloidal DPD simulation results.

SIMULATION MODEL The basic simulation method is derived from that described by Groot and Warren (10), which was later extended to a two-component mixture by Jones et al. (11). The brief description given here is designed to highlight the main differences of our proposed model to those workers. We envisage two types of generalized DPD momentum carriers or “particles” with “core” radii a1 and a2 and a DPD “coat” of thickness rc . The cutoff rc is set to unity and used as the length scale for the simulation. For two particles of types n, m with center– center distance r , a “surface-to-surface” separation h is then h mn = r − am − an .

[1]

The DPD interaction is defined in terms of this parameter rather than the usual center–center separation. Thus a weighting function ωi j is defined by ωi j = (1 − h i j /rc ),

[2]

106

where rc is a potential cutoff beyond which forces vanish. For the DPD fluid phase considered here the radius a2 was set to zero, although programming is easier to check if the general case is retained in principle. In this instance the type 2 particles revert to “normal” DPD particles and the weighting function is unity when the particle centers overlap. For DPD interactions between type 1 (colloidal) particles and type 2 (fluid) particles the weighting function is then unity when the type 2 center coincides with the type 1 nominal surface or when the surfaces of two type 1 “colloidal” particles touch. The weighting function

107

NOTE does not preclude the overlap of colloid particles (where the weighting function would exceed unity), but this is prevented by the colloidal interaction described below. Similarly a conservative interaction term, α12 , prevents significant penetration of the fluid particles below the colloidal particle surface. From this point our model follows Jones et al. (11) The total force acting on particle i is the sum of three contributions from the conservative FiCj , dissipative FiDj , and fluctuation FiRj forces: fi =

X j6=i

FiCj + FiDj + FiRj .

[3]

The DPD conservative interaction is usually assumed to be FiCj (h) = αωi j ei j ,

[4]

where ei j is a unit vector in the direction ri j . The fluctuation force is defined by β FiRj (h) = √ ζi j ωi j ei j , 1t

[5]

where 1t is the time step, β is a scaling parameter, and ζi j is a Gaussian random variable with zero mean and unit variance. The dissipative force scales with the parameter γ according to FiDj (h) = −γ ωi2j [ei j · (vi − v j )]ei j ,

[6]

where vi is the velocity of particle i. The fluctuation–dissipation theorem requires that β 2 /2γ = kB T,

[7]

where kB is the Boltzmann constant and T the temperature. p With kB T = 1 defining the unit of energy, the unit of time becomes t0 = mrc2 /kB T for particles of mass m. Following Groot and Warren (10) our test simulations have used a reduced time step of 0.05. N1 “colloidal” particles were set up with a central radius a1 and a DPD “coat” of thickness rc . The surrounding fluid was composed of N2 = (N − N1 ) DPD particles with a2 = 0. The mass, m 2 , of these particles was set to unity and defines the unit of time. The mass of the “colloidal” particles, m 1 , was adjusted to match the density of the “fluid.” To update the particle positions we used the modified velocity–Verlet algorithm described by Groot and Warren (10) as implemented by Jones et al. (11),

and Warren we performed most of our simulations at a DPD fluid density ρ2 = N2 /V2 = 3.0. The repulsion parameter α is related to the compressibility of the fluid and we choose α22 = 25kB T appropriate to water for this DPD number density. The noise level was normally set at β = 3.35 (11), which constrains the dissipative term to γ22 = 5.6 through Eq. [7]. A pure DPD fluid composed entirely of these type 2 particles equilibrated to give a radial distribution function equivalent to that shown by Groot and Warren (10). The “colloidal” particles interact with each other through the repulsive interaction α11 and the dissipative term γ11 . The latter would represent the energy dissipation mechanisms associated with deformable particles such as emulsion droplets. For the simulation of emulsion droplets these quantities would be expected to have values similar to those of the suspending fluid. In practice we have found that nonzero values of γ11 introduce an overshoot into the stress response to a step in strain rate. This represents a viscoelastic response excited by collisions of the colloidal particles. In our earlier simulations of particle gels such overshoots were associated with the time needed for restructuring the dispersion at the higher shear rate (4). In the present case, both of these parameters were set to zero since, in these simulations, we envisage the colloidal particles as relatively hard spheres with high compressibility. In place of α11 we introduce a more realistic colloidal interaction acting only between particles of type 1 of the form φ(r ) = E 0 exp(−κh);

F(r ) = F0 exp(−κh);

h < rc .

[9]

Here h is the surface–surface separation, φ is the interaction potential, F0 = E 0 /κ is the force magnitude and κ is a range parameter. The potentials used are sketched in Fig. 1 and could represent an electrostatic or steric repulsion. Colloidal interactions are frequently much shorter ranged than the potentials used here but the use of more rapidly changing potentials may necessitate a shorter time step and longer computation times despite the larger mass of the colloidal particles. The remaining terms α12 and γ12 represent interactions and dissipation between the fluid phase and the particles. The effect of γ12 was not clearly revealed in our results, but a nonzero value would represent interfacial dissipation in the simulation of particles with a surface polymer coat. The interaction between fluid and colloidal particle represented by α12 was found to affect the height of the first peak in the radial distribution function defined between fluid and colloidal particle. A value of α12 = 12.5 reduced this almost to 1.0 and so was adopted. The resulting model is very similar to that recently presented by Dzwinel and Yuen (12). With all DPD parameters set to zero for the colloidal particles the dissipative terms in the algorithm vanish and for these particles the update is identical to a normal molecular dynamics Verlet scheme.

µ · ¶¸ 1 Fi (t) = f ri (t), vi t − 1t , 2 µ ¶ 1 1 vi (t) = vi t − 1t + Fi (t)1t, 2 2 ri (t) = ri (t − 1t) + vi (t)1t + µ ¶ 1 1 vi t + 1t = vi (t) + Fi (t)1t, 2 2

1 Fi (t)1t 2 , 2

[8]

where the function f gives the total force on the particle i at time t according to Eq. [3] and includes the dependence on relative particle velocities. A shear field can be introduced through the Lees–Edwards moving boundary condition (15). In this scheme the periodic images in the z–direction are moved with velocity ± G L/2, where G is the required shear rate and L the box length. The resulting stress propagates through the DPD fluid, setting up a shear field. In the early stages this generates a typical sigmoid velocity profile which eventually settles down to a mean linear steady state. Unlike the homogeneous shear method this shear field is responsive to local inhomogeneities in the system (5). For systems including two types of particle, the quantities α and γ each have three independent values α11 , α12 , α22 and γ11 , γ12 , γ22 (11). Following Groot

FIG. 1. Simulations of colloidal particle scattering trajectories using DPD. The inset shows the interaction force F(h) used. F0 /mrc t0−2 = 75.0, G/t0 = 0.1. κrc = 5.0(s); 7.5; (m) 10.0; (n) 19.0. (· · ·). Classical Hydrodynamics (—).

108

NOTE

RESULTS

Colloidal Scattering Experiments We have performed two types of simulation to test this model. The first involves only two colloidal particles and is comparable to a colloidal particle scattering (CPS) experiment (13). In such an experiment a particle is fixed to a wall, and the motion of a second, mobile particle is tracked as it passes near the fixed particle in a shear field. For a purely hydrodynamic interaction the trajectories are perfectly symmetrical about the central particle (13). However, the symmetry is broken by the presence of a colloidal interaction. For our present purposes the presence of the wall is an unnecessary complication. Instead we arrange for one “colloidal” particle to be immobilized at the origin while a boundary-driven shear field drives the motion of the second. It is possible to do this under conditions of kB T = 1, but Brownian motion destroys the detail in the results. In our studies of CPS using a standard deterministic hydrodynamic model (13, 14) Brownian motion is absent unless deliberately introduced. For comparative purposes only, it is possible also to set kB T = 0 in this DPD simulation. In this case we have β22 = 0, and so there is no random motion, but the DPD particles still act to dissipate the shear field since γ22 remains nonzero. We retain the same energy and force scale that we would have with kB T = 1, but this has to be reexpressed in terms of the time t0 . Figure 1 shows some results for this experiment, including a repulsive colloidal interaction with four different values of the range parameter. The purely hydrodynamic trajectory obtained by the methods of Ref. (13) is also included for comparison. Since we have set α11 = 0, this colloidal interaction is entirely responsible for the excluded volume of the particles. The magnitude of the interaction potential, represented by F0 , is just sufficient to lead to repulsion, but weaker forces will allow the colloidal particles to overlap and coalesce. For this set of simulations the colloidal particle size was set at a1 = 2rc and the mobile particle was made to move in a fluid of N2 = 5998 DPD particles (with a2 = 0) by applying the shear field. The longest ranged interaction (κrc = 5.0) leads to a repulsive-type collision of the mobile particle characterized by a positive deflection. The two intermediate range values (κrc = 7.5, κrc = 10.0) give almost equivalent trajectories which are more nearly symmetrical, such as would normally be characteristic of a purely hydrodynamic collision. However, the shortest ranged interaction (κrc = 19.0) leads to capture of the mobile particle in the stable downstream position, as we would expect for an appreciably attractive colloidal force (15). This behavior is also found for even shorter ranged interactions examined but not shown in the figure. There is not a smooth transition to this behavior as the range parameter is increased; instead it appears to switch suddenly from a near-symmetrical response. Responses (10.0 < κrc < 19.0) are indeterminate in that they can give either behavior. Since there is no introduced noise in this system, we surmise that the behavior in this region ultimately depends on the exact arrangements of DPD particles throughout the simulation. These results are further confirmatory evidence for the presence of a depletion interaction in these systems. Boek and van der Schoot (6) first described this effect for their simulations of colloids composed of “frozen” DPD particles. It is also notable that, compared with the purely hydrodynamic path, the initial phase of all the trajectories is relatively flat and even moves toward slightly lower values of y as the fixed particle is approached. In this region, the mobile particle is beyond the influence of the colloidal interaction, which explains why the initial trajectories for all four interactions are essentially the same. For this reason we would expect purely hydrodynamic behavior. The observed negative deviation from this expectation is again likely to be due to the depletion effect, which is the same for all cases.

Viscosities Boek and van der Schoot (7) have presented viscosities for DPD systems in which the colloidal particles were represented by sets of “frozen” DPD particles. In order to compare with their results we choose to work with a volume fraction of big particles, φ, close to 0.3. These simulations used a total of N = 6 × 104 particles with the colloidal particle size set to a1 = 3.5rc . For these parameters we have N1 = 46, m 1 = 538.8 m 2 and φ = 0.293. Viscosities were obtained by computing the stress in response to an imposed strain over runs of 104 time

FIG. 2. Stress response Sxy to a step in strain rate to G = 0.05 with and without colloidal particle dissipative terms. (s) γ11 = γ12 = γ22 = 5.6, α11 = 0.0, α12 = α22 = 25.0; (d) γ11 = γ12 = 0.0, γ22 = 5.6, α11 = 0.0, α12 = 12.5, α22 = 25.0. steps. The stress is obtained from the usual interparticle average hSαβ i =

N −1 X 1 NX (rαi j rβi j /ri j )ei j · fi j V i=1 j>1

[10]

In Fig. 2 the stress responses are shown for systems with and without the dissipative terms, γ11 and γ12 , for the colloidal particles themselves and their interaction with the DPD “fluid.” As we have already mentioned, the inclusion of these dissipative terms leads to a significant stress overshoot that mainly originates with γ11 . For this reason γ11 is set to zero. We have performed simulations with and without γ12 but find its effect more ambiguous. Without these terms the stress tends to show some early relaxation toward dynamic equilibrium and viscosities are therefore obtained from an average taken beyond any such initial response. In order to obtain relative viscosities, the viscosity η0 of the singlephase fluid was first obtained by the same methods. The relative viscosities ηr obtained for the dispersed systems are reported in Tables 1 and 2 for several choices of parameters. Shear rates for the dispersed phases are given in terms of the dimensionless Peclet number Pe defined by Pe =

6π η0 a 3 G kB T

[11]

The values shown in Table 1 are obtained for F0 /mrc t0−2 = 75.0, and κrc = 5.0—the same interaction that was found to give the best symmetry in the colloidal scattering results. However, these results are higher by some 20% than those reported by Boek and van der Schoot (7) for these conditions. By

TABLE 1 Relative Viscosities ηr for Different Shear Rates Defined by Peclet Number Pe Pe

γ12

ηr

g11 (r )max

0.0 1.83 5.12 10.3 18.3 18.3

5.6 5.6 5.6 5.6 5.6 0.0

0.0 4.52 ± 2.6 5.11 ± 2.0 5.06 ± 1.5 4.12 ± 0.83 3.36 ± 0.72

22.0 21.7

Note: φ = 0.293, α11 = 0.0, α12 = 12.5, α22 = 25.0, γ22 = 5.6, F0 / mrc t0−2 = 75.0, κrc = 5.0.

109

NOTE

TABLE 2 Relative Viscosities ηr for Different Shear Rates Defined by Peclet Number Pe Pe

ηr

g11 (r )max

0.0 1.83 5.12 10.3 18.3

0.0 1.72 ± 1.0 2.40 ± 0.95 2.30 ± 0.54 2.69 ± 0.60

10.0 8.98 8.65 6.35 4.30

Note: φ = 0.293, α11 = 0.0, α12 = 12.5, α22 = 25.0, γ22 = 5.6, F0 / mrc t0−2 = 200.0, κrc = 5.0.

setting F0 /mrc t0−2 to 200 and considerably increasing the repulsive strength of the interaction, viscosities roughly 50% of those reported by Boek and van der Schoot (7) are obtained. For some of these results the radial distribution function for the large particles, g11 (r ), has been computed. Since few particles are available this was averaged over 104 steps to obtain good statistics. The zero-shear peak values for both interaction strengths are much higher than the expected value of ca. 2.5 for a 30% suspension of hard spheres. The results in Table 2 also show that the g11 (r ) peak height reduces with increasing shear rate. These observations can be explained on the basis of aggregation taking place. Notably, the peak in g11 (r ) is highest for the weaker interaction (Table 1). This is consistent with the presence of a strong underlying depletion interaction due to the coarse-grained nature of the suspending fluid.

Computational Timings and Viability Using linked neighbor lists the simulation ran at ∼0.54 s/time step for a total of 6000 particles and ∼6 s/time step for 6 × 104 particles on a single Silicon Graphics Origin2000 R10000 processor. The scaling with system size is thus roughly linear as we would expect for the link cell method (15). The volume fraction most commonly used for our Brownian dynamics gel simulations was 0.05 and the closest comparable runs using this DPD model were for 11 colloidal particles with radius 3.5 and a volume fraction 0.0899. Because the interesting behavior of gels is linked to their long-range structure, we require at least 1000 colloidal particles to perform a meaningful simulation. For these volume fractions we would therefore require >5 × 106 particles and would expect computation times of ∼500 s/step. Since a minimum of 104 steps would normally be required in a simulation run this is almost certainly prohibitive and parallel code would be needed on current machines to even approach viability. On this point we are therefore in agreement with Tanaka and Araki (16) and we believe that numerically more efficient methods, such as the one they describe, must be sought for this problem.

significant aggregation. Thus, the CPS experiment for two particles indicates that the combined effect of this interaction and the depletion effect is neutral or slightly repulsive whereas the dispersed system displays evidence that the combined forces are attractive. Taken together, these results may therefore be some evidence of cooperative effects in the depletion interaction. DPD simulations have frequently used sets of “frozen” DPD particles to simulate a dispersed phase. This has the advantage of being able to model particles of any shape, but it incurs some computational overhead. Also the resulting particles retain interactive dissipative terms and a colloid–colloid potential arising from a sum over the DPD conservative interactions. The resulting interaction potential is not necessarily appropriate to colloids and the construction of a realistic potential for these systems must compensate accordingly. The scheme presented here does not suffer from this complication. Our results further suggest that the inclusion of dissipative terms for the colloidal particles may not always be desirable since they lead to viscoelastic response.

ACKNOWLEDGMENTS E. D. acknowledges financial support from the Biotechnology and Biological Sciences Research Council (UK). M.W. thanks E. S. Boek for a stimulating discussion.

REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.

Heyes, D. M., and Mitchell, J., J. Chem Soc. Faraday Trans. 90, 1931 (1994). Whittle, M., and Dickinson, E., Mol. Phys. 90, 739 (1997). Dickinson, E., J. Colloid Interface Sci. 225, 2 (2000). Whittle, M., and Dickinson, E., J. Chem. Phys. 107, 10191 (1997). Whittle, M., and Dickinson, E., J. Chem. Soc., Faraday Trans. 94, 2453 (1998). Koelman, J. M. V. A., and Hoogerbrugge, P. J., Europhys. Lett. 21, 363 (1993). Boek, E. S., and van der Schoot, P., Int. J. Mod. Phys. C 9, 1307 (1998). Boek, E. S., Coveney, P. V., Lekkerkerker, H. N. W., and van der Schoot, P., Phys. Rev. E. 55, 3124, (1997). Flekkoy, E. G., and Coveney, P. V., Phys. Rev. Lett. 83, 1775 (1999). Groot, R. D., and Warren, P. B., J. Chem. Phys. 107, 4423 (1997). Jones, J. L., Lal, M., Ruddock, J. N., and Spenley, N. A., Faraday Discuss. 112, 129 (1999). Dzwinel, W., and Yuen, D. A., J. Colloid Interface Sci. 225, 179 (2000). Whittle, M., Murray, B. S., Dickinson, E., and Pinfield, V. J., J. Colloid Interface Sci. 223, 273 (2000). Whittle, M., Murray, B. S., and Dickinson, E., J. Colloid Interface Sci. 225, 367 (2000). Allen, M. P., and Tildesley, D. J., “Computer Simulation of Liquids.” Oxford University Press, Oxford, 1987. Tanaka, H., and Araki, T., Phys. Rev. Lett. 85, 1338 (2000).

CONCLUSIONS The DPD simulations that we describe here are a generalized extension of earlier approaches that could be used to model colloidal and emulsion systems. They have much in common with the recent simulations reported by Dzwinzel and Yuen (12). However, these authors studied aggregation using an attractive Lennard–Jones interaction, and this may have obscured the effects of the depletion interaction which is responsible for aggregation here. The simulations of colloidal particle scattering experiments show primary evidence of a strong effective attractive interaction between the particles that can be attributed to the coarse-grained DPD fluid. Indirect evidence is also found in the measurements of relative viscosity and the radial distribution function for the colloidal particles in bulk dispersed systems. Interestingly, the same interaction that produced a nearly symmetrical trajectory in the CPS experiment (F0 /mrc t0−2 = 75.0, κrc ∼ 5.0) gave a system with a very strong peak in the radial distribution function indicating

Martin Whittle1 Eric Dickinson The Procter Department of Food Science University of Leeds Leeds LS2 9JT, UK Received January 30, 2001; accepted June 1, 2001; published online August 16, 2001

1 To whom correspondence should be addressed at present address: Department of Information Studies, University of Sheffield, Western Bank, Sheffield S10 2TN, UK.