Brownian Dynamics of a Suspension of Particles with Constrained

Article pubs.acs.org/Langmuir

Brownian Dynamics of a Suspension of Particles with Constrained Voronoi Cell Volumes John P. Singh,*,† Stuart D. C. Walsh,‡,† and Donald L. Koch*,† †

School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, New York 14853, United States Department of Geology and Geophysics, University of Minnesota-Twin Cities, Minneapolis, Minnesota 55455, United States

‡

ABSTRACT: Solvent-free polymer-grafted nanoparticle fluids consist of inorganic core particles fluidized by polymers tethered to their surfaces. The attachment of the suspending fluid to the particle surface creates a strong penalty for local variations in the fluid volume surrounding the particles. As a model of such a suspension we perform Brownian dynamics of an equilibrium system consisting of hard spheres which experience a many-particle potential proportional to the variance of the Voronoi volumes surrounding each particle (E = α(Vi−V0)2). The coefficient of proportionality α can be varied such that pure hard sphere dynamics is recovered as α → 0, while an incompressible array of hairy particles is obtained as α → ∞. As α is increased the distribution of Voronoi volumes becomes narrower, the mean coordination number of the particle increases and the variance in the number of nearest neighbors decreases. The nearest neighbor peaks in the pair distribution function are suppressed and shifted to larger radial separations as the constraint acts to maintain relatively uniform interstitial regions. The structure factor of the model suspension satisfies S(k=0) → 0 as α → ∞ in accordance with expectation for a single component (particle plus tethered fluid) incompressible system. The tracer diffusivity of the particles is reduced by the volume constraint and goes to zero at ϕ ∼ 0.52, indicating an earlier glass transition than has been observed in hard sphere suspensions. The total pressure of the suspension grows in proportion to (αkBT)1/2 as the strength of the volume-constraint potential grows. This stress arises primarily from the interparticle potential forces, while the hard-sphere collisional contribution to the stress is suppressed by the volume constraint. share of the fluid volume tethered to its surface, it must exclude exactly one neighbor leading to S0 = 0. In a suspension of hard spheres of radii a, the effective volume occupied by each particle is equal to the sum of its own volume (equal to (4/3)πa3) plus the surrounding void space− which depends on the volume fraction and exact location of the particle’s nearest neighbors. For most suspensions subject to Browian motion, the effective volume of each particle fluctuates as the particles move around due to the thermal motion. Nonetheless, this volume or the free volume outside the hard cores is often used in models of the dynamics of suspensions (for example, see Frankel and Acrivos,6 Krieger and Dougherty,7 and Batchelor8 among others). In this paper, we describe a model of particulate suspensions that conserves the effective local volume occupied by each particle and characterize the equilibrium structure and dynamics of such a model suspension. This system may be regarded as a colloidal suspension in which the nanoparticles carry a constant volume of fluid around them. In the model of NOHMs developed here, we will assume that the fluid volume tethered to a particle is contained within

1. INTRODUCTION Recent advances in the synthesis of hairy nanoparticles have led to the development of many rheologically interesting materials. One such class of materials consists of self-suspended nanoparticles.1,2 These materials have densely grafted oligomers on the surface of inorganic nanoparticles and are devoid of any other solvent, hence the term “self-suspended”. Also called Nanoparticle-Oligomer Hybrid Materials (NOHMs) − these materials have been shown to exhibit a range of rheological properties.3 Some systems exhibit a yield stress and soft glassy rheology, while others relax to an equilibrium state as indicated by the presence of a zero shear-rate Newtonian viscosity. Density functional based theory4 predicts and molecular dynamics simulations5 confirm that such a system is essentially incompressible and the structure factor goes to zero with decreasing wavevector. The structure factor at zero wavenumber S(k = 0) = S0 = 1 + n

∫ [g(r) − 1]d r

(1)

measures the excess number of particles in a region surrounding a test particle. When S0 = 0 as is the case in a solvent-free nanoparticle fluid, each particle occupies the same volume due to the incompressible nature of the core and the oligomeric fluid tethered to its surface. Since each particle carries its own © 2015 American Chemical Society

Received: January 22, 2015 Revised: May 14, 2015 Published: June 2, 2015 6829

DOI: 10.1021/acs.langmuir.5b00274 Langmuir 2015, 31, 6829−6841

Article

Langmuir the Voronoi volume surrounding the particle or there is an equal exchange of oligomers between a particle and its neighbors’ Voronoi cells and will therefore apply an energy penalty for variations in the Voronoi volume. This model is more restrictive than the actual tethering process because the oligomer brushes of neighboring NOHMs particles can overlap and allow variations of the Voronoi volume due to exchange of fluid volume by near neighbors. However, the typical radius of gyration of the oligomers is smaller than the core radius, and so the exchange of volume is local and modest.3 We consider a “particle” as the entity comprising both the core nanoparticle and the fluid contribution of the oligomers attached to its surface. Such a system of particles necessarily has constraints which restrict its motion as particles fill space completely and are always in contact with each other. A single particle cannot be moved without disturbing or changing the volume of the neighbors. In other words, there are no loose particles or “rattlers” in the system. The most dramatic effect can be observed in a hypothetical case of a one-dimensional system of such particles. In this case, the length of each particle is its volume, and the core-to-core distance between the particles is equal to its length. In such a scenario, motion of any particle has to be accompanied by moving all the other particles in the same direction, otherwise no motion is possible. In higher dimensions, similarly, motion of one particle is expected to be accompanied by concerted motion of at least some of its nearest neighbors. At any given time, the constraint of volume conservation may impose restrictions on motion of particles in certain directions, while they may move without constraint in other directions within the 3N-dimensional phase space of the configurations of the N particles. While a crystal or periodic array would exhibit a constant Voronoi volume, electron micrographs (see Figure 1) of NOHMs1−3,9

Figure 2. Two dimensional illustration of the model (these images are from an actual 2D implementation of the model). a) NOHMs (Nanoscale Organic Hybrid Materials) are represented as a set of spherical nanoparticles, embedded in a fluid comprised of their surface oligomers. Nanoparticles are prevented from overlapping by a hard sphere constraint, while their attached oligomers are assumed to occupy the Voronoi volume associated with the nanoparticle’s center minus the region occupied by the particle’s spherical core (illustrated by close-up in b).

