Mesoscopic Dynamics of Colloids Simulated with Dissipative Particle Dynamics and Fluid Particle Model

Witold Dzwinel¹, ÝDavid A.Yuen^2* Ýand Krzysztof Boryczko¹

¹AGH University of Mining and Metallurgy, Institute of Computer Science, KrakÛw, Poland,

Ý²Department of Geology and Geophysics and Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55415-1227, USA,

J.Molecular Modeling, submitted for publication, July 2001

Abstract

In this electronic communication we report results taken from numerical simulations of complex fluids, using a combination of discrete-particle methods. Our molecular modeling reportoire comprises three simulation techniques:Ý molecular dynamics (MD), dissipative particle dynamics (DPD) and fluid particle model (FPM).Ý This type of model can depict multi-resolution molecular structures found in complex fluids ranging from single micelle, colloidal crystals, large-scale colloidal aggregates up to the mesoscale processes of hydrodynamical instabilities in the bulk of colloidal suspensions. We can simulate different colloidal structures in which the colloidal beds are of comparable size to the fluid particles. This undertaking is accomplished with a two-level discrete particle model consisting of the MD paradigm with a Lennard-Jones (L-J) type potential for defining the colloidal particle system and DPD or FPM for modeling the solvent. We observe the spontaneous emergence of spherical or rod-like micelles and their crystallization in stable hexagonal or worm-like structures, respectively. The ordered arrays obtained by using the particle model are similar to the 2-D colloidal crystals observed in laboratory experiments. The micelle shape and its hydrophobic or hydrophilic character depend on the ratio between the scaling factors of the interactions between colloid-colloid to colloid-solvent. Unlike the miscellar arrays, the colloidal aggregates involve the colloid-solvent interactions prescribed by the DPD forces. Different from the assumption of equilibrium growth, the two-level particle model can display much more realistic molecular physics, which allow for the simulation of aggregation for various types of colloids and solvent liquids over a very broad range of conditions. We discuss the potential prospects of combining MD, DPD and FPM techniques in a single three-level model. Finally, we present results from large-scale simulation of the Rayleigh-Taylor instability and dispersion of colloidal slab in 2-D and 3-D.

Keywords: complex molecular fluids - fluid particle model - dissipative particle dynamics - molecular dynamics simulations - cross-scale simulations - phase separation

^* To whom correspondence should be addressed

1. Introduction

2. Fluid Particle Model

3. Three-level model

4. Colloidal arrays

5. Multiresolution molecular structures developed in mixing induced by the Rayleigh-Taylor instability

6. Concluding Points

7. References

Introduction

Today there exist multitudinous numerical methods for modelling physical and chemically reacting phenomena in complex molecular fluids. They include molecular dynamics, lattice Boltzmann gas and cellular automata [Rothman and Zaleski, 1997; Chopard and Droz, 1998,], diffusion limited aggregation [Kremer, 2000, Meakin, 1998] or those employing finite element simulation, e.g., for Cahn-Hillard fluids [Cahn and Hillard, 1958]. Our model falls under the category of the discrete-particle paradigm and comprises three distinct kinds of numerical techniques: molecular dynamics (MD), dissipative particle dynamics (DPD) and fluid particle model (FPM) techniques. The dissipative particle dynamics method [Hooggerbruge and Koelman, 1992] and its most recent generalization ñ the fluid particle model (FPM) [EspaÒol, 1998] ñ can provide a bridge between the basic domains of complex fluids, i.e., micelles and colloidal particles and large-scale structures such as droplets, crystals, colloidal agglomerates and self-organized patterns from molecular fluid instabilities.

The distinct advantage held by DPD over the other methods is the direct possibility for matching the scales of the discrete-particle simulation to the dominant spatio-temporal scales of the entire system. The real cross-scale endeavor, where different kinds of timesteps can be used, reduces then to a proper definition of DPD interactions by using bottom-up [Flekk¯y and Coveney, 1999] or top-down approaches [EspaÒol et. al, 1999]. We may then work on clusters of atoms. As shown in Fig.1, the relatively narrow gap between the smallest structures and large-scale structures in complex fluids can be connected by using reasonable number of particles (a hundred thousand in 2D and few millions in 3-D).

