# Direct numerical simulations for non-Newtonian rheology of concentrated particle dispersions

## Abstract

The non-Newtonian behavior of a monodisperse concentrated dispersion of spherical particles was investigated using a direct numerical simulation method, that takes into account hydrodynamic interactions and thermal fluctuations accurately. Simulations were performed under steady shear flow with periodic boundary conditions in the three directions. The apparent shear viscosity of the dispersions was calculated at volume fractions ranging from 0.31 to 0.56. Shear-thinning behavior was clearly observed at high volume fractions. The low- and high-limiting viscosities were then estimated from the apparent viscosity by fitting these data into a semi-empirical formula. Furthermore, the short-time motions were examined for Brownian particles fluctuating in concentrated dispersions, for which the fluid inertia plays an important role. The mean square displacement was monitored in the vorticity direction at several different Peclet numbers and volume fractions so that the particle diffusion coefficient is determined from the long-time behavior of the mean square displacement. Finally, the relationship between the non-Newtonian viscosity of the dispersions and the structural relaxation of the dispersed Brownian particles is examined.

###### pacs:

82.70.-y, 82.20.Wt, 47.50.-d, 83.50.Ax## I Introduction

The links between macroscopic rheological properties and microstructures in colloidal dispersions have been extensively investigated for many systems, including dispersions of sterically-stabilized and charged-stabilized particles in host fluids (1); (2).

The behavior of a monodisperse dispersion composed of solid particles immersed in a Newtonian host fluid strongly depends on the volume fraction of the dispersed particles and the shear rate . When the shear rate is zero (), the shear viscosity of the dispersion is referred to as the zero-shear viscosity . In the dilute limit (), the zero shear viscosity is well approximated by Einstein’s formula (3):

(1) |

where is the shear viscosity of the host fluid. In a concentrated dispersion, theoretical difficulties become rather severe, since the behavior of the dispersed particles is complicated by the interactions between the particles and thermal fluctuations. In particular, the solvent-mediated many-body hydrodynamic interactions (HI) between the particles complicate the dynamical behavior. A number of experiments for concentrated dispersions have been performed to reveal the origin of the non-Newtonian behavior of these dispersions, and several semi-empirical formulas for have been proposed to characterize the experimental results. For example, the Krieger-Dougherty relationship (4):

(2) |

where is the intrinsic viscosity and is the packing volume fraction at which the viscosity diverges, is often used for fitting the experimental data for uniform colloidal spheres suspended in non-aqueous media.

When a dispersion is subjected to shear (), the flow properties of the dispersion show a variety of non-Newtonian behaviors, such as shear-thinning and shear-thickening. These non-Newtonian behaviors are associated with the changing microstructures of the dispersion. Several physical mechanisms for these peculiar behaviors have been proposed; for example, shear-induced order-disorder transitions (5); (6); (7), formations of dynamic clusters of the particles (8); (9). However, a full understanding of the relationships between the rheological properties and the microstructure has not yet been obtained, despite extensive studies.

Computer simulations are very powerful tools in the direct investigation of the dynamics of individual particles in concentrated dispersions. The Stokesian dynamics (SD) method (10) has been widely used to measure the rheology of dispersions and provides valuable information regarding the non-Newtonian behavior of flowing dispersions (14); (15); (17); (16). The SD method, however, is based on the Stokes approximation, which assumes that all relaxation times associated with fluid motions are short as compared with those of the particle, i.e., , where and . Here, and are the density of the fluid and the particle, respectively, and is the radius of the particle. With this assumption, the HI is treated as the Ronte-Prager-Yamakawa (RPY) tensor and the lubrication correction. Furthermore, it is assumed that the relaxation time associated with the particle’s inertia is zero (). These approximations are valid for the motion of the particles at time scales much greater than the time scales of the relaxation of both the fluid and the particle inertia. Therefore, the short-time motion of Brownian particles over the kinematic time scale cannot be described by simulation methods based on the Stokesian approximation. On the kinematic time scale, the dynamic coupling between the fluid motion and the particle motion remains strongly.

Sheared dispersions have a non-dimensional parameter that includes the particle’s inertia, e.g., the particle Reynolds number , where represents the characteristic time due to shear. For most colloidal dispersions, the particle Reynolds number is very small. The SD method has a particle Reynolds number of zero and, therefore, cannot be applied to problems with finite particle Reynolds numbers. Typical examples include the motions of dispersions subjected to strong shear or composed of large particles; in these cases, the particle Reynolds number has a relatively large finite value.