the tethering of oligomers to the surface of a particle slows their segmental relaxation as shown by dielectric spectroscopy measurements in NOHMs.10 This slowing would be further accentuated if the tethered molecules were long entangled polymers as in a star polymer.11,12 In contrast, we model the frictional resistance to the motion of the cores by a simple Stokes drag law. As captured by density functional theory13 and confirmed by experiments3 as well as molecular dynamics,14 the increasing oligomer configurational entropy penalty for filling interstitial spaces as the oligomer molecular weight is reduced suggests that NOHMs with short oligomers can be more viscous and have slower dynamics. Finally, the detailed chain structure captured by coarse-grained15 or atomic16 molecular dynamics of NOHMs influences their structure and dynamics. The intent of the present study is to capture through the simplest possible model the effects of one essential feature of NOHMs, the requirement that the fluid tethered to the particle surfaces remain within the vicinity of the particle leading to a suppression of long-range density fluctuations. While this feature is conceptually simple, its implementation in a simulation is made challenging by the many-body nature of the constraint. It is important to note that if each particle conserves the total effective volume (core + oligomers) at all times, the structure factor at zero wavenumber goes to zero (eq 1). One way in which a system can achieve this is by crystallizing. As we will show, the model system studied here does not show signs of any long-range order. In fact, a density functional theory model of NOHMs that includes a soft spring representation of the oligomer interactions remains disordered beyond the volume fraction 0.5 leading to transition to an ordered state in a hard sphere system.17 While the details of the model in Yu and Koch17 differ from those in the present model, the primary origin of the delayed transition, the inability of a solvent-free nanoparticle system to exhibit coexisting ordered and disordered phases with different volume fractions is present in both models. In this paper, we first outline the details of the model based on Voronoi tesselation. We introduce an energy penalty based on the Voronoi volume of each particle which is controlled by a parameter, α. For α = 0 we obtain a hard sphere suspension, and α → ∞ results in a suspension with each particle having a fixed Voronoi volume. We then develop equilibrium Brownian

Figure 1. TEM of 14 nm diameter silica based nanoparticle-organic hybrid materials (Bourlinos et al.2). Reproduced with permission from Wiley Publishers.

indicate that they have uniform interparticle spacing but are disordered. Fernandes et al.9 have measured the Voronoi volumes of the 2D TEM scans of the experimental system and report a relatively narrow distribution of the interparticle distances between the particles with a median much larger than the core size. We have found that the proposed Voronoivolume-preserving colloidal model also exhibits a disordered structure. For example, see Figure 2 for an example of a particle configuration produced by a two-dimensional version of the model. While the present study is motivated by self-suspended nanoparticle fluids, it is not intended to capture all the factors influencing the structure and dynamics of a particular system such as nanoparticle-oligomer hybrid materials. For example, 6830


Article

Langmuir

dynamics code (see Foss and Brady,24 for example). Hard sphere suspensions have been well characterized in the literature. It has been shown that hard spheres generally have a distribution of Voronoi volumes peaked at the average volume of each particle in the suspension. For the Voronoi cell properties of equilibrium hard sphere suspensions see Kumar and Kumaran,21 Richard et al.,25 and Kumar and Kumaran.20 To evolve the suspension from a hard sphere configuration to one having a monodisperse volume, we introduce a contribution to the potential energy whenever the particle has a volume not equal to the mean volume. Specifically, we write

dynamics simulations of this system in Section 2.1. In Section 3 we report the structural information on the system, including mean coordination number, the structure factor and pair distribution function, and the equilibrium dynamics of the system, including the long time diffusivity and the total pressure.

2. MODEL DESCRIPTION In this section, we describe the model used to represent the NOHMs particle assembly. For simplicity, the nanoparticles are assumed spherical in shape and equal in size. It is assumed that the grafting density of oligomers on each particle is the same, thus each particle occupies the same volume. A hard sphere constraint is imposed on the inorganic cores to prevent them from overlapping. We use the convenient framework of the Voronoi tessellation to define the nearest neighbors of a particle and the volume occupied by the particle. Each particle’s oligomers are filling the space given by the Voronoi cell of the nanoparticle’s center minus the region occupied by the spherical core (Figure 2). Voronoi tessellations have several unique properties which make them useful for our purpose. They are by definition space filling, or, in other words, there are no void spaces. The faces of the Voronoi cell are defined by the planes which are the perpendicular bisectors of the lines joining the centers of two cores. The smallest closed polyhedron around any core formed by these faces is the Voronoi cell of the particle. All points inside the Voronoi cell of a particle are by definition closer to that particle than to any other particle, and the points on the faces of a Voronoi cell are equidistant from two particles. Conversely, Voronoi vertices of a nondegenerate tessellation are equidistant from four particles and edges are equidistant from three particles. So-called degenerate tessellations occur if five or more particles are equidistant from a Voronoi vertex. When such a degeneracy was encountered, small random perturbations (with a maximum displacement 10−6 times the particle radii) were given to the particles. The dual graph of a Voronoi tessellation is a Delaunay triangulation. For given positions of the particles in 3D, a Delaunay triangulation divides the domain into space-filling tetrahedra with particle positions for vertices. The tetrahedra have the property that the circumsphere of any tetrahedron does not contain any other particle. Moreover, the circumcenter of each tetrahedron in the Delaunay triangulation defines a vertex of the dual Voronoi tessellation. A review of the general properties of Voronoi tessellations and Delaunay triangulations is provided in the classical paper by Aurenhammer.18 Voronoi tessellations are frequently employed to characterize the geometrical or structural properties of hard sphere packings, foams, suspensions, etc. (see Finney,19 Kumar and Kumaran,20 Richard et al.,21 Kraynik et al.,22 and Troadec et al.23 among others). In this work, we use Voronoi as a tool to generate a desired property of a suspension, namely, to preserve the volume of each particle. Using the method developed in this work we can develop suspensions which contain monodisperse particles in terms of the volume occupied. Kraynik et al.22 have previously created monodisperse foams by first generating hard sphere distribution at relatively high concentration (near random close packing). The method developed in this work is more amenable to studying Brownian suspensions which preserve the volume during the time evolution of the suspension. The idea behind the generation of constant volume suspensions is simple. The initial configuration is taken as the hard sphere equilibrium generated by the typical Brownian

N

Ev =

N

∑ Evi = ∑ α(Vi − V0)2 i=1

(2)

i=1

where the summation over i is over all the particles, Vi is the current volume of particle i, and V0 is the mean volume of the particles. α is a constant controlling the strength of the constraint. Evi is the contribution to the total energy due to the volume of particle i. In other words, we hypothesize a quadratic dependence of the energy on the volume yielding a force which linearly increases with the volume difference from the mean. The quadratic energy function is an approximation locally for any sufficiently smooth energy function. If we assumed that each particles’ tethered oligomers were strictly constrained to its own Voronoi volume, the energy penalty could be considered to result from a bulk elastic modulus of the tethered fluid. If we realize that core particles may swap oligomers across their Voronoi cell boundaries, the energy penalty could include the configurational entropic penalty for the oligomers from one particle to stretch and help fill the larger Voronoi volume of a neighboring particle. Note that the Voronoi volume is determined by the position of a particle and all its nearest neighbors. As a result the potential energy of the system does not consist of a sum of pairwise potential energies. Nonetheless the direct potential interactions do remain local; this is a reasonable assumption given the limited length of the oligomers. In this work, Brownian Dynamics is employed to study the equilibrium properties. As demonstrated in Section 3, by using this approach the width of the volume distribution of particles is reduced by 2 orders of magnitude compared to a hard sphere distribution with equivalent values of the particle volume fraction ϕ. 2.1. Numerical Implementation. An introduction to Brownian Dynamics simulations for hard spheres is given in Foss and Brady,24 Rastogi et al.,26 and Heyes and Melrose27 among others. We will discuss the implementation of the constraint given by eq 2 in the context of Brownian Dynamics. At every time step, each particle i undergoes a displacement Δx i =

1 FIiΔt + X i 6πμa