One disadvantage of DPD is the lack of a drag between central particle and the second one orbiting in a circumference around the first one. The non-central force, introduced in fluid particle model (FPM), is proportional to the difference between particle velocities and eliminates this defect. However, this makes the model more complex and suggests its use for longer length-scales, where the FPM particle is large enough in that it interacts only with its closest neighbors. The other problem with DPD is that DPD particles cannot simulate "solid" granular material with attractive forces. The colloidal particles are just such the solid grains immersed in a solvent. For simulating suspensions consisting of small number of colloidal beds in which the hydrodynamic interactions between colloidal particles play the dominant role, we can employ the fluid particle dynamics (FPD) approach [Tanaka and Araki, 2000]. In FPD, the colloidal beds are MD particles interacting with the Lennard-Jones forces. The hydrodynamic interactions are derived from the global velocity field, which is computed directly by integrating the Navier-Stokes fluid dynamical equations. For systems with colloidal beds of a similar size as the complex fluid microstructures (e.g., polymeric clusters, large blood cells) or not much larger, the bed-solvent particle interactions become important. These molecular interactions are responsible for creating colloidal microstructures, such as micelles and colloidal crystals, which cannot be simulated within the framework of FPD.

In [Dzwinel and Yuen, 2000, 2] we attack this problem by employing a two-level particle model in which the colloidal beds are simulated with MD particles, while the solvent is represented by the DPD model. Because both the colloidal beds and fluid particles areÝ comparable in size, we have employed a uniform time-step in the integration of the Newtonian equation of motion.

Figure 1 ÝMultiresolution in a complex fluid.

We will demonstrate the cross-scaling results by simulating complex fluids with a three-level method for which colloidal beds are represented by MD particles but fluid particles - simulated with DPD and/or FPM particles ñ represent a solvent. At first, the three-level method is presented starting from the most general fluid particle model. By introducing two-level model comprising MD and DPD (or FPM) forces, we can simulate miscellar solutions, which undergoes crystallization into stable hexagonal or worm-like structures. By merely changing the character of the particle-particle interactions, we can go up the spatio-temporal scale simulating large colloidal aggregates, which involve as many as 20 million particles. On the largest scale, we focus our attention on the Rayleigh-Taylor [Taylor, 1950] instability, which induces mixing of the colloidal mixtures and dispersion of colloidal slab. These simulations are carried out with a three-level method (MD, DPD and FPM) and with a combined total number of particles in the neighborhood of several million. Finally, we summarize the results and discuss the prospects of the discrete particle method as a computational tool for modeling mesoscopic dynamical phenomena of complex fluid.

Fluid Particle Model

A theoretical framework for fluid particle model can be found in [EspaÒol, 1998]. The fluid particles defined by its mass m_i, position r_i, and velocity v_i interact with each other. The particle can be viewed as a ìdropletî consisting of liquid molecules with an internal structure and with some internal degrees of freedom. We use the two-body, short-ranged DPD force F as postulated in [EspaÒol, 1998]. This type of interaction consists of a conservative force F^C, two dissipative components F^T and F^R and a Brownian force, which are defined by:

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ ÝÝ(1)

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ Ý(2)

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ(3)

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (4)

Ý (5)

where:

r_ij ñ is a distance between particles i and j, r_ij= r_i ñ r_j is a vector pointing from particle i to particle j and e_ij=r_ij/r_ij,

D ñ is spatial dimension,

T ñ is the temperature of particle system,

k_B ñ the Boltzmann constant,

dt ñ is the timestep,

g - scaling factor for dissipation forces,

w - angular velocity,

dW^S, dW^A, tr[dW]1Ý - are symmetric, antisymmeric and trace diagonal random matrices of independent Wiener increments defined in [EspaÒol, 1998].

A(r), B(r), C(r), , F(r) ñ functions dependent on the interparticle separation distance r=r_ij.

T ñ is a dimensionless matrix given by:

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (6)

1 ñ is the unit matrix.

As proved in [EspaÒol, 1998], the one-component FPM system yields the Gibbsí distribution as the equilibrium solution to the Fokker-Planck equation under the detailed balance (DB) Ansatz. Consequently it satisfies the fluctuation-dissipation theorem. According to the fluctuation-dissipation theorem the normalized weight functions are chosen such that:

ÝÝÝÝ (7)

The non-central force in FPM, which is proportional to the difference between particle velocities, introduces an additional drag lacking in the DPD model. The central particle and the second one orbiting about around the first one, do not produce any force in DPD algorithm, thus making it inconsistent. The non-central force results also in an additional rotational friction given by Eq.(4). The temporal evolution of the particle ensemble obeys the equations:

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (8)

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ Ý ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (9)

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ(10)

where the torques in Eq.(11) are given by:

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (11)

ÝOne can verify easily that the total angular momentum is conserved.

The main purpose of this model is to generalize both the smoothed particle hydrodynamics method (SPH) [Libersky et. al, 1993] ñ the particle based algorithm used for simulations in macroscale - and dissipative particle dynamics. FPM can predict precisely the transport properties of the fluid, thus allowing one to adjust the model parameters according to the formulas of kinetic theory. Unlike the SPH, the angular momentum is conserved exactly in FPM. The FPM model can thus be interpreted as a Lagrangian prescription of the non-linear fluctuating hydrodynamic equations [EspaÒol, 1998].

Three-level model