Another problem that cannot be treated by the SD method is the short-time motion of the particle, in which the fluid inertia becomes significant even if the particle Reynolds number is very small. For example, the characteristic time scale is for a neutrally buoyant particle of 1 in water, and for . The effects of the fluid inertia appear as memory effects of the particles. For a complete understanding of these flowing dispersions, full time-dependent HIs are required.

In recent years, several numerical methods have been developed in order to accurately simulate dispersions in a variety of situations, including those described above. These methods of dispersion modeling are based on the same approach, which involves resolving the fluid motion simultaneously with the particle motion. We refer to this approach as the direct numerical simulation (DNS) approach. This approach enables us to accurately treat the full time-dependent HI. The numerical methods differ mainly in the approach used to resolve the HIs between the fluid and particle motions (18); (19); (20); (21); (22); (23); (24); (25); (26).

In this work, we apply a direct numerical scheme based on the smoothed profile method (SPM) (26); (27); (28); (29); (30) to a monodisperse concentrated dispersion of repulsive and neutrally buoyant Brownian particles in a shear flow. In the SPM, the Navier-Stokes equation for the fluid motion is discretized on a regular grid, and the Newtonian equations for the particle motion are solved simultaneously with the fluid motion. This developed scheme accounts for thermal fluctuations, a shear flow, and memory effects. The SPM can reproduce the correct short-time behavior of a single Brownian particle in a shear flow on the kinematic time scale, and the numerical tests for a single Brownian particle in a shear flow were reported in (30).

Although some groups have studied flowing dispersions composed of repulsive spherical particles by using a DNS method based on the lattice Boltzmann method (31); (32); (33), all the simulations have been applied to a dispersion of non-Brownian particles and ignore the thermal fluctuations of the particles. Systematic analyses have not yet fully carried out for concentrated dispersions of Brownian particles in a shear flow at finite Reynolds numbers.

In the present study, we examine the non-Newtonian rheology of concentrated dispersions in shear flow at finite Reynolds and Peclet numbers using a DNS method that takes into account hydrodynamic interactions and thermal fluctuations accurately. Note that the simulated situations are different from those of the prevailing experiments for colloidal dispersions, since even the lowest particle Reynolds number in these simulations is several hundred times larger than those of the experiments. We first present the simulation method and the manner in which the apparent shear viscosity of the dispersions is calculated. Three-dimensional simulations are performed with periodic boundary conditions, and the non-Newtonian behavior of the shear viscosity is obtained for several volume fractions. Both the high and low shear limiting viscosities obtained from the simulations are then compared with the Krieger-Doughty relationship. Moreover, we investigate the short-time motions of Brownian particles in a sheared concentrated dispersion on the kinematic time scales. The mean square displacement (MSD) of Brownian particles is monitored at several Peclet numbers and volume fractions, and the long-time diffusion coefficient is determined from the long-time behavior of the MSD. We also suggest a simple relationship between the non-Newtonian viscosity of the dispersions and the structural relaxation time of the dispersed Brownian particles.

## Ii Simulation Method

A direct numerical scheme that implements both a shear flow and thermal fluctuations is briefly explained. A more detailed explanation is given in a previous publication (30). We consider a monodisperse dispersion of repulsive spherical particles of diameter in a Newtonian host fluid. The particles interact via a truncated Lennard-Jones (LJ) potential:

(3) |

where is the distance between two particles and the parameter characterizes the interaction strength. The position of the th dispersed particle is , the translational velocity is , and the rotational velocity is . The time evolution of the th particle with mass and moment of inertia is governed by Newton’s equations of motion:

(4) | ||||

(5) |

where and are the hydrodynamic forces and torques exerted by the host fluid on the particle. is a repulsive force arising from the potential of Eq.(3), which prevents the particles from overlapping. and are random forces and torques, respectively, due to thermal fluctuations. These random fluctuations are assumed to be Markovian (white or time delta-correlated) and determine the particle temperature . The procedure for determining the temperature is described in (28); (29); (30).

In the SP method, the velocity and pressure fields, and , are defined on three-dimensional Cartesian grids, which consist of fluid and particle domains. In order to distinguish the particle and fluid domains on the grids, a smoothed function , which is equal to in the particle domains and in the fluid domains, is introduced. These domains are separated by a thin interfacial domain of thickness . The system size is .

The time evolution of the velocity field is governed by the Navier-Stokes equation with the incompressibility condition :

(6) |