(3)

where the first term is a deterministic displacement due to the interparticle interaction force FIi acting on particle i, a is the particle radius, μ is the fluid viscosity, and Xi is a stochastic displacement with zero mean

⟨X i⟩ = 0

(4)

and a variance reflecting the particle’s diffusion ⟨X iX i⟩ = 2ID0(Δt )1/2 6831

(5) DOI: 10.1021/acs.langmuir.5b00274 Langmuir 2015, 31, 6829−6841

Article

Langmuir with D0 = kBT/(6πμa) being the Stokes−Einstein diffusivity. For simplicity, we have modeled the viscous drag on the particle as a simple Stokes drag without hydrodynamic interactions. The interparticle force FIi = FIvi + FIHSi includes a contribution due to the Voronoi volume constraint energy given by eq 2, i.e.

of Voronoi volumes for hard-sphere suspensions (see Section 3 and ref 25), and the gradual reduction of this variance as α is increased. The reproducibility of the results of the dynamics for these particle numbers is reflected in the similarity of the diffusivity obtained from two simulations with different time steps shown in Figure 9.

n

FIvi = −∇i Ev = −2α((Vi − V0)∇i Vi +

3. RESULTS We have examined the equilibrium structure, self-diffusion, and stress (total pressure) of suspensions with a range of volume fractions ϕ from 0 (point particles) to 0.50. We present results with lengths nondimensionalized by the particle radius a, time by the diffusion time a2/D0, and energy by kBT. In addition, we characterize the strength of the volume constraint using αv = α⟨(V−V0)2⟩HS which represents the energy per particle that would arise due to the volume constraint if the particles had a hard sphere distribution. Here ⟨(V−V0)2⟩HS is the variance of the Voronoi volume in a hard sphere distribution. We expect the Voronoi volume constraint to begin to alter the structure when αv ≥ 6(1) and the Voronoi energy in a hard sphere system is comparable with the thermal energy. The Voronoi volume distribution for the equilibrium hard sphere suspension has been computed previously by Kumar and Kumaran25 using a particle swelling algorithm, and the variance of the Voronoi volume is summarized for the range of volume fractions considered in this study in Figure 3.

∑ (Vj − V0)∇i Vj) j=1

(6)

where ∇i is the gradient operator based on the position xi of the ith particle. The first term results from the change in particle i’s Voronoi volume due to its own displacement, and the second term is a sum over the n Voronoi neighbors of particle i of the change in the energy associated with particle j’s Voronoi volume due to the displacement of particle i. Equation 6 does not represent pairwise additive potential forces because the Voronoi volumes of particles i and j depend on the positions of their other Voronoi neighbors in addition to their own positions. Equation 6 will lead to a coordinated motion of the particles in a neighborhood toward a state of equal Voronoi volumes V0. It is important to note here that V here refers to Voronoi volume of the particles and not interaction potential. For the details of numerical implementation of the volume constraint, please see the Appendix. The interparticle potential force acting on particle i that prevents overlap of the hard core particles involves a sum over the m particles with which particle i collides in a given time step m

FIHSi =

∑ FIHSij (7)

j=1

For pure hard spheres this force is discontinuous and poses problems for implementation in computer simulations models. Generally, hard sphere interactions are modeled using a steep function which gives rise to a repulsion at a small distance between the particles. This can occasionally give rise to very large forces and limits the time step which can be used in the Brownian Dynamics simulation. Heyes and Melrose27 have developed a potential free version of the hard sphere force. In their model the hard spheres are displaced according to the Brownian displacements imposed, and the particles are allowed to overlap in any given time step. Then the overlapping particles are moved back to their contact point by displacing them by equal and opposite distances along the line joining their centers. This motion of the particles back to their contact point can be thought of as a resultant of a force acting along the vector joining the centers with the magnitude given by FIHSij = 6πμa

Figure 3. Variance of the volume distribution as a function of particle volume fraction for hard sphere suspensions: present simulations (squares); results of Kumar and Kumaran25 (diamonds).

Our results using BD are also shown for comparison. In general, good agreement between the two results indicates validity of the implementation of our numerical technique. We varied the time step between Δt = 10−4 and 10−5 and found that Δt = 5 × 10−5 was sufficient to give results independent of time step at most volume fractions, while Δt = 10−5 was used for highly concentrated suspensions (ϕ ≥ 0.45). We first generate hard sphere configurations at a given concentration and equilibrate the hard spheres using Brownian dynamics. Then we switch on the volume constraint. Initially, the value of αv was kept low. When the equilibrium is reached for a given value of αv, we gradually increase αv and reequilibrate and continue the process until the desired value of the constraint energy coefficient is reached. 3.1. Structure. Figure 4 shows the variance of the Voronoi volume (normalized by the square average volume occupied by each particle) of the particles for various values of 1/αv at ϕ = 0.134 and 0.262. The plot clearly shows that as the constraint

Δx ij 2Δt

(8)

where −Δxij is the maximum overlap between the particles, and Δxij/2 is the displacement of each particle required to bring the particles back to point of contact. According to Heyes and Melrose, this pairwise approach of removing overlaps can resolve all hard sphere contacts even at high densities. A similar approach was also used later by Foss and Brady.24 The simulations were performed with 250 particles in a periodic cell for ϕ = 0.262 and 1000 particles for ϕ = 0.134. The larger number of particles for the low volume fraction case were required in order to accurately produce the large variance 6832


Article

Langmuir

Figure 4. Plot of variance of the Voronoi volume distribution normalized by the square of the average volume per particle as a function of 1/αv. (◊) at ϕ = 0.134 and (○) at ϕ = 0.262. Lines are linear fits to the data for αv ≥ 1. Inset shows more clearly the data for 1/αv ≤ 0.5.

is increased, the width of the distribution decreases, as expected. The lines are linear fits to the data for the range 1/αv < 0.4. We observe that the variance of the Voronoi volume shows much weaker dependence on 1/αv when αv < 1. This corroborates our expectation that αv ≥ 6(1) is required for the volume constraint to affect the structure and dynamics of the suspension. The linear relationship between the variance of the volume distribution and the strength of the constraint or ⟨(V − V0)2 ⟩ ∼ 1/αv ∼ 1/α

Figure 5. Histogram of Voronoi volume with bin size dV = 0.6 × St. dev(V) as a function of (V−V0)/St. Dev(V). (a) ϕ = 0.134; (□) are at αv = 6, (○) at αv = 30, and (◁) at αv = 60. (b) ϕ = 0.262; (○), (*), (□), and (▷) are at αv = 0.6, 3.5, 7, and 70, respectively. Also included is the hard sphere distribution (black diamonds). Lines are the corresponding normalized Gaussian distribution.

(9)