The FPM model represents a generalization not only of DPD but also of the molecular dynamics (MD) technique. It can be used as DPD by setting the non-central forces to zero or MD, by dropping the dissipative and Brownian components. The fluid particle model holds an advantage over DPD but only for larger scales where the fluid particles are adequately large and can interact only with their closest neighbors. In this situation DPD is less efficient because many more particles than for FPM should be involved for creating a drag between circumvented DPD particles. DPD is computationally more efficient than FPM at smaller scales, for which the interaction range (r_cut) of the potential must be longer

The three-level numerical model is designed for simulating complex fluids. Three types of particles are defined accordingly by:

Colloidal particles (CP), with an interaction range „2.5¥l, where l is a characteristic length, equal to the average distance between particles. The CP-CP interactions can be simulated by a soft-sphere, energy-conserving potential with an attractive tail. The CP-CP forces conservative in nature and are given by Eq.(1).

Dissipative particles, the ìdroplets of fluidî represented by solvent particles (SP) located in the closest neighborhood of the colloidal particles with an interaction range „2.5¥l. The SP-SP and CP-SP forces represent only the two-body central forces given by Eqs.(1-5).Ý

Fluid particles (FP), the ìlumps of fluidî represented by the particles in the bulk solvent, with an interaction range £1.5¥l. Non-central forces are included within this framework.

In conventional DLVO theory (Derjaguin, Landau, Verwey, Overbeek) [Grier and Behrens, 2001] the long-range electrostatic interactions between colloidal spheres can be modelled by a screened-Coulomb repulsion [Daniel and Audebert, 1999]. However, some experimental findings [Larsen and Grier, 1997, Grier and Behrens, 2001] show that like-charged macroions have been attracted to one another by short ranged forces. This phenomenon cannot be explained by conventional theories. The recent simulation results [Prausnitz and Wu, 2000] show that the fluctuations of the charge distribution by the small ions result in the attraction between microions. The mean force is a combination of hard-sphere and electrostatic force. As shown in Fig.2, the Lennard-Jones (L-J) force [Haile, 1992] is close to the mean force obtained from large-scale Monte-Carlo calculations performed for a real colloidal mixture [Prausnitz and Wu, 2000]. The experimentally measured potential [Larsen and Grier, 1997] between a pair of 0.65m-diameter polystyrene sulphate spheres, which lie very close to a electrically charged glass, is also very similar to the L-J interactions. A better approximation of the colloid-colloid interactions is possible by adding to the L-J force a very steep force with a soft core (see Fig.2). However, because of simplicity, we will use the Lennard-Jones force as a sufficiently accurate approximation for the effective force F(r_ij) (see Eq.(2)).

Figure 2 The model of the two-body conservative forces representing the CP-CP interactions.

Solvent particles are represented by DPD particles and they are employed for simulations involving a larger cut-off radius. In the case of three-level computations in which the interaction range of CP particles is „2.5¥l, DPD particles are defined only within the nearest neighborhood of CP particles. We employ the FPM model for simulating the rest of the fluid system. Due to the comparable size for the free types of particles, the timestep for integrating the Newtonian equations of motion is uniform.

Thus the three-level system consists of three different procedures, each representing a particular technique of interparticle interactions.Ý The main assumptions in this model are as follows:

1) We consider here a two dimensional isothermal system, which consists of M particles. The number of species can easily be generalized for simulating multiphase systems [Coveney and EspaÒol, 1997],

2) The particle system is simulated inside a rectangular box with periodic or hard wall boundary conditions. Particles of various types can be scattered randomly in the box, i.e., this multi-component system can be perfectly mixed initially, or separated by a sharp interface (stratified, circle, rectangular, random shape).

3) The weight functions (Eq.(6)) satisfy the conditions imposed. Due to the freedom allowed by the model in selecting the weight functions, we have assumed

ÝÝ (12)

where, r_cut ñ is a cut-off radius, which defines the range ofÝ interactions. For r_ij>r_cut, F_ij=0.

The equations of motion are integrated by using fourth-order time-stepping O(Dt⁴) scheme [Boryczko et. al., 2001]

From the continuum limit equation Eqs.(14) :

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (13)

[Hoogerbrugge, Koelman, 1992], we compute the value ofÝ P - the scaling factor for the conservative interactions (Eq.12). The g parameter of dissipative forces are adjusted according to the formula from kinetic theory [EspaÒol, 1998]:

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (15)

where: n_b ñ is a bulk kinematic viscosity, D ñ is dimensionality, c²= k_BT/m,

ÝÝÝÝÝÝÝÝÝ (16)

Colloidal arrays

With the two-level model (DPD and MD) we can study the molecular dynamics of binary fluids. This system consists of colloidal particles (CPs) immersed in a solvent represented by an ensemble of dissipative, solvent particles (SPs). The colloidal particles are scattered randomly in the box, i.e., this binary system is initially well mixed.