where the stress tensor , and is an external force field that is introduced to enforce a simple shear flow on the system. The flow is imposed in the direction and the external force is introduced as a constraint force so that the velocity field satisfies at and at , where denotes the distance in the velocity gradient direction. A simple shear flow with a shear rate of is then approximately produced over a range from to . represents a body force that ensures the rigidity of the particles and the appropriate non-slip boundary conditions at the fluid/particle interface, which is further elaborated upon in reference (26); (27).

The unit of length is taken to be the lattice spacing , and the unit of time is . Unless otherwise stated, we set , , , , , , , and . Assuming dispersions of neutrally buoyant particles of radius m in water at room temperature, our unit length and time correspond to m and s, respectively The second-order Runge-Kutta algorithm is used to integrate the Newtonian equations. The Navier-Stokes equation is discretized with a Fourier spectral scheme in space and with a second-order Runge-Kutta scheme in time. The discretized time step is . Although this value is choosen from the stability condition of the Navier-Stokes equation, it can be used safely for the particle’s equations motion because is much smaller than the Lennard-Jones time unit . The phase diagram of the present (36:18) LJ system depends weakly on the system temperature as well as the volume fraction of the particles. To avoid the delicate issue of crystallization, all the simulations in the present paper were carried out at and with which the system is in amorphous states.

To measure the rheological properties of the particle dispersion in shear flow, we can calculate the apparent stress of the dispersions in the following manner. The momentum equation for the dispersion is formally written as:

(7) | ||||

(8) |

where denotes the stress tensor of the dispersion including the inertia, the pressure and the viscous terms. The full stress tensor of the flowing dispersion is then defined by introducing a convective momentum-flux tensor explicitly as:

(9) |

where represents momentum transport by the bulk flow of the dispersion.

Although cannot be calculated directly, we can obtain the apparent stress by using the local stress :

(10) | ||||

(11) | ||||

(12) | ||||

(13) | ||||

(14) |

with a volume . In the derivation of Eq.(11), we use an second rank identity, . In the steady state,

(15) |

where denotes time averaging over the steady state. The time-averaged apparent shear stress of the dispersion can then be written as:

(16) |

Then one can obtain the apparent shear stress under steady shear flow from the external force imposed in the Navire-Stokes equation.

## Iii Results and discussion

Simulations were performed in a three-dimensional cubic box, whose side length is , with periodic boundary conditions. Most present simulations were done with , where the number of particles are , , , , and for , , , , and , respectively. The index of axis , , and represent the flow, velocity gradient, and vorticity directions, respectively. For a spherical particle, , . The temperature was determined by equilibrium calculations before the shear flow was imposed. The temperature is . The initial configuration of the particles is set to be a random distribution. A large number of simulations were performed for volume fractions of and shear rates of . These systems have particle Reynolds numbers ranging from to 1.6.

For dispersions composed of Brownian particles in the steady state, the apparent shear viscosity of the dispersion is defined as:

(17) |

where denotes the -components of . To examine the system size effects in the present simulations, we calculated the apparent shear viscosity for three different system sizes for a constant , i.e., with , with , and with . We found only negligibly small differences among the values of the apparent shear viscosity calculated for the three systems within a statistical error.

Figure 1 shows the dependence of the apparent shear viscosity on the Peclet number for several volume fractions. Here, the Peclet number is defined as . For the lowest concentration (), the viscosities remain nearly constant and exhibit Newtonian behavior. For higher concentrations (), the dispersions show non-Newtonian behavior. As increases, the shear-thinning behavior is clearly observed from the higher plateau region for of order to a lower plateau region for of about 10. From both of the plateau values of the viscosity curve, we can obtain the low shear limiting viscosity (identified as ) and the high shear limiting viscosity for each volume fraction.

In order to evaluate and , we fit our simulation data into the following simple empirical function of and :

(18) |

where is a fitting parameter. This empirical function is plotted in Fig. 1 with solid lines for . One finds that it coincides with the simulation data reasonably well. It is also seen that the onset of shear thinning appears at smaller Peclet numbers with increasing volume fraction. We note that the particles are randomly distributed in the dispersions in the present simulations, i.e. all the viscosity data are taken in situations at which no shear-induced crystallization occurs.

Figure 2 shows the dependence of the low shear limiting viscosity and the high shear limiting viscosity on the volume fraction. Both shear limiting viscosities increase monotonically with the volume fraction. The simulation results for both of the shear limiting viscosities agree well with the semi-empirical relations of Eq.(2) by Krieger-Dougherty, where for and for . The dotted line represents Einstein’ s formula of Eq.(1).