implies that the energy of the constraint E ∼ α⟨(V−V0) ⟩ is constant for sufficiently large values of αv. In addition, the relationship suggests that in the limit α → ∞, the variance in the Voronoi volume goes to zero. The volume constraint parameter α is most closely related to the bulk modulus of the NOHMs suspension, which has not been measured experimentally. However, a typical value of the shear modulus of NOHMs with ϕ = 0.12 is G = 11 kPa.3 Assuming that the bulk and shear moduli are comparable and therefore taking an estimate of the constraint parameter in dimensional form as α = nG, we obtain dimensionless αv ∼ 160 which is sufficiently large to lead to a fully constrained system. In Figure 5 the Voronoi volume distribution function is shown for various values of αv. Figure 5 (a) is at ϕ = 0.134, and Figure 5 (b) is at ϕ = 0.262. The distribution is plotted as a function of (V−V0)/((⟨(V−V0)2⟩)1/2). Also, shown on the plots are the normalized Gaussian defined by 2

Ng =

⎛ (V − V )2 ⎞ 0 ⎟ exp⎜ − ⎝ 2Var(V ) ⎠ 2πVar(V )

are far enough from the mean to lie in the tail of the hardsphere Voronoi volume distribution and so the Voronoi volumes populated by thermal energy are in the region where the distribution of available V values is nearly Gaussian. Second, in the limit of small volume fraction, one can expect the volume distribution to follow a Boltzmann distribution given by exp(−E/(kBT)) with the Voronoi energy E = α(V−V0)2 yielding a Gaussian distribution. At moderate volume fractions and moderately high αv the probability of realizing an available configuration follows the Boltzmann/Gaussian distribution and the number of available sites (as indicated by the hard-sphere results) is also a Gaussian function of V in the populated region of V−V0. Figure 6 shows the effect of αv on the mean coordination number (or, the number of nearest Voronoi neighbors of a particle) and the variance of the coordination number. The coordination number is equal to the number of faces of a particle’s Voronoi cell. The behavior at both concentrations (ϕ = 0.134 and 0.262) is very similar. The mean coordination number increases modestly from about 15.4 to almost 15.55 for the relatively dilute case (ϕ = 0.134) and from 15.14 to 15.26 for the denser suspension, remaining in a range characteristic of disordered suspensions. In comparison, a slightly perturbed FCC crystal has approximately 14 nearest Voronoi neighbors (see Kumar and Kumaran21 and Richard et al.20). The variance of the coordination number decreases with increasing αv and reaches a plateau. For the dilute suspension

1

(10)

where Var(V) is the variance of the Voronoi volume distribution, and V0 is the mean volume of each particle. We observe that for large enough αv the Voronoi volume follows a Gaussian distribution within the noise of the simulations. Kumar and Kumaran25 have shown that the distribution of Voronoi volumes in a hard sphere suspension is nearly Gaussian in a region surrounding the mean volume but exhibits a tail for large Voronoi volumes. The Gaussian nature of the distribution for large αv can be attributed to two factors. First, large αv yields a large energy penalty for Voronoi volumes that 6833


Article

Langmuir

Figure 6. a) Mean coordination number; b) variance of the coordination number for volume fraction ϕ = 0.134 (◊) and 0.26 (○).

(ϕ = 0.134), the plateau value is around 2.5, and for the relatively concentrated suspension (ϕ = 0.262) the value is around 1.8. The effect of the constraint is such that the particles try to create relatively uniform spacing. The particles which are almost touching each other in the hard sphere suspension will be pushed away from each other. As will be shown later, this pushes the first peak of pair-distribution function to a slightly higher r, so that more particles can be accommodated on average within the first shell of near neighbors. The Voronoi volume distribution of a hard sphere suspension also becomes somewhat more uniform with increasing volume fraction. However, such narrow distributions of the coordination number are only obtained for hard spheres when the particle volume fraction is above 0.5.21 The decrease in the variances of both Voronoi volumes and coordination numbers indicate that in some sense the effect of the constraint is to make the suspension more homogeneous compared to hard spheres. This trend is further corroborated by the pair probability density, g(r), and the structure factor, S(k). Figure 7 a, b, and c, respectively, show g(r) for ϕ = 0.134, 0.262, and 0.4 at various values of αv. The solid lines indicate the Percus−Yevick predictions for the hard-sphere distribution, while red squares are at intermediate and blue stars at high values of αv. In all cases, the first peak in g(r) is reduced with the introduction of constraint. In addition, the first peak shifts to larger radial separations for ϕ = 0.134 and 0.262. Pairwise

Figure 7. Effect of volume constraint on the pair distribution function, g(r) for various volume fractions (ϕ = 0.134, ϕ = 0.262, and ϕ = 0.4). Solid lines are Percus−Yevick hard sphere solutions; red squares are at intermediate values of αv and blue stars are at higher values. For ϕ = 0.134, αv = 0.6 and 60. For ϕ = 0.262, αv = 0.7 and 70. For ϕ = 0.4, αv = 1 and 50.

repulsive interparticle potentials used to model soft colloids in solvents also lead to an outward shift of the first peak.28 However, such purely repulsive potentials greatly enhance the size of the first peak and in some cases a second peak, while the pair distribution function for the Voronoi model returns to a form quite similar to the hard sphere distribution at large r. The most dramatic effects are observed in dilute suspensions. For ϕ = 0.134 the contact value of pair distribution is reduced by almost a factor of 2, while the first peak whose magnitude is less than 1.1 for αv = 60 is shifted outward from the close 6834


Article

Langmuir contact position to r = 1.46*2a. Touching particles will necessarily lead to smaller Voronoi volumes which are disfavored by the constraint force. The particles displaced from the close contact region are pushed outward forming a more uniform shell of nearest neighbors separated by intervening tethered fluid. The second nearest neighbor shell will also become more uniform, enhancing the second peak. Molecular dynamics simulations of NOHMs have shown a similar outward shift of the first peak and a stronger second peak.5 The height of the peaks of g(r) in the molecular dynamics simulations is much stronger for the same core volume fraction than that in the Voronoi simulation. However, much of the enhanced structure might be attributed to an effective hard sphere diameter that includes the first bead in the oligomer chain which is fixed to the particle surface.29 The influence of the Voronoi volume constraint on g(r) is smaller for higher volume fraction, and at ϕ = 0.4 a slight decrease in g for nearly touching particles is the only evident change due to the constraint. The corresponding changes in the structure factor are depicted in Figure 8 a for ϕ = 0.134, b for ϕ = 0.262, and c for ϕ = 0.4. The most evident influence of the constraint is to reduce the structure factor at small k. For sufficiently large αv, the structure factor is reduced nearly to zero for sufficiently small k so that S(k = 0) = 0 and the long-range density fluctuations are completely suppressed. For a smaller value, αv = 6 in Figure 8 a, S(k) reaches a finite plateau value that is much smaller than that for a hard sphere suspension. This case is analogous to the molecular dynamics simulations and density functional theory for a compressible tethered oligomer fluid achieved at elevated temperature in ref 4 where small but finite long-range density fluctuations arise in the compressible solvent-free nanoparticle fluid. To compare the structure factor obtained from the Voronoi model with that seen in more detailed models of NOHMs, we plot in Figure 8 (c) the results of molecular dynamics (MD) simulations5 as well as density functional theory (DFT) based calculations.29 The MD simulations are for a core volume fraction of 0.169. However, if we consider an effective hard particle radius that includes the first bead of the oligomer chain which does not move relative to the core, an effective volume fraction of ϕ = 0.46 is obtained which is similar to but slightly larger than the value ϕ = 0.4 for the Voronoi model simulations. The DFT calculations are for ϕ = 0.46 and consider an incompressible fluid of tethered soft-spring oligomers. It can be seen that all three models of NOHMs show a similar reduction of the structure factor relative to that of the hard sphere suspension for k ≤ 2. The peaks of the MD model are stronger than those of the DFT model, and this may be attributed to the influence of the structuring of the oligomer chains on the core− core interactions. The differences in the location and the height of first peak between the Voronoi volume constraint model and the DFT based calculations can be mostly attributed to the higher concentration in the DFT model. 3.2. Long Time Diffusivity. In this section we will examine the long time self-diffusivity of particles experiencing hard sphere repulsion and a Voronoi volume constraint force. To begin with, we might ask whether the limit αv → ∞ in which the volume constraint becomes a hard constraint would prevent particles from diffusing. Each particle’s Voronoi cell is constrained leading to N constraints in a system of N particles. If the particles exist in a one-dimensional space this constraint would prevent any relative motion of the particles. A onedimensional system of hard rod particles is known to undergo a growth of its mean-square displacement whose slope decreases