The collision operator for the whole particle system consists of DPD and Lennard-Jones components is given by the formula:

A(t)=0,Ý W_ij^t = 0Ý for r_ij>r_cut,ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (17)

Thus, the CP particles interact with the Lennard-Jones potential with both CPs and SPs particles. The dual character of SP particles interactions ñ dissipative particle dynamics with other SP beads and molecular dynamics with L-J interaction and CP particles ñ can be explained by assuming that the SP ensemble can simulate a complex fluid. For example, the SP particles can represent parts of polymeric chains whose mutual interactions are modeled by using DPD forces, while the conservative and repulsive-attractive forces are responsible for gluing the chains to the hard-core nano-particles. The fundamental parameters of the particle system are presented in Table 1 in dimensionless units.Ý The size of the system is scaled-up to the cut-off radius. The timestep and the mass of colloidal particles are set to 1. The energy unit (e_Unit) is set arbitrarily as a reference point for scaling both the L-J CP-CP and CP-DP interactions. It coresponds to the e parameter of the Lennard-Jones force.

ÝÝÝÝÝÝÝÝÝÝÝ In [Dzwinel and Yuen, 2000, 2] we show that, depending on e_CP-CP/e_CP-SP aspect ratio, different colloidal structures (2-D colloidal arrays) can be simulated. Simultaneously we have repeated the calculations by using FPM with r_cut=2.5¥l instead of DPD (for FPM with r_cut=1.5¥l we do not observe the creation of micelles). The simulations were carried out at different concentrations of CP particles, ranging from 10% to 30%.

Table 1 Principal parameters of the particle system

Parameter

Value (in dimensionless units)

Entire particle

system

r_c ñ unit of length

Dt ñ unit of time

M=M_CP=M_SP ñ unit of mass

n ñ particle density

system size in r_c Ýunits

e_Unit ñ unit of energy

T₀ ñ temperature of the system

6.37

0.4

4.75¥10^-5

0.1¥e_Unit

DPD fluid

P - scaling factor for conservative

Ýforces

10, 100

3.8¥10^-3

CP-CP interactions

e_CP-CPñ L-J well depth

s_CP-CPñ L-J parameter

0.2-1.6¥e_Unit

0.4

CP-SP interactions

e_CP-SPñ L-J well depth

s_CP-SPñ L-J parameter

0.1-1.6¥e_Unit

0.4