The fluid and particle inertia contribution to the apparent shear stress can be written as

(19) |

which is similar to the Reynolds stress of a uniform fluid in turbulence. This stress represents the strength of the hydrodynamic instability. The inertia contribution to the viscosity, , is about two order of magnitude smaller than the apparent shear viscosity. We conclude that the inertia contribution is negligible in the range .

We next examine the dynamical motion of Brownian particles in a shear flow. We analyze the mean square displacement (MSD) in the vorticity direction () for the Brownian particles,

(20) |

Figure 3 (a) shows the MSD of the dispersion with the lowest Peclet number for each volume fraction. These results are considered to be the MSD at thermal equilibrium. As the volume fraction is increased, the dynamical behavior of the MSD varies greatly, as compared with that of the hydrodynamic analytical solution for a single Brownian particle in a shear flow (34). The analytical solution is derived from the generalized Langevin equation with memory effects. With increasing volume fraction, MSD increases more slowly with time. For , we found that a plateau region starts to appear around a time scale of order , which is slightly greater than the kinematic time . This is a typical behavior of colloidal dispersions in glassy states, e.g., colloid glasses.

Figure 3 (b) shows the MSD of the dispersion at the highest Peclet number for each volume fraction. Since here, the shear force is much stronger than the thermal force. The volume fraction dependence of the MSD in this figure is opposite to the previous case at the lowest Peclet number shown in Fig. 3 (a). At short times, the MSD grows more rapidly in time with increasing volume fraction because HIs between particles are more enhanced at higher volume fraction. At later times, the MSDs tend to exhibit diffusive motions with a volume fraction independent diffusion coefficient.

We then calculated the long-time diffusion coefficient in the vorticity direction via

(21) |

and examined the system size effects in . In contrast to the negligibly small system size effects observed in the apparent shear viscosity , shows notable increase with increasing system size for all volume fractions. It is, however, confirmed that the diffusion coefficients behave as

(22) |

with , where represents the value of diffusion coefficient extrapolated for . The system size effects observed in the present simulations are essentially the same as those reported in earlier papers (11); (12); (13). We calculated from extrapolation of our simulation data with finite . Figure 4 shows the volume fraction and the Peclet number dependencies of normalized by . For the low Peclet numbers, the diffusion coefficient decreases with increasing volume fraction, and this dependence is similar to that of the dispersion at thermal equilibrium. On the other hand, for the highest Peclet number, the volume fraction dependence of is almost constant. We can see that the diffusion coefficients at high volume fractions are strongly affected by shear. The dynamics of the Brownian particles differs considerably, depending on whether the thermal or shear force is dominant. The volume fraction dependence is in qualitative agreement with the numerical results achieved based on the SD method, which is valid for long time scales (14).

Finally, the relationship between the non-Newtonian viscosity of the dispersions and the structural relaxation of the dispersed Brownian particles is examined. The structural relaxation time is defined as

(23) |

where is the diffusion coefficient for each volume fraction at zero shear. This relaxation time represents a time needed for a particle to diffuse away a distance comparable to its radius. Figure 5 shows versus for . Here, we used at the lowest Peclet number as , and can be understood as a reduced Peclet number . One can see that the data lies on a single master curve fairly well, where and is a fitting parameter. The shear-thinning starts at around . This behavior is very similar to the non-Newton rheology of supercooled liquids (35).

We confirmed that the present numerical method of introducing the thermal fluctuation successfully reproduces the fluctuation dissipation theorem for time scales longer than the so-called the Brownian time (28); (29). Although it was confirmed also that the present method works quite well for the volume fractions considered in the present study, further careful tests must be needed for highly concentrated dispersions, where the fluctuation of a tagged particle tends to correlate with motions of surrounding particles.

## Iv Conclusions

