Next: Integration of the equations of motion Up: Main Previous: Introduction

Basic DPD Equations

I start by briefly reviewing the basic equations of DPD. In DPD, as in molecular dynamics, the evolution of the position, r_i, and momentum, p_i > m v_i, of particle i with mass, m, and velocity v are described by

where F_ij is the force on particle i due to particle j and the dot indicates a time derivative. Interparticle forces are typically represented as three types: conservative F, dissipative F, and random F so that

The conservative force is simply a central force, derivable from some effective potential φ_ij. The dissipative force is proportional to the difference of velocity, v_ij = v_i–v _j, between particles and acts to slow down their relative motion, producing a viscous effect. The random force (usually based on a Gaussian random noise) helps maintain the temperature ofthe system and provides an additional viscous effect. The three forces are given below

The distance between the DPD particles i and j is given by r_ij, ê_ij is a unit vector pointing from particle j to particle i, w_R(r _ij) and w_D(r_ij ) are weight functions and _ij is a randomly fluctuating variable described by Gaussian statistics. It can be shown that, in order to maintain a well defined temperature by way of consistency with the fluctuation-dissipation theorem [Español and Warren (1995)], coefficients describing the strength of the dissipative (γ) and random (σ) forces must be coupled, that is

where k_b is the Boltzmann constant and T is the temperature. Further, so that the DPD fluid system possess a Gibbs–Boltzmann equilibrium state, the following relation must hold (detailed balance for an infinitesimal time step) [Español and Warren (1995)]:

In this study, the choice of parameters and weight functions closely follow that described in Groot and Warren (1997). Here, w_R (r_ij) = 1–r_ij for (r_ij<1) and w_R(r _ij) = 0 for r_ij ≥ 1. All lengths described in this paper are defined in units of the cutoff radius, r_c = 1, of the DPD interaction. The conservative force is taken to be F = F_m (1–r_ij)ê_ij. For all the simulations in this paper, σ ≥ 40 and F_m = 75k_b T/ρ where ρ is the global density of DPD particles. Units of k_bT are chosen such that k_b T = 1 and F_m was chosen so that the DPD fluid has the same compressibility of water [see the discussion in Groot and Warren (1997)].

An important parameter that characterizes suspensions under shear is the Peclet number, Pe. Peclet number is a dimensionless number describing the competition between viscous and Brownian forces and, for spheres, is given by Pe = 6πµa³k_bT. Here, µ is the viscosity, a is the sphere radius, and is the shear rate. Also, for spheres under shear, the Reynolds number is given by Re = ρa²/µ. In general Re ≈ O(1) or smaller in this study. Depending on the simulation, system sizes were 15³, 45³, and 90³ in our units. Finally, because the DPD interactions are short range the code parallelized in a fairly efficient manner. For example, a spatial decomposition version of our code scaled nearly linearly up to about 16 processors on a Linux cluster. For more information on the parallelization of this code see Sims and Martys (2004).