The results shown below are produced by assuming a concentration of CP particles of 20%.Ý For a low aspect-ratio configuration we find that the ìhydrophilicî circular micelles can self-organize and be formed spontaneously. The two-dimensional hexagonal colloidal arrays (see Fig.3) are produced as a consequence. These structures are very common [Taupin, 1995] and are observed in laboratory experiments, e.g., see [http://www.lrsm.upenn.edu/~dinsmore/ surfxtal.html, Dinsmore, Yodh, 2000]. By increasing the ratio (see Fig.4 and movie), the hexagonal phase undergoes a self-organized transition to a lamellar metastable phase [Taupin, 1995, Poon, 2000]. The rod-like micelles form worm-like arrays, as shown on the URL http://chemeng.stanford.edu/~gastgrp/images/dendritic-Xsmall.jpg]. The lamellar metastable phase evolves into ìhydrophobicî hexagonal arrays, as shown in Fig.5, which then produce larger aggregates.

The clustering of micelles and small nano-particles forming colloidal aggregates plays a very prominent role in the development of new materials with a scale size ranging from nano to mesoscale. Fractal aggregates represent very fragile mechanical structures, which can be easily torn apart as a result of adequately strong external forces. Therefore, aggregating additives are used for controlling the rheology of paints and other coating systems. At low shear rates these shear-thinning non-Newtonian fluids have a high viscosity and a low viscosity at high shear rates.

Figure 3 Colloidal structures of circular îhydrophilicî micelles simulated with two-level model. An hexagonal order can be observed. (e_CP-CP/e_CP-SP =1, M=5¥10³, number of timesteps N=10⁶). See the supplementary materials for zoom-in.

ÝÝÝ

Figure 4 Worm-like structures in colloidal suspension. On the right, the zoom-in of the lamellar and coexistence of the ìhydrophobicî phases. The ratio e_CP-CP/e_CP-SP = 5, the number of particles is M=2¥10⁴, the number of timesteps is N=10⁶_. See the supplementary materials for movie.

There are many numerical techniques devised for simulating the colloidal aggregates. Most of these methods employ the Smoluchowski principle of coagulation according to a given reaction scheme [Meakin, 1988, Meakin, 1998]. These methods are still far from achieving reality. They are adequate for investigating static fractal structures of large agglomerate by assuming a low initial concentration of colloidal particles. For denser systems, the rheological properties of solvent and the mechanisms of aggregation vary tremendously with the particle concentration, which cannot be predicted with simple composite theories [Larson, 1999].

Figure 5 Colloidal structures of circular îhydrophobicî micelles simulated with two-level model. When we zoom-in we can observe the hexagonal order (e_CP-CP/e_CP-SP =16, M=2¥10⁴, number of timesteps N=10⁶). See the supplementary materials for zoom-in.

As shown in Fig.5, the two-level method fits in very well in fine-grained structures during the initial aggregation. However, the method is too computationally demanding for investigating large structures. More reliable statistics are required for monitoring the global properties of the growth process. We propose an alternative approach for simulating large aggregates, which will entail modifying the collision operator given in Eq.(17).

Going up to a larger spatio-temporal scale, we can assume CP particles can be represented not by hard-core beads but by the micelles. The SP particles are the DPD ìdropletsî of a complex liquid solvent of the size of CPs. At this time the CP particles interactions assume a dual character.Ý In order to avoid fluidization of the colloidal particle system and allow them to aggregate, we have insisted that the colloid-colloid interactions should possess hard-sphere core with a very short-ranged adhesive part [Gast and Russel,1998]. The soluble additives are excluded from a region, which is comparable to their own size near the particle surface, thus producing what is commonly called a depletion force [Yaman et.al, 1998], which is of entropic origin. For the sake of simplicity, we have assumed here that CP-CP interactions can be modeled by a Lennard-Jones 6-12 potential, as it was in the colloidal crystals case. In order to tune the system better we can employ more realistic depletion potentials [Daniel and Audebert, 1999] or tabularized experimentally measured potentials [Larsen and Grier, 1997]. Unlike the case of miscellar arrays, the CP-SP interactions are simulated by employing DPD forces.Ý This can be justified on the grounds that now the SP beads represent ìdropletsî of complex fluid, but not the fragments of polymeric chain. Because the CP ìdropletî contains a hard core, we have modified the repulsive portion of the conservative F^C CP-SP forces. This makes the potential to be steeper than that for the SP-SP interactions.

Figure 6 Colloidal aggregates simulated with two-level model (M=1.2¥10⁶, number of timesteps N=10⁵_,33% colloidal particles 67% fluid). See the supplementary materials for zoom-in.

In Fig.6 we present a snapshot from MD-DPD simulation of a developed colloidal aggregate. By zooming into a particular region of this frame, we can see clearly the hexagonal structure of the aggregated body.

Employing the two-level model, we have studied the scaling properties of mean cluster size S(t), expressed in number of particles versus time, by assuming a high concentration of colloidal particles in the system. For the cases of non-cohesive systems, with a low concentration of colloidal beads, the asymptotic growth for tÆ€ of the mean cluster size S(t) is given by:

ÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝÝ (18)

where k is the scaling-law index.

We show that in dissipative solvent of high concentration of colloidal particles, the growth of mean cluster size can be described by the power law S(t)µt^k (see Fig.7a). We have found the intermediate DLA regime, for which k=‡. It spans for a relatively long time. The intermediate regime depends on physical properties of the solvent such as the viscosity, temperature and partial pressure. The character of clusters growth varies with time and the exponent k shifts at longer timescales from ‡ to ª1. This result agrees with the theoretical predictions for diffusion-limited cluster-cluster aggregation, which state that for tÆ€ the value of k=1 for a low colloidal particle concentration. Instead, we show in Fig.7 that this process cannot be asymptotic with time for larger CP concentrations (see Fig.7b).

Figure 7 The a) mean cluster size S(t) and b) the largest cluster size Smax in time for different CP concentrations. The M=1.2¥10⁶ particles were simulated. Linear fits with kª0.5, kª1 are depicted in a).

Due to rotational diffusion in FPM, application of MD-FPM two-level model for high concentration of the CP particles should accelerate the initial agglomeration process. However, the time required for formation of the seeds is very short (about 1000 timesteps, comparing to total simulation time about 10⁵ timesteps). It is much shorter than the time needed for producing the miscellar phase. Therefore, we may expect that for the same resolution the MD-FPM model would be not as efficient as the combined MD-DPD scheme.

Multiresolution structures developed in mixing induced by the Rayleigh-Taylor instability

The great advantages of the three-level model become obviously clear in the simulation of multiscale phenomena associated with the mixing induced by the Rayleigh-Taylor instability. This process involves the CP particles in the solvent.Ý Let us assume that the upper part of the computational box is filled up with heavy (H) CP particles and the lower part with lighter (L) SP solvent particles. The heavy layer is ten times thinner than the lighter one. The collision operator is given by Eq.(17).Ý