Figure 8. Effect of volume constraint on the structure factor, S(k). Solid line is Percus-yevick solution for hard spheres. (a) ϕ = 0.134; (◊) for αv = 6 and (○) for αv = 60. (b) ϕ = 0.262; (◊) for αv = 0.6 and (○) for αv = 60. (c) ϕ = 0.4; (○) for αv = 50. Also shown in (c) are the results of MD simulations (dotted line) and DFT based calculations (dashed blue line) with an effective volume fraction of ϕ = 0.46 for comparison. HS Percus−Yevick predictions for ϕ = 0.4 and 0.46 are shown using solid red and black lines, respectively.

with time as more and more particles must be set into motion to allow one particle to diffuse increasing the mass or resistance to the motion.30 However, with the hard Voronoi constraint, the resistance diverges at the initial instant that any particle moves. In the three-dimensional space that we consider, however, there remain 2N unconstrained degrees of freedom. This suggests that the system may remain diffusive at least at sufficiently low volume fractions. The diffusivity was determined by computing the meansquare displacement of the particles in simulations for times up 6835


Article

Langmuir to 10 times the diffusivity time. The diffusivity was estimated as the value of D = ⟨r2⟩/(6t) achieved at long times. In most of the simulations a plateau of this value was achieved long before t = 10. Figure 9 shows the diffusivity as a function of the

Figure 9. Long time diffusivity of particles as a function of particle volume fraction and the volume constraint coefficient αv. (○) is for point particles; diamonds (◊ and ⧫) for ϕ = 0.134; squares (□ and ■) for ϕ = 0.262; and triangle (△) for ϕ = 0.4. All the computations are conducted using time step of Δt = 5E−5 (shown using open symbols) except for (⧫) which were at Δt = 1E−5 and (■) which were at Δt = 1E−4. Inset shows the same data limited to αv ≤ 1, for clarity.

constraint energy parameter αv for ϕ = 0 (point particles with no hard sphere repulsion), ϕ = 0.134, 0.262, and 0.4. The results for αv = 0 are in good agreement with previous Brownian dynamics simulations of hard spheres by Foss and Brady.24 The diffusivity for point particles decays from 1 with no constraint to a value of about 0.6 as αv → ∞. This value is close to the value 2/3 that one might infer based on the fact that two-thirds of the degrees of freedom in the system remain unconstrained. For example, in a porous medium in which the fluid space consists of straight parallel planar channels with random orientation and a channel thickness large enough to allow for continuum diffusion in the pores, the solute would be restricted from diffusing in one of the three directions in each channel. Upon averaging over the channel orientation, the diffusivity in the porous medium would be 2/3 times the product of the bulk diffusivity and the porosity. This situation is not unlike that in the Voronoi simulations, where particles can diffuse much more easily along the faces of the Voronoi cells than perpendicular to the faces but the orientation of the faces varies throughout the suspension. As the concentration increases, the long time diffusivity is reduced from its point particle value as is observed in pure hard sphere suspensions. As the value of αv is increased the diffusivity decreases, as expected. We note that the time step for the range considered does not affect the diffusivity value indicating convergence with respect to the time step. For time steps between Δt = 1 × 10−4 to 1 × 10−5 the long time diffusivity value changed by less than ±1%. The behavior of the diffusivity for large constraint coefficients can be understood in terms of the mean-square Voronoi volume fluctuations for large but finite αv It was shown in the previous section that ⟨(V−V0)2⟩ ∼ 1/αv as αv → ∞. This indicates that particles moving in one of the constrained directions can achieve a displacement Δr ∼ 1/(√αv) in each time step, resulting in a contribution to the diffusivity of (Δr)2/(δt) ∼ 1/αv.

Figure 10. (a) Linear fits to diffusivity data as a function of inverse αv. Cross (×), filled triangle (▲), filled square (■), filled diamond (⧫), and star (*) at ϕ = 0, 0.134, 0.26, 0.4, and 0.5, respectively. (b) The diffusivity in the limit of α → ∞ (D∞) (shown as filled diamonds) is compared with hard sphere diffusivity (open circles), as reported by Foss and Brady24 as a function of ϕ.

Consistent with this observation, we find in Figure 10 (a) that the diffusivity is a linear function of 1/αv, and an extrapolation of this linear behavior to 1/αv → 0 can be used to obtain the diffusivity D∞ for particles with a hard volume constraint αv → ∞. The value of D∞ is plotted as a function of volume fraction in Figure 10 (b) along with diffusivities for hard sphere suspensions as reported by Foss and Brady.24 The diffusivity with the Voronoi constraint is consistently smaller than that for hard spheres and decreases with increasing particle volume fraction. It seems to be approaching zero at a volume fraction slightly above 0.5, a concentration at which the diffusivity of the hard sphere system remains finite. It is known that at high concentrations hard sphere suspensions exhibit a transition from a fluid system with a finite diffusivity to a glassy system with zero diffusivity. Early estimates of the volume fraction for the glass transition of hard spheres were ϕ0 = 0.58. There have been many attempts in the literature (for example, see Pusey and Van Megen,31 Brambilla et al.,32 and Doliwa and Heuer33) to verify this value and to characterize the dynamics of the system close to glass transition. We will examine the question of whether a Voronoi constrained system may have glassy dynamics at a smaller volume fraction than that for hard spheres. To characterize the dynamics of the suspension at the volume fraction giving a small nonzero diffusivity, we plot in Figure 11 the mean-square 6836


Article

Langmuir

Figure 11. Mean squared displacement of particles as a function of volume fraction for hard spheres, αv = 3.15 and αv = 35 at ϕ = 0.5. The straight line indicates a slope of 1.

