Understanding the flow properties of complex fluids like suspensions (e.g., colloids, ceramic slurries, and concrete) is of importance to industry and presents a significant theoretical challenge. The computational modeling of such systems is also a great challenge because it is difficult to track boundaries between different fluid/fluid and fluid/solid phases. Recently, a new computational method called dissipative particle dynamics (DPD) [2] has been introduced which has several advantages over traditional computational dynamics methods while naturally accommodating such boundary conditions. In structure, a DPD algorithm looks much like molecular dynamics (MD), where atomistic particles move according to Newton's laws. However, the DPD "particles" are a mesoscopic description of the fluid, and do not represent individual atoms or molecules, but loosely correspond to "lumps" of fluid or clusters of molecules. As a result, the interactions between the DPD particles are not directly based on a Lennard-Jones potential, but are typically subject to three types of forces, namely, conservative forces, dissipative forces, and a random force. All of the forces conserve momentum and mass. The conservative force is simply a central force, derivable from some potential. The dissipative force is proportional to the difference in velocity between particles and acts to slow down their relative motion. The dissipative force can be shown to produce a viscous effect. The random force (usually based on a Gaussian random noise) helps maintain the temperature of the system while producing a viscous effect. It can be shown that, in order to maintain a well defined temperature by way of consistency with a fluctuation-dissipation theorem [3], coefficients describing the strength of the dissipative and random forces must be coupled. By mapping of the DPD equations of motion to the Fokker-Planck equation [4], it has been demonstrated that the DPD equations can recover hydrodynamic behavior consistent with the Navier-Stokes equations.
As in MD, the forces on each particle are computed in each time step. The particles are then moved and the forces recomputed. In DPD the interparticle interactions are chosen to allow for much larger time steps so that physical behavior, on time scales many orders of magnitude greater than that possible with MD, may be studied. The original DPD algorithm [2] used an Euler algorithm for updating the positions of the free particles (which represent "lumps" of fluids), and a leap frog algorithm for updating the positions of solid inclusions (rigid bodies). Our algorithm QDPD [5], for quarternion- based dissipative particle dynamics, is a modification of DPD that uses the velocity-Verlet algorithm of Groot and Warren [6] to update the positions of both the free particles and the solid inclusions. The velocity-Verlet algorithm for DPD [5] is chosen because it is less sensitive to variation in time step size than the Euler algorithm. The solid inclusion motion is determined from the quaternion-based scheme of Omelayan [7] (hence the Q in QDPD).
QDPD in its present form is being used to study the steady-shear viscosity of a suspension of solid inclusions (such as ellipsoids) in a Newtonian fluid. The model consists of N particles moving in a continuum domain of volume V. As in MD the system is completely defined by specifying all N positions ri and momenta pi (i = 1, ..., N). To model a rigid body inclusion in a fluid, a subset of the DPD particles are initially assigned a location in space so that they approximate the shape of the object [8]. The motion of these particles is then constrained so that their relative positions never change. The total force and torque are determined from the DPD particle interactions and the rigid body moves according to the Euler equations. As mentioned above, our simulations use a quaternion-based scheme developed by Omelayan and modified by Martys and Mountain [5] for a velocity-Verlet algorithm to integrate the equations of motion. Finally, we use a Lees-Edwards boundary condition [9] (pp 246-247) to produce a shearing effect akin to an applied strain at the boundaries.
The basic idea is to compute all of the forces on each particle (which accounts for the momenta change in the collision phase) during each time step, and then move the particles (propagation phase). The forces are short-range and are a sum of contributions over pairs of particles. The interaction decays rapidly with separation, which means that only particles closer than some cutoff distance rc need be considered. Several methods are available for identifying the nearest neighbors of a particle, i.e, those within the cutoff distance. QDPD uses an implementation of the link-cell method of Quentrec et al. [10] described in Allen and Tildesley’s book [9] (pp. 149-152). Here, the simulation box is partitioned into a number of cells. For example, see Fig. 1, which depicts a 2D system. To find the particles within the cutoff distance rc of the particle shown in the central cell, it is sufficient to only consider particles within the central cell and each of its eight nearest neighbor cells (where rc is ≤ the cell widths in X and Y, lx and ly). The use of Newton's third law makes it possible for us to only have to consider half of the nearest neighboring cells, which are cross-hatched (lines parallel to the right-diagonal in the cell) in the figure. Generalizing this to all particles in the system, a linked list of all the particles contained in each cell is constructed every timestep. Then, for each particle, the selection of all particles within the cutoff is achieved by looping over one half (considering Newton’s third law) of all nearest neighbor cells, and considering only the particles within these cells. We show this schematically in Fig. 2.
Fig. 1. Schematic diagram of link-cell algorithm for a two dimensional system (after Tildesley, Pinches, and Smith [11]).
Fig. 2. Schematic 2-D representation of the link-cell algorithm.
The (forces calculation) search scheme involves an outer loop over all 25 link-cells. In this outer loop, each particle in a link-cell interacts with all particles within its link-cell that are within rc of the particle. Then there is an inner loop over four of the eight nearest neighbor link-cells, and each particle interacts with all of the particles within the chosen neighbor link-cells that are within rc of the particle. For example, particles in cell 13 interact with other particles in 13 plus particles in 17, 18, 19, and 14 that are within the cutoff distance of the chosen particle. Note that to account for the forces on particles in edge cells, periodic boundaries are used to have, for example, 25 interacting with the appropriate nearest neighbor periodic cells 4', 5', 1", and 21' (more on this later). The program for figuring out nearest neighbor cells is easy to set up. Introducing cell indices Ix and Iy for the 2D grid in Fig. 2, each cell's index in the 2D grid can be computed from
where MOD is the function which returns the modulo of its arguments and Mx and My are the number of cells in X and Y (Ix = {1, Mx}, Iy = {1, My}). For each cell, one-half of the nearest neighbor cells are given by
which correctly gives the cell neighbors of 13 to be 17, 18, 19, and 14. Now we can explain the treatment of Newton's third law. A particle in cell 13 interacts with the particles in 8 neighboring cells, but the algorithm only checks particles in 17, 18, 19, and 14. Interactions of particles in 13 with particles in 12, for example, are treated when particles in cell 12 are the focus of attention, and similarly for 7, 8, and 9. Note that this formula also gives the nearest neighbors of particles in cell 25 to be those in 4, 5, 1, and 21 (not the periodic cells 4', 5', 1", and 21' ). This will be explained later. Because of their regular arrangement, the list of neighboring cells is fixed and may be precomputed once and for all at the beginning of the program (in subroutine MAPS).