By using MD-DPD model for a medium sized particle system comprising of M=10⁶ particles of both types, we need one IBM-SP processor running a week for integrating 10⁵ timesteps.Ý However, by using MD-DPD model (with r_cut=2.5¥l) more than 90% DPD particles simulating the solvent will not contribute to the small-scale phenomena occurring on the CP-SP interface. Therefore, it is reasonable to employ FPM with a shorter cut-off radius (r_cut=1.3¥l) in order to reduce greatly the computational time spent for calculating SP-SP interactions. We assume that SP particles will be treated as DPD (with r_cut=2.5¥l) being only in the closest neighborhood of CP particles. The parameters of DPD and FPM interactions were matched by using the kinetic theory equations [EspaÒol, 1998] for the same transport coefficients. The same solvent particle can assume dual roles (DPD or FPM), depending on whether CP particle lies within the interaction range.

Figure 8 The formation of hexagonal structures of circular micelles is observed in front of the spikes.Ý See the supplementary materials to observe the difference in scales. On the right, is shown the number of micelles in time for two particle systems at different temperatures.

In Fig.8 and its zoomed-out portion we show both the large-scale R-T patterns and the small-scale nucleation of crystallites in the front of CP spikes. As displayed in Fig.8b, in the beginning the number of micelles increases slowly on the CP-SP interface due to diffusion. The nucleation process is accelerated in the non-linear regime of the Rayleigh-Taylor instability. The length of the interface separating two immiscible DPD fluids during the mixing induced by the Rayleigh-Taylor instability increases as t² [Dzwinel and Yuen, 2000, 3]. The nucleation process becomes faster. The number of micelles increases with time as t^kwith k approximately equal to 5/2. The particular value of k does not depend on the temperature of the system.

2-D ÝÝ3-D

Figure 9. A snapshot from simulation of an accelerated colloidal cluster falling in a long box in 2-D (about 1 million particles are simulated) and 3-D (two million particles are simulated). The MD-FPM models were employed.

Since the CP-SP interactions do not have any singularities and can be described by FPM-like forces, we can investigate the dispersion of colloidal agglomerate in fluid over the mesoscale. Unlike in the liquid-liquid case in which the mixing starts spontaneously by the thermal fluctuations, the solid-liquid system does not mix so well. This is due to the energy dissipated in colloidal systems and the resistance offered to the mixing by the attractive molecular forces. The fragmentation occurs provided that the colloidal system is permeable and undergoes excessive wetting. Three types of fragmentation structures are found: rupture, erosion and shatter [Ottino et. al., 2000]. The microstructures produced during mixing are completely different from the bubbles and spikes detected in classical Rayleigh-Taylor instability [Taylor, 1950]. Instead, one can observe multi-scale structures consisting of large fishbone clusters made of a long threads or more oblate agglomerates.

The snapshots from simulations involving more than 1 million particles are displayed in Figs.9, 10. One can discern the places where more vigorous flow occurs, producing slender thread-like structures. The fishbone fragmentation of large clusters, as shown in Fig.9, is caused by the flow dynamics. We depict in Fig.10 that the threads go along with the flow streamlines, while the oblate aggregates are found in stagnant regions. After a period of vigorous fragmentation, the energy of flow damps out and the threads shrink and then agglomerate. The aggregation of colloidal beds occurs due to collisions between the aggregates, which are perturbed by the flow, or due to cohesion forces. The fishbone structures in this mixing process come from a quasi equilibrium situation due to the dynamical encounter between the colloidal particles and the overall flow induced by the global fluid dynamics.

Ý2-D

3-D

Figure 10. A snapshot from simulation of the mixing induced by the Rayleigh-Taylor instability in 2-D and fragmentation of a colloidal slab for MD-DPD two-level model in 3-D. Fragmentation and agglomeration patterns can be easily observed. For the 2-D picture the Huygens System 2.2.1 [www.svi.nl/huygens2.html] was implemented for data extraction. The densest regions are in dark blue and green, and the regions of the lowest density in dark red. The colloidal slab accelerated in 3-D cubic box is fragmented due to the shatter mechanism. See the supplementary materials for movie.

Concluding Points

Our discrete-particle simulations have clearly revealed the prospects of the three-level model comprising of: molecular dynamics (MD), dissipative particle dynamics (DPD) and fluid particles method (FPM). We have delineated the many significant advantages of the three-level model over other mesoscopic methods used for complex fluid simulation. These advantages are:

1) The three-level model is an homogeneous, particle-based model, which can operate over diverse spatio-temporal scales, which extend from a micelle, colloidal crystals, droplets, to large scale features.

2) This model is consistent with the nonlinear fluctuating hydrodynamic equations [Breuer and Petruccione, 1993]. Since thermal noise is introduced consistently, it can be used for studying non-equilibrium thermal fluctuations in hydrodynamic systems.

3) The particles evolve in a gridless fashion in real-time, thus allowing for a real-time, realistic visualization.