We investigated the rheological properties of monodisperse concentrated dispersions with repulsive spherical particles by using a DNS method that accounted for a shear flow and thermal fluctuations. Three-dimensional simulations were performed at Peclet numbers ranging from to and at particle Reynolds numbers from to . The apparent viscosity for is almost constant over the Peclet number change. For , the viscosities decrease from the plateau region at the low Peclet number to the plateau region at the high Peclet number. From the viscosity versus Peclet number curves, we can obtain both the low and high limiting viscosities, and these results are in good agreement with the Krieger-Doughty relationship. The inertia contribution to the apparent shear viscosity is very small throughout the entire range of Peclet numbers and volume fractions examined in the present study. As the volume fraction increases, the behavior of the MSD in the vorticity direction deviates from the analytical solution for a single Brownian particle in a shear flow. For the lowest Peclet number, the MSD develops more slowly in time with increasing volume fraction. At , an onset of glassy dynamics was observed. On the other hand, for the highest Peclet number, the MSD develops more rapidly at short times with increasing volume fraction. Finally, the volume fraction dependence disappears at long times. The diffusion coefficient was calculated from the long-time behavior of the MSD in the vorticity direction for different Peclet numbers and volume fractions. For the lowest Peclet number, the diffusion coefficient decreases with increasing volume fraction, while it is almost constant over the volume fraction change for the highest Peclet number. We estimated the structural relaxation time from the diffusion coefficient at the lowest Peclet number. The present non-Newtonian viscosity data agrees well with a simple scaling function for , similar to the non-Newtonian viscosity of a model supercooled liquid (35).

### References

- W. B. Russel, D. A. Saville, W. R. Schowalter, Colloidal dispersions (Cambridge University Press, Cambridge, UK, 1989).
- R. G. Larson, The Structure and Rheology of Complex Fluids (Oxford University Press, Oxford, 1999).
- A. Einstein, Ann. Phys. (Leipzig) 17, 459 (1905).
- I. M. Krieger and T. J. Doughty, Trans. Soc. Rheol. 3, 137 (1957).
- R. L. Hoffman, Trans. Soc. Rheol. 16, 155 (1972).
- B. J. Ackerson and P. N. Pusey, Phys. Rev. Lett. 61, 1033 (1988).
- M. D. Haw, W. C. K. Poon, and P. N. Pusey, Phys. Rev. E 57, 6859 (1998).
- H. Watanabe, M.-L. Yao, T. Shikata, H. Niwa, Y. Morishima, N.P. Balsara, and H. Wang, Rheol. Acta 37, 1 (1998).
- Brent J. Maranzano and N. J. Wagner, J. Chem. Phys. 117, 10291 (2002).
- J. F. Brady, G. Bossis, Ann. Rev. Fluid Mech. 20, 111 (1988).
- B. Dnweg and K. Kremer, J. Chem. Phys. 99, 6983 (1993).
- I.-C. Yeh and G. Hummer, J. Phys. Chem. B 108, 15873 (2004).
- D. M. Heyes, J. Phys: Condens. Matter 19, 376106 (2007).
- D. R. Foss and J. F. Brady, J. Fluid Mech. 407, 167 (2000).
- A. Sierou and J. F. Brady, J. Fluid Mech. 506, 285 (2004)
- D. I. Dratler and W. R. Schowalter, J. Fluid Mech. 325 1 (1996).
- D. I. Dratler, W. R. Schowalter, and R. L. Hoffman, J. Fluid Mech. 353 1 (1997).
- A. J. C. Ladd, Phys. Rev. Lett. 70, 1339 (1993).
- A. Malevanets and R. Kapral, J. Chem. Phys. 112, 7269 (2000).
- P. Ahlrichs and B. Dünweg, J. Chem. Phys. 111, 8225 (1999).
- H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 (2000).
- J. T. Padding and A. A. Louis, Phys. Rev. E 74, 031402 (2006).
- P. J. Atzberger, P. R. Kramer, and C. S. Peskin, J. Compt. Phys. 224, 1255 (2007).
- N. Sharma and N. A. Patankar, J. Compt. Phys. 201, 466 (2004).
- M. Fujita and Y. Yamaguchi, Phys. Rev. E 77, 026706 (2008).
- Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707 (2005).
- Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E 26, 361 (2008).
- T. Iwashita, Y. Nakayama and R. Yamamoto, J. Phys. Soc. Jpn. 77, 074007 (2008).
- T. Iwashita, Y. Nakayama and R. Yamamoto, Prog. Theor. Phys. Suppl. 178, 86 (2009).
- T. Iwashita and R. Yamamoto, Phys. Rev. E 79, 031401 (2009).
- A. Shakib-Manesh, P. Raiskinmaki, A. Koponen, and M. Kataja, and J. Timonen, J. Stat. Phys. 107, 67 (2002).
- J. Kromkamp, D. T. M. van den Ende, D. Kandhai. R. G. M. van der Sman, and R. M. Boom, J. Fluid Mech. 529, 253 (2005).
- P. M. Kulkarmi and J. F. Morris, Phys. Fluids 20, 040602 (2008).
- H. J. H. Clercx and P. P. J. M. Schram, Phys. Rev. A 46, (1992) 1942.
- R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515 (1998).