displacement versus time for a suspension with ϕ = 0.5 for αv = 0, 0.35, and 3.5. It can be seen that in addition to reducing the magnitude of the mean-square displacement, the constraint leads to an extended period of subdiffusive dynamics in which the mean-square displacement grows slower than linearly with time. Such a period of subdiffusive dynamics is seen in hard sphere suspensions close to their glass transition.33 We will now attempt to estimate the volume fraction corresponding to the glass transition for the constrained system with αv → ∞. It has been found34 that the diffusivity near the glass transition volume fraction ϕ0 follows a power law D ∼ (ϕ−ϕ0)β. However, to accurately determine the exponent β would require simulations spanning a few orders of magnitude in the diffusivity or ϕ−ϕ0. We have only been able to obtain diffusivities accurately for D/D0 > 0.04. We therefore adopted a strategy of obtaining the glass transition for the fully constrained system through two approaches. First we estimate the glass transition ϕ0 for several values of the constraint parameter αv by plotting, in Figure 12 (a), D/D0 versus volume fraction and assuming a linear relationship D/D0 = A(ϕ−ϕ0). The resulting glass transition volume fraction is plotted as a function of αv in Figure 12 (b). The glass transition volume fraction decreases with increasing αv and reaches a value of 0.522 at the highest value αv = 141 considered. An alternative estimate of the glass transition by extrapolating the fully constrained diffusivity D∞ using a linear function of ϕ−ϕ0 for ϕ > 0.4 yields ϕ0 = 0.524. The close agreement between these two methods of estimating the glass transition suggests that the Voronoi constrained suspension has a glass transition near ϕ = 0.52 and below the value for a hard sphere suspension. 3.3. Isotropic Stress. In this subsection, we calculate the pressure in the model suspension. Since the suspension is at equilibrium, the stress is isotropic, and the pressure is obtained as one-third of the trace of the stress tensor. The stress in a system of particles with periodic boundary conditions can be calculated as35 N

τ=

Figure 12. (a) D/D0 as a function of suspension volume fraction for various values of αv. Lines are linear fits to the data. αv = 1.7, 28.2, and 141, shown using (△), (◊), and (○), respectively. (b) Estimated concentration at which the glassy behavior is approached as a function of the volume constraint.

There are two contributions to the interparticle force Fij (and therefore the stress) in the Voronoi constrained suspension, a hard-sphere collisional force FIHSij and the Voronoi constraint force FIvij. The hard sphere collisional force defined by eq 8 is pairwise additive. Although the Voronoi constraint force acting on particle i depends on the coordinates of all of its Voronoi neighbors in a coupled fashion, we can define an interparticle force between particles i and j as

FIVij = −∇ij Ev

where ∇ij is the gradient operator with partial derivatives taken with respect to the relative position rij of particles i and j holding the relative positions rik of particle i and its other Voronoi neighbors k ≠ j fixed. It should be noted that FIVij is not a pairwise force and depends on the coordinates of the other Voronoi neighbors. A similar form for the stress in a nonpairwise additive potential has been used by Zimmerman et al.36 in a molecular dynamics simulation. Figure 13 shows the hard-sphere and Voronoi constraint contributions to the stress as a function of αv for various particle volume fractions. The constraint contribution to the stress grows steeply with αv at large αv. The collisional component to the stress decreases and reaches a high αv plateau. This is consistent with the observation that the contact value of the pair-probability distribution decreases due to the volume

N

1 ⟨ ∑ ∑ rijFij⟩ Vs i = 1 j = 1

(12)

(11)

where rij is the minimum distance between particles i and j, Fij is the interparticle force, and Vs is the suspension volume. 6837


Article

Langmuir

Table 1. Contact Value of Pair Distribution Function from the Collisional Stress and from Extrapolating the Pair Probability Obtained by Binning Particle Pair Positions ϕ

αv

gc from PHS

gc from g(r)

0.134 0.262 0.40

60.0 70.0 93.2

0.68 1.33 2.78

0.71 1.29 2.88

⎛ k T ⎞1/2 δV ∼⎜ B ⎟ ⎝ α ⎠ V0

Since the constraint acts primarily along the lines of centers of two neighboring particles, the particles have a small standard deviation of relative position δr in the normal direction while being free to translate through distances of order V1/3 in the 0 two tangential directions. Thus, δr ∼ δV/V2/3 0 . Approximating the potential locally as a harmonic well in the constrained direction, the interparticle force scales as

Figure 13. Hard sphere (open symbols and left axis) and volume constraint (filled symbols and right axis) contributions to the pressure for ϕ = 0.134 (□), 0.262 (◊), and 0.4 (△).

constraint but remains finite. There are fewer configurations with touching particles that can satisfy the volume constraint, reducing the frequency of hard sphere collisions. The kinetic theory of dense gases predicts that the collisional contribution PHS to the total pressure can be related to the pair probability at particle contact gc by PHS = nkBT 4ϕgc

(14)

F∼

kBT (αkBT )1/2 ∼ δr V01/3

(15)

leading to a Voronoi constraint pressure of Pv ∼ nFNcV01/3 ∼

(13)

ϕNc(kTα)1/2 a3

(16)

where n is the number of particles per unit volume. Using this relationship one can infer the contact pair distribution function from the collisional part of the total pressure. The estimates of the contact value of the pair distribution function as a function of αv are shown in Figure 14. These values compare well with

The primary result of this scaling analysis is that the Voronoi constraint contribution to the pressure grows as α1/2 as the constraint coefficient α approaches infinity. In Figure 15, we compare the simulation results with this prediction. The simulated volume-constraint pressure contributions

Figure 14. Contact value of the pair distribution function inferred from the collisional stress as a function of αv. (△) at ϕ = 0.4, (□) at ϕ = 0.262, and (◊) at ϕ = 0.134.

Figure 15. Volume constraint contribution to the stress plotted on a log scale as a function of α for three different core volume fractions. The dark solid line has the slope of 0.5. The lines going through the data are power law fits to the data at high α. We obtain the power law index to be 0.53, 0.54, and 0.55 for ϕ = 0.134, 0.262, and 0.4, respectively.

those obtained by binning the particle pair configurations (see Figure 7) and extrapolating the pair distribution to r = 2a. See the comparison in Table 1. The monotonic rise in the constraint contribution to the pressure with increasing α can be understood in terms of the increasing stiffness of the potential forces that enforce the constraint. By supplying the thermal energy kBT to the degrees of freedom associated with the Voronoi constraint, one obtains a fluctuation δV of each Voronoi volume that scales as

have a power law dependence on α at large α with the best fits having exponents ranging from 0.53 to 0.55. A solid bold line with slope 0.5 is also shown for comparison. At small α, the deviation from the steep harmonic potential energy leads to smaller potential contributions to the pressure than predicted by the scaling. It may also be noted that for a given α the constraint pressure is lower at higher particle volume fractions. This may be attributed to the fact that particles whose motion is limited by the hard sphere 6838


Article

Langmuir