The principal results obtained in the modelling of complex fluids with discrete particles can be summarized as follows :

1) Using two-level model consisting of DPD and MD, we have observed the creation of colloidal crystals with different phases: ìhydrophilicî, lamellar and ìhydrophobicî. By changing the balance of the scaling factors in the CP-SP and CP-CP interactions, we can discern clearly the transition of one phase to the other. The metastable regimes with two-phases coexistence are detectable.

2) By changing the character of CP-SP interactions from conservative to dissipative, we can investigate the dynamics of colloidal aggregates. We show that in the dissipative solvent with high concentration of colloidal particles, the growth of mean cluster size can be described by the power law S(t)µt^k. In the intermediate DLA regime, for which k=‡ this lasts for a long time. The character of clusters growth varies with time and the exponent k shifts at longer times from ‡ to ª1. We show that this process cannot approach asymptotic limit in time for larger CP concentrations.

3) The three-level (MD-DPD-FPM) simulation of mixing associated with the Rayleigh-Taylor instability of colloidal particles in a solvent shows that the colloidal arrays may appear not only in the well mixed colloidal suspension but also on the mixing front. An increase in the number of micelles in time during mixing occurs faster than the increase in the length of the CP-SP interface.

4) In simulating the dispersion of colloidal slab, three types of fragmentation structures can be developed: rupture, erosion and shatter.

However, there are also some drawbacks of DPD approach.

1) The viscosities computed through the kinetic theory may not coincide with the input values [Marsh, et. al., 1997]. For the fluid particle model (and also DPD) transport coefficients has been given in the limit where no conservative forces are present.

2) The pressure is the key element in the phase separation process involving colloidal crystals, colloidal aggregation and mixing. The pressure and the sound velocity c are coupled with scaling factor P of conservative component of the DPD collision operator. As shown in [Dzwinel and Yuen ,2000, 1], this value places an upper limit on the mass S of a single particle in DPD fluid.

3) Both DPD and FPM cannot simulate the free-surface problem in fluid dynamics, because of the repulsive nature of the DPD and FPM potentials.

Both the DPD and FPM methods face a conceptual problem. This dilemma is closely related to the lack of a precise definition o the actual physical time and length scale, which DPD fluid model can describe. Introduction of volume particles as a relevant dynamical variable [EspaÒol et. al,Ý 1999] may possibly solve the problem of matching [Groot and Warren, 1997].

Acknowledgements

Support for this work has been provided by the Energy Research Laboratory Technology Research Program of the Office of Energy Research of the U.S. Department of Energy under subcontract from the Pacific Northwest National Laboratory and Polish Committee of Scientific Research (KBN).

References

[Beazley and Lomdahl, 1996]

Beazley, D.M..Lomdahl, P.S Gronbech-Jansen, N. Giles, R. Tomayo, P., Parallel Algorithms for Short Range Molecular Dynamics, in World Scientificís Annual Reviews of Computational Physics III, 119-175, (1996).

[Boryczko et. al., 2001]

Boryczko, K., Dzwinel, W., Yuen, D.A., Parallel implementation of fluid particle model for simulating complex fluids in the mesoscale, Concurrency: Practice and Experience, submitted for publication June 2001

[Breuer and Petruccione, 1993]

Breuer, H.-P. Petruccione, F., A master equation description of fluctuating hydrodynamics, Physica A, 192, 569-588 (1993).

[Cahn and Hillard, 1958]

Cahn, J.W. and J.E. Hillard, Free energy of nonuniform system. I. Interfacial free energy, J. Chem. Phys., 28, 258-267, (1958).

[Chopard and Droz, 1998]

Chopard B, Droz M. Cellular Automata Modelling of Physical Systems. Cambridge University Press, (1998).

[Coveney and EspaÒol, 1997]

Coveney, P. V., EspaÒol, P, Dissipative particle dynamics for interacting multicomponent systems, J.Phys.A:Mathematical and General, 30, 779-784, (1997).

[Daniel and Audebert, 1999]

Daniel, J.,C., and Audebert, R., Small Volumes and Large Surfaces: The World of Colloids, in Soft Matter Physics, Eds. Daoud, M., and Williams C.E., Springer Verlag,Ý 320, (1999).

[Dzwinel and Yuen, 2000, 1]

Dzwinel, W., Yuen, D.,A., Matching macroscopic properties of binary fluid to the interactions of dissipative particle dynamics, International Journal of Modern Physics C, 11/1, 1-25 (2000).

[Dzwinel and Yuen, 2000, 2]

Dzwinel, W., Yuen, D.A., A two-level, discrete-particle approach for simulating ordered colloidal structures, J. Colloid Interface Science, 225/1,179-190, (2000).

[Dzwinel and Yuen, 2000, 3]

