The flow properties of suspensions (e.g., colloids, ceramic slurries, and concrete) are of fundamental interest and play an important role in a wide variety of technological processes crucial to industry [Larson (1999)]. There have been many theoretical advances in understanding the rheological properties of simple suspensions (e.g., very dilute and semidilute suspensions, suspensions composed of spheroidal objects), however, understanding the flow of more complex suspensions (e.g., dense suspensions, random shaped particles, suspensions composed of particles that interact) remains a great challenge. Here, computational modeling can play an important role in investigating the properties of such systems. One possible approach is to apply standard computational fluid dynamics methods. This involves considerable effort in carrying out the difficult task of tracking boundaries between different fluid and solid phases, usually involving various meshing, moving grid, and interpolation schemes to account for motion of the rigid bodies. A second approach, based on the lattice Boltzmann method [Ladd (1997); Nguyen and Ladd (2002)], involves calculation of the momentum transfer which results from particles that "bounce" off the rigid body. The kinetics of the momentum transfer has to be carefully evaluated as the rigid body's surface can be located at any point and with any orientation between the lattice nodes from which the particles "propagate." A third, and perhaps best known approach, is called Stokesian dynamics (SD) [Brady and Bossis (1988); Phung and Brady (1996); Sierou and Brady (2001)]. In many respects, Stokesian dynamics serves as a standard benchmark as it was the first computational method to properly incorporate long range hydrodynamic interactions, Brownian forces, and lubrication forces for modeling suspensions composed of hard spheres. Some of its successes include the demonstration of shear induced ordering and shear thickening in dense hard sphere systems [ Foss and Brady (2000)]. Recently, a new computational method, called dissipative particle dynamics (DPD) [Hoogerbrugge and Koelman (1992); Koelman and Hoogerbrugge (1993)] has shown promise for modeling a variety of complex fluid systems. Further, DPD may potentially have some advantages over some computational fluids dynamics based approaches in that DPD can naturally accommodate many boundary conditions while not requiring meshing (or remeshing) of the computational domain. On the surface, DPD looks very much like a molecular dynamics algorithm where, in that case, particles subject to interatomic forces move according to Newton's laws. However, the particles in DPD are not atomistic but, more so, a mesoscopic representation of the fluid.
One can take several "philosophical" views of DPD. Ideally, one would like to think of DPD as a consequence of the systematic coarse graining of atomistic or microscopic domains. Indeed, there has been much effort in this direction [Flekkøy et al. (2000)]. While this view provides a general framework for understanding the structure of the DPD equations, there are still several gaps in bridging the microscopic and macroscopic domains. For example, it is necessary to impose constitutive relations (e.g., stress-strain rate relations) at some point. Hence, further work is needed to make such scale-up procedures clearer. A second view is that DPD belongs to a class of Lagrangian formulations of the Navier–Stokes equations {e.g., smooth particle hydrodynamics [Monaghan (1992)]}. Related but more sophisticated models [ Serrano and Español (2001); Espñol (1998)] utilize Voronoi cells to establish a grid that fills space or associate a time dependent volume parameter [Español and Serrano (1999)] to each DPD particle. From this perspective, DPD does not conserve volume in a proper sense, making the implementation of arbitrary equations of state difficult [Español and Serrano (1999); Español and Revenga (2003)]. In this respect, one can think of DPD as a "poor man's" Lagrangian formulation of the Navier–Stokes equations having sacrificed some rigor for computational expediency. For a third view, and what was probably the original intent, one can think of DPD as a somewhat abstract cellular-automata-based construct that, in certain regimes, recover hydrodynamics consistent with the Navier–Stokes equations {similar, in a way to how lattice gas and lattice Boltzmann methods were originally thought of [Rothman and Zaleski (1994)]}. Indeed, it has been shown, by mapping the DPD equations to an equivalent stochastic differential equation (the Fokker Planck equation) [Español and Warren (1995)] and applying a Chapman–Enskog analysis [Marsh et al. (1996, 1997)], DPD does produce hydrodynamic behavior consistent with the Navier–Stokes equations to second order in the Chapman–Enskog expansion. Thus, the challenge is to carefully connect solutions obtained from DPD to physical regimes of interest {cf. [Groot and Warren (1997); Dzwinel and Yuen (2000, 2002)]}. Hopefully, universal features of both the cellular automata approach and the "real" physical system can be exploited to help gain insight into the system of interest. Regardless of what computational approach or philosophical view one takes, it is extremely important to validate the computational method, especially if it is going to be used as a predictive computational tool.
Previous papers [Koelman and Hoogerbrugge (1993); Boek et al. (1997)] have demonstrated the potential of DPD to model colloidal suspensions including hard sphere and spheroidal objects [Boek et al. (1997)]. However, comparisons with experiments and theory have been more qualitative rather than quantitative and there were no comparisons with other approaches. In this paper, a DPD based approach for modeling suspensions is examined with an emphasis on comparing simulation results to well known theoretical predictions concerning simple flow scenarios and the rheology of dilute to semidilute suspensions. Tests include comparisons with the intrinsic viscosity prediction of Einstein and the Huggins coefficient for dilute and semidilute suspensions, respectively. At higher volume fractions, the DPD algorithm had to be modified to include lubrication forces as the usual DPD interactions are too weak to prevent overlap of the rigid bodies. Results are compared with previous studies concerning the flow of dense suspensions based on the Stokesian dynamics method and experimental data. As an alternative to the commonly used Lees–Edwards boundary condition [ Allen and Tildesley (1987)], which can, roughly be thought of as a constant applied strain rate, simulations were also carried out using a constant applied stress. It was found that use of a constant stress to drive the system helped mitigate large temporal fluctuations in the derived viscosity which occurred in the constant strain rate case. Interestingly, rheological measurements at higher volume fractions are often carried out using a constant applied stress. Spheroidal rigid bodies are also considered. The rotation period of a single prolate spheroid under shear is consistent the predictions of Jeffery (1922); Eirich (1967). Studies of the flow of spheroids at higher volume fractions produce an apparent nematic phase [Larson (1999)]. An example of application of the DPD algorithm to model flow in other geometries like that encountered in the flow and placement of concrete is given. To contrast this work with previous DPD based simulations of suspensions [Boek et al. (1997)], it should be noted that in this paper lubrication forces are explicitly included in the simulations of dense sphere suspensions. In addition, size polydispersivity, Jeffery's orbits and the onset of an apparent nematic phase were studied. Finally, flow in alternate geometries and under applied stress instead of applied strain were examined.