stress. Indeed oscillatory shear experimental measurements indicate that NOHMs have very large viscosities and elastic moduli.3 The self-diffusivity is reduced by the volume constraint. In the point particle limit, the diffusivity is reduced by a factor of 0.6, which is similar to the 2/3 factor that one might expect based on the fact that 1/3 of the degrees of freedom of particle motion have been constrained. As the core volume fraction increases the diffusivity decreases and remains lower than the hard-sphere suspension value. The simulation results suggest that a fully constrained system, α → ∞, experiences a glass transition at a volume fraction ϕ0 ≈ 0.52 that is smaller than the value for the hard-sphere glass transition. The self-diffusivity of NOHMs has not yet been measured experimentally. However, coarse-grain molecular dynamic simulations15 that resolve the structural and dynamical effects of freely jointed chains of 10−15 beads tethered to nanoparticles indicate that the diffusivity is reduced from the value for nanoparticles in untethered oligomer fluids by much larger factors (3−30) at ϕ = 0.02−0.2 than those predicted here. A part of the reduction in the diffusivity in molecular dynamics may come from the slower relaxation time of the tethered oligomers as seen experimentally by Agarwal et al.10 It is interesting to note, however, that the molecular dynamics like the present Brownian dynamics suggest that the diffusivity remains finite at low to moderate volume fractions as long as the chains are long enough to provide fluidity. In particular, NOHMs with 5-mer chains yielded a diffusivity that vanished at a glass transition ϕ0 = 0.2, while NOHMs with 10- and 15-mer chains retained a finite diffusivity at all volume fractions accessible with the given grafting density and oligomer volume.15 The present Brownian dynamics estimate ϕ0 ≈ 0.52 may be considered at upper bound on the glass transition for self-suspended nanoparticles, and the actual value may be decreased if a portion of the tethered molecules are immobilized and act as an effective portion of the core or the chain configurational entropy exerts stronger constraints than the requirement of equal Voronoi volumes.

constraints are not experiencing the steep Voronoi constraint potential. The resulting reduction in the constraint-induced pressure is most marked at high volume fraction, and surprisingly it can lead to a smaller total pressure at higher particle volume fractions.

4. CONCLUSION One of the most basic and unique features of self-suspended nanoparticle fluids such as nanoparticle-organic hybrid materials (NOHMs) is the requirement that the tethered molecules that fluidize the inorganic cores remain in the vicinity of the cores. In this paper, we have developed and implemented a conceptually simple approach to prescribing this constraint in a Brownian dynamics simulation of the structure, self-diffusivity, and pressure of equilibrium suspensions. The core particles experience, in addition to hard core constraints, a potential energy proportional to the variance of the Voronoi volumes of the particles. As the coefficient α that enforces this constraint is increased, the particles begin to fill space more uniformly. This is measured by a decrease in the variance of the Voronoi volumes and a reduction of the structure factor at zero wavenumber, a measure of long-range particle density fluctuations. The suspension remains disordered, however, for all α. The effect of the constraint on the pair distribution function is most obvious when the particle-core volume fraction is small. The first peak of the pair distribution function is shifted to higher r and is suppressed as the particles try to fill space “more uniformly” than the hard sphere distribution, features seen in molecular dynamics simulations of NOHMs.15 The coordination number in constrained suspensions with moderate volume fractions ϕ = 0.262 reaches values greater than 15.5, which would be seen in hard-sphere random suspension only at much higher volume fraction ϕ > 0.5. The contact value of pair distribution is reduced, and this leads to a reduction in the hard-sphere contribution to the total pressure. The overall pressure, however, increases monotonically with increasing α. At high α, the Voronoi constraint limits the particle motion along the lines of centers to a thin δr ∼ 1/α1/2 region. This constraint arises from the Boltzmann distribution of the Voronoi volumes. The constraint energy can be approximated as a steep harmonic well leading to the prediction that the potential contribution to the pressure grows in proportion to α1/2. The pressure from the Brownian dynamics results agree with this prediction. We may anticipate that when the suspension is deformed this large isotropic stress would be transformed into a large deviatoric

■

APPENDIX In this Appendix we will discuss the implementation of the interparticle force given by eq 6 for computer simulation. We start by calculating the Delaunay triangulation for the given position of the particles. There are many subroutines available freely for such a calculation. We used GEOMPACK which implements a series of local transformations of tetrahedra to

Figure 16. a) A schematic representing a face of the Voronoi cells around particle P1 and P2. O is midway on the line joining the two points. OABP1 form a tetrahedron. b) A schematic showing the geometrical relationship between particles P1, P2, and P3 which are all nearest neighbors of each other. 6839


Langmuir

■

achieve Delaunay triangulation.37 Then using the relationship between Voronoi and Delaunay we determine the Voronoi vertices by calculating the circumcenter of each tetrahedron in the Delaunay triangulation. Evaluation of the Voronoi cell volume given the position of the particle and the Voronoi cell vertices can be achieved by tesselating the cell into nonoverlapping but space filling tetrahedra and then calculating the volume of each tetrahedron. The volume of any tetrahedron is given by | a · ( b × c )| 6

Vtet =

n

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by Award No. KUS-C1-018-02 made by the King Abdullah University of Science and Technology (KAUST).

(17)

■

m(j) k=1

(18)

and n

∇Vi =

m(j)

∑ (∑ ∇Vktet ) j=1

k=1

(19)

where n is the number of nearest neighbors of particle i, and m(j) is the number of edges on the face between particle i and its jth neighbor. The sum on the index k is over all the tetrahedra formed between particle i and its jth neighbor. Note that when a particle is perturbed it causes the Voronoi cell vertices to be moved in a nontrivial way which can be determined using the fact that the Voronoi vertices are circumcenters of the simplices of the Delaunay triangulation. In 3−D the circumcenter of the tetrahedron is given by co =

REFERENCES