Dzwinel, W., Yuen, D.A., Mixing driven by Rayleigh-Taylor Instability in the Mesoscale Modelled with Dissipative Particle Dynamics, International Journal of Modern Physics C , 12, 1, 91-118, (2001)

[EspaÒol, 1998]

EspaÒol, P., Fluid particle model, Physical Review E, 57/3, 2930-2948, (1998).

[EspaÒol et. al,Ý 1999]

EspaÒol, P., Serrano M., Ottinger, H.Ch., Thermodynamically Admissible Form for Discrete Hydrodynamics, Phys. Rev.Lett, 83/22, 4542-4545, (1999).

[Flekk¯y and Coveney, 1999]

Eirik G. Flekk¯y and Peter V. Coveney, From Molecular Dynamics to Dissipative Particle Dynamics, Phys. Rev. Lett. 83, 1775-1778, (1999)

[Gast and Russel 1998]

Gast, A.,P., Russel,W.,B., Simple Ordering in Complex Fluids, Physics Today, December, 24-30, (1998).

[Grier and Behrens, 2001]

Grier, D., G., and Behrens, S.,H., Interactions in colloidal suspensions, in. Electrostatic Effects in Soft Matter and Biphysics, Ed. Holm C., Keikchoff P., Podgornik R., Kluwer, 2001.

[Groot and Warren, 1997]

Groot R.D. and Warren P.B., Dissipative Particle Dynamics - bridging the gap between atomistic simulation and mesoscopic simulation, J. Chem. Phys., 107, 4423 (1997).

[Haile, 1992]

Haile, P.,M., Molecular Dynamics Simulation, Wiley&Sons, New York, (1992).

[Hoogerbrugge and Koelman, 1992]

Hoogerbrugge, P., J., Koelman,J.,M.,V.,A., Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics, Europhysics Letters, 19, 3, 155-160, (1992).

[Kremer, 2000]

Kremer, K., Computer simulations in soft matter science, in Soft and Fragile Matter, Eds. Cates, M.E., and Evans, M.,R., Proceedings of the Fifty Third Scottish Universities Summer School in Physics, A NATO Advanced Study Institute, 393 p., (2000).

[Libersky et. al, 1993]

Libersky,L.,D., Petschek, A.G., Carney,T.C., Hipp,J.,R., Allahdadi, F.A., High Strain Lagrangian Hydrodynamics, Journal of Computational Physics, 109/1, 67-73, (1993).

[Larson, 1999]

Larson R.G., The Structure and Rheology of Complex Fluids, Oxford University Press Inc., New York (1999).

[Marsh, et. al. 1997]

Marsh, C., Backx, G., Ernst, M.H., Static and dynamic properties of dissipative particle dynamics, Physical Review E, 56, 1976, (1997).

[Meakin, 1998]

Meakin, P., Fractals, scaling and growth far from equilibrium, Cambridge University Press,Ý (1998).

[Meakin, 1988]

Meakin, P., Fractal aggregates, Advances in Colloid and Interface Science, 28, 249-331, (1988).

[Ottino et. al., 2000]

Ottino, J.,M., De Roussel, P., Hansen, S., Khakhar, D.,V., Mixing and dispersion of Viscous Liquids and Powdered Solids, Advances in Chemical Engineering, 25, 105-204, 2000.

[Poon, 2000]

Poon, W., A day in the life of a hard-sphere suspension, in Soft and Fragile Matter, Eds. Cates, M.E., and Evans, M.,R., Proceedings of the Fifty Third Scottish Universities Summer School in Physics, A NATO Advanced Study Institute, 393 p., (2000).

[Prausnitz and Wu, 2000]

Prausnitz, J.,Ý Wu, J., Resolving mysterious particle attraction in colloids,Ý En Vision, 16/4, 18-19, 2000.

[Rothman and Zaleski, 1997]

Rothman, D.H., Zaleski, S., Lattice-Gas Cellular Automata: Simple Models of Complex Hydrodynamics, Cambridge University Press, (1997).

[Taupin, 1999]

Taupin, C., Physicochemistry of surfactants, in Soft Matter Physics, Eds. Daoud, M., and Williams C.E., Springer Verlag,Ý 320, (1999).

[Tanaka and Araki, 2000]

Tanaka, H., and Araki, T., Simulation Method of Colloidal Suspensions with Hydrodynamic Interactions: Fluid Particle Dynamics, Phys. Rev.Lett., 85/6, 1338-1341, (2000).

[Taylor, 1950]

Taylor, G.I., The instability of liquid surfaces when accelerated in a direction perpendicular to their plane, Proc. R. Soc. London, Ser. A, 201, 192-201, (1950).

Ý[Yaman, Jeppesen and Marques, 1998]

ÝYaman, K., Jeppesen, C., Marques, C., M., Depletion forces between two spheres in rod simulation, Europhys.Lett, 42 (2), 221-226, (1998)