At higher volume fractions (φ > 0.4), it became increasingly difficult to carry out simulations without having sphere overlaps occurring. This problem worsened as the Peclet number increased. One reason for overlaps is that the interactions between individual DPD particles are "soft" allowing for some penetration. A simple fix to the code was attempted by including a very steep repulsive interaction between spheres. While such forces greatly suppress the overlaps, it was found that the relative viscosities were, at high Pe ≈ 10 000, low by a factor of 2 or more when compared to SD or experimental data (Fig. 4). Indeed, these simulation results are roughly akin to the extrapolation of relative viscosity data, from the low to high Pe number limit, without consideration of lubrication forces. From lubrication theory [Kim and Karrila (1991)], it is well known that the force between approaching spheres, to lowest order, scales as (VA–VB)/sA B where VA and VB are the velocities of spheres labeled A and B andsAB is the distance between the nearest points of the respective two sphere surfaces {the reader is referred to the literature for further details on lubrication forces [Kim and Karrila (1991)]}. Clearly, as smaller and smaller distances between spheres are probed, lubrication forces are not properly accounted for by the usual DPD interactions, in part because the spatial resolution required is impractical. Hence, it was necessary to directly incorporate the lubrication forces into the DPD code. Here, analytical expressions for the lubrication forces were kept up to first order, (including terms that scale as 1/sAB, ln sAB, and sAB ln sAB). Unfortunately, lubrication theory makes the assumption that the distance between spheres is much smaller than the radius a, so it is not precisely clear when to turn off the DPD interactions between spheres and when to turn on the lubrication forces. For simulation results presented in this paper, only spheres where sAB<a were evaluated for lubrication forces. Also, the velocity dependent DPD interparticle interactions between spheres were turned off and an empirical functionSf was introduced to smoothly incorporate the lubrication forces into the algorithm. For this study, the lubrication forces were multiplied by the following function
Fig. 4. Comparison of simulation predictions and experimental data for high volume fraction effective hard sphere systems. Here φ = 0.49. The open circles are results from Stokesian dynamics simulations [Foss and Brady (2000)], open triangles are from DPD simulations without lubrication forces and open squares are from DPD simulations with lubrication forces. The open diamonds are results from the DPD simulation where a constant stress was applied instead of the Lees–Edwards boundary condition. The solid and dashed lines are experimental data from sheared suspensions of silica particles [Bender and Wagner (1996)].
While a "best" choice of smoothing function still needs further study, the
form chosen was fairly simple and allows for a close approach to the "true"
lubrication force when sAB 
a. Indeed, some other forms of
smoothing functions were tested, but there was no significant difference in the results.
For such dense systems, all neighboring spheres nearly touch and, as a result, the force
between them is dominated by its singular nature.
Even with the introduction of the lubrication forces it was difficult, at values of Pe ≈ 1000 and greater to avoid some overlap. This was a result of using a constant time step that was not sufficiently small to account for forces when the spheres were in very close proximity to each other. This issue has been noted elsewhere in the literature [Ball and Melrose (1995a, 1995b)]. It is interesting that earlier simulations using SD [Foss and Brady (2000)], with a constant time, allowed for overlaps of about 1% of the sphere radius. When this occurred, a very small separation was assumed (of order 10–8 the sphere radius) and the simulation was allowed to progress. This approximation was probably not unreasonable because, at the length scales probed, the forces between the spheres cannot be described by lubrication theory alone. Also, as the spheres approach each other, a slip velocity may become apparent since the mean free path of the fluid atoms will be of order the spacing between sphere surfaces. Consequently, an assumption underlying the derivation of lubrication forces, no slip at fluid/surface boundary, would need to be modified and there may no longer be a singularity in the force as the spheres touch for actual physical systems. Further research is needed on this issue.
An attempt was made to avoid overlaps by including a dispersive short range interaction potential with an adjustable decay width, λ. Here the hope was that introducing a repulsive force would disperse the spheres enough to avoid overlaps from taking place. Then, by decreasing the decay width, we could probe smaller and smaller distances between spheres to see the effect of the lubrication forces. The following form of a repulsive force, similar to the construction used by Foss and Brady (2000), was chosen:
where Z is a constant, λ is the decay width, and
AB
is a unit vector pointing from the center of sphere A to sphere B. Not
unexpectedly, it was found that the viscosity was sensitive to λ and increased with
decreasing λ. For one set of simulations, 27 spheres were used with φ = 0.477
and Pe ≈ 1000. Here the suspension was sheared with a strain equivalent to over
ten simulation cells for cases of λa ≈ 8.0×10–5,
2.0×10–5, and 4.0×10–6. The relative
viscosity was 8.58, 10.1, and 11.1, respectively. In this case, no overlaps occurred as
the spheres managed to squeeze by each other, although coming quite close, with the ratio
of the distance between sphere centers to diameter equal to 1.000 000 01 and smaller.
However a second set of simulations was carried out with the sphere radius about 1%
larger, making φ = 0.49. The same set of λa was used. This time, despite
the relatively small increase in sphere radius, the time step needed for the simulation to
proceed without overlaps, was too small to be practical.
Instead of relying on a dispersive force to help separate the spheres, it was then decided that the best route would be to incorporate a variable time step into the code. A simple modification was made to the algorithm such that, as the spheres approached each other, the time step was reduced by a factor of 5 if the sphere's projected trajectory appeared close to creating an overlap. A suspension of 663 spheres was simulated with volume fraction 0.49. It was found that at Pe = 1000 the system evolved in a relatively smooth fashion but at higher Pe ≈ 10 000 large fluctuations were found in the stress. As the viscosity is related to stress [see Eq. 13] it would also appear as if the viscosity was dramatically fluctuating (see Fig. 5). Recently, in experimental studies of constant strain rate driven dense suspensions [Lootens et al. (2003)], large fluctuations in the stress have been observed. The large fluctuations have been related to the onset of a jamming transition. What was not clear from the simulations was whether the onset of the larger fluctuations in the measured viscosity was a consequence of the constant strain rate boundary condition when employing the Lees–Edwards boundary condition. As an alternative to the Lees–Edwards boundary condition, a simulation was set up so that an applied stress was used to drive the system. Here two narrow bands of spheres were constrained to move in parallel planes having a spacing of about four sphere diameters between each other. A force was applied on the spheres, in opposite directions in each separate plane, so that a shearing motion was established. The simulation cell contained 340 spheres. Figure 6 shows the average velocity difference ΔV between the top and bottom bands of spheres as a function of integrated strain. In this case, the resulting strain rate is no longer constant with the average velocity varying about 5%–10%. Clearly, temporal fluctuations in the measured viscosity were greatly reduced for the constant stress case (see Fig. 5). On the other hand, the average viscosity determined from the stress controlled simulation was about 10%–30% higher than the strain controlled simulation in this high Pe regime. Since the gap between plates was about four sphere diameters, finite size effects could have made the relative viscosity appear higher. A related observation was made by Boek and van der Schoot (1998) concerning finite size and resolution effects. Here it was found that, at low Pe, estimates of the relative viscosity improved (when compared to experimental values) when the colloidal spheres were made sufficiently small. Unfortunately, it was not clear if this was a consequence of having smaller spheres relative to the simulation box size alone (a finite size effect), or, in part,the result of a repulsive force between particles helping better disperse the spheres and reduce overlaps (a resolution effect). Such finite size effects will be the subject of future research.
Figure 5. Calculated values of relative viscosity as a function of integrated strain rate. The dashed line corresponds to data from a constant stress driven system. The solid line is from a simulation with constant strain rate (Lees–Edwards boundary condition). Note the large temporal fluctuations in relative viscosity for the constant strain rate case as spheres must respond to an unyielding motion resulting from such boundary condition.
Figure 6. Constant stress driven shear. The velocity difference of the two parallel regions where the force is applied is given by ΔV. To set the velocity scale, ΔV = 40 corresponds to Pe ≈ 10 000. For this system φ ≈ 0.49. There were moderate fluctuations in ΔV as the simulation progressed.
For Pe = 10 000 the relative viscosity was, for the constant strain rate case, about 10% higher than that of previous Stokesian dynamics simulations (see Fig. 4). Although within the statistical uncertainty, it is possible that the correction for the overlaps in the SD simulations was in part the cause of the discrepancy as can be seen from the above study of dispersed spheres. Allowing for the smaller distances between spheres (i.e., as λa is decreased) would have probably increased the viscosity in the high Pe case. Finally, using a variable time step does not guarantee, for some unique configuration of spheres and flow history, that the time step may again become too small to be practical. However, this problem is not necessarily the same as jamming, where the system cannot move without overlaps, but more a case of developing a reasonable strategy for updating the sphere positions.
Figure 4 contains experimental data from a stress controlled measurement of a sheared silica suspension due to Bender and Wagner (1996). The agreement with our stress controlled simulation is good in regards to capturing trends. However, one needs to take great care when comparing simulations of such an idealized system to experiments (and vice versa). Interestingly, the silica particles were slightly polydisperse so that one might think the experimental measurements of viscosity would be a bit low, perhaps up to 10% or so at this volume fraction, relative to a monosize sphere case (see section on polydispersivity). Unfortunately, such corrections for polydispersivity would make agreement worse. Second, the data shown and other comparable experimental data {e.g., spherical silica particles [Bender and Wagner (1996)] and poly(methylmethanylate) (PMMA) [D'Hane et al. (1993); Phan et al. (1996)]} are based on measurements of suspension with particles approximately 100–1000 nm in diameter. Again, consider the earlier set of simulations where a dispersive force was introduced. As the width of the potential λa ranged from 0.0001 to 0.000 004 the viscosity was not quite at its asymptotic limit (also note Pe ≈ 1000). Probing the experimental particles at similar λa (and at higher Pe) would put one at atomic scales and smaller. Clearly at such length scales the silica and PMMA particles are not exactly hard spheres and the embedding fluid can no longer be represented as a continuum. So it is not clear if the experiments in the high Pe > 1000 can be modeled as hard sphere fluids without consideration of these features. A related observation was made by Ball and Melrose (1995b) who described their simulation results as unphysical with respect to modeling colloidal systems when such small distances were probed.