(1) Bourlinos, A.; Chowdhury, S.; Herrera, R.; Jiang, D.; Zhang, Q.; Archer, L.; Giannelis, E. Functionalized nanostructures with liquid-like behavior: Expanding the gallery of available nanostructures. Adv. Fun. Mater. 2005, 15, 1285−1290. (2) Bourlinos, A. B.; Herrera, R.; Chalkias, N.; Jiang, D. D.; Zhang, Q.; Archer, L. A.; Giannelis, E. P. Surface-Functionalized Nanoparticles with Liquid-Like Behavior. Adv. Mater. 2005, 17, 234−237. (3) Agarwal, P.; Qi, H.; Archer, L. The Ages in a Self-Suspended Nanoparticle Liquid. Nano Lett. 2010, 10, 111−115. (4) Hsiu-Yu, Y.; Koch, D. Structure of Solvent-Free NanoparticleOrganic Hybrid Materials. Langmuir 2010, 16801−16811. (5) Chremos, A.; Panagiotopoulos, A.; Yu, H.; Koch, D. Structure of solvent-free grafted nanoparticles: Molecular dynamics and densityfunctional theory. J. Chem. Phys. 2011, 135, 114901. (6) Frankel, N.; Acrivos, A. On the viscosity of a concentrated suspension of solid spheres. Chem. Eng. Sci. 1967, 22, 847−853. (7) Krieger, I. M.; Dougherty, T. J. A mechanism for non-Newtonian flow in suspensions of rigid spheres. Trans. Soc. Rheol. 1959, 3, 137− 152. (8) Batchelor, G. The effect of Brownian motion on the bulk stress in a suspension of spherical particles. J. Fluid Mech. 1977, 83, 97−117. (9) Fernandes, N. J.; Akbarzadeh, J.; Peterlik, H.; Giannelis, E. P. Synthesis and Properties of Highly Dispersed Ionic Silica-Poly (ethylene oxide) Nanohybrids. ACS Nano 2013, 7, 1265−1271. (10) Agarwal, P.; Kim, S. A.; Archer, L. A. Crowded, Confined, and Frustrated: Dynamics of Molecules Tethered to Nanoparticles. Phys. Rev. Lett. 2012, 109, 258301. (11) McLeish, T. C.; Milner, S. T. Branched Polymers II; Springer: 1999; pp 195−256. (12) Briels, W.; Vlassopoulos, D.; Kang, K.; Dhont, J. K. Constitutive equations for the flow behavior of entangled polymeric systems: Application to star polymers. J. Chem. Phys. 2011, 134, 124901. (13) Yu, H.-Y.; Koch, D. L. Self-diffusion and linear viscoelasticity of solvent-free nanoparticle-organic hybrid materials. J. Rheol. 2014, 58, 369−395. (14) Goyal, S.; Escobedo, F. A. Structure and transport properties of polymer grafted nanoparticles. J. Chem. Phys. 2011, 135, 184902. (15) Chremos, A.; Panagiotopoulos, A. Z.; Koch, D. L. Dynamics of solvent-free grafted nanoparticles. J. Chem. Phys. 2012, 136, 044902. (16) Hong, B.; Panagiotopoulos, A. Z. Molecular dynamics simulations of silica nanoparticles grafted with poly (ethylene oxide) oligomer chains. J. Phys. Chem. B 2012, 116, 2385−2395. (17) Yu, H.-Y.; Koch, D. L. Predicting the Disorder-Order Transition of Solvent-Free Nanoparticle-Organic Hybrid Materials. Langmuir 2013, 29, 8197−8202. (18) Aurenhammer, F. Voronoi diagrams: survey of a fundamental geometric data structure. ACM Computing Surveys (CSUR) 1991, 23, 345−405. (19) Finney, J. Random packings and the structure of simple liquids. I. The geometry of random close packing. Proc. R. Soc. London, Ser. A 1970, 479−493. (20) Kumar, V.; Kumaran, V. Voronoi neighbor statistics of harddisks and hard-spheres. J. Chem. Phys. 2005, 123, 074502.

∑ (∑ Vktet ) j=1

AUTHOR INFORMATION

Corresponding Authors

where one vertex of the tetrahedron is at the origin, and the other three are located at positions a, b, and c. Now, the major work involves calculating the change in volume of a particle’s Voronoi cell and those of its neighbors, when it is perturbed. It is important to note that each Voronoi face is shared by two neighboring particles, each edge by three particles and each Voronoi vertex is shared by four particles. Moreover each face is a convex polygon. Figure 16 (a) depicts such a face between particles P1 and P2. O is a point midway between the particles and lies on the Voronoi face. We form a tetrahedron OABP1. The whole Voronoi cell around point P1 can be divided into such tetrahedhra or Vi =

Article

1 (p ·p (p × p3) + p2 ·p2(p3 × p1) + p3·p3(p1 × p2)) 12V 1 1 2

(20)

where co is the circumcenter of tetrahedra Op1p2p3, and O is the origin of the coordinate system. The volume of the tetrahedron is denoted by V. We apply eq 19 to calculate the second term of eq 6 involving the force on particle i due to the change in volume of its neighbors. As shown in Figure 16 (b) P3 shares an edge AB with P1 and P2. P2 and P3 also share a face ABWXYZ with O′ point midway on the line joining P2 and P3. When P1 is perturbed the volume of P2 changes due to the change in the volume of tetrahedron OABP2 as well as that of O′ABP2. We maintain lists of Voronoi cell vertices of each particle, nearest neighbors, and particles contributing to each Voronoi vertex. Using these lists we can easily construct the Voronoi face between two neighbors and find other particles which share a particular Voronoi vertex with any two given neighboring particles. 6840


Article

Langmuir (21) Richard, P.; Oger, L.; Troadec, J.; Gervois, A. Geometrical characterization of hard-sphere systems. Phys. Rev. E 1999, 60, 4551− 4558. (22) Kraynik, A.; Reinelt, D.; van Swol, F. Structure of random monodisperse foam. Phys. Rev. E 2003, 67, 31403. (23) Troadec, J.; Gervois, A.; Oger, L. Statistics of Voronoi cells of slightly perturbed face-centered cubic and hexagonal close-packed lattices. Europhys. Lett. 1998, 42, 167. (24) Foss, D.; Brady, J. Brownian dynamics simulation of hard-sphere colloidal dispersions. J. Rheol. 2000, 44, 629. (25) Kumar, V.; Kumaran, V. Voronoi cell volume distribution and configurational entropy of hard-spheres. J. Chem. Phys. 2005, 123, 114501. (26) Rastogi, S.; Wagner, N.; Lustig, S. Rheology, self-diffusion, and microstructure of charged colloids under simple shear by massively parallel nonequilibrium Brownian dynamics. J. Chem. Phys. 1996, 104, 9234. (27) Heyes, D.; Melrose, J. Brownian dynamics simulations of model hard-sphere suspensions. J. Non-Newtonian Fluid Mech. 1993, 46, 1− 28. (28) Rissanou, A.; Vlassopoulos, D.; Bitsanis, I. Thermal vitrification in suspensions of soft colloids: Molecular dynamics simulations and comparison with experiments. Phys. Rev. E 2005, 71, 011402. (29) Yu, H.-Y.; Srivastava, S.; Archer, L. A.; Koch, D. L. Structure factor of blends of solvent-free nanoparticle-organic hybrid materials: density-functional theory and small angle X-ray scattering. Soft Matter 2014, 10, 9120−9135. (30) Rallison, J. Brownian diffusion in concentrated suspensions of interacting particles. J. Fluid Mech. 1988, 186, 471−500. (31) Pusey, P.; Van Megen, W. Phase behavior of concentrated suspensions of nearly hard colloidal spheres. Nature 1986, 320, 340− 342. (32) Brambilla, G.; El Masri, D.; Pierno, M.; Berthier, L.; Cipelletti, L.; Petekidis, G.; Schofield, A. Probing the equilibrium dynamics of colloidal hard spheres above the mode-coupling glass transition. Phys. Rev. Lett. 2009, 102, 85703. (33) Doliwa, B.; Heuer, A. Cage effect, local anisotropies, and dynamic heterogeneities at the glass transition: A computer study of hard spheres. Phys. Rev. Lett. 1998, 80, 4915−4918. (34) Götze, W. Complex dynamics of glass-forming liquids: A modecoupling theory; Oxford University Press: USA, 2009; Vol. 143. (35) Allen, M.; Tildesley, D. Computer simulation of liquids; Clarendon Press: 1999. (36) Zimmerman, J.; WebbIII, E.; Hoyt, J.; Jones, R.; Klein, P.; Bammann, D. Calculation of stress in atomistic simulation. Modell. Simul. Mater. Sci. Eng. 2004, 12, S319. (37) Joe, B. Construction of three-dimensional Delaunay triangulations using local transformations. Comput.-Aided Geometric Des. 1991, 8, 123−142.

6841


Brownian Dynamics of a Suspension of Particles with Constrained

Recommend Documents