# A direct numerical simulation method for complex modulus of particle dispersions

## Abstract

We report an extension of the smoothed profile method (SPM) [Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E 26, 361 (2008)], a direct numerical simulation method for calculating the complex modulus of the dispersion of particles, in which we introduce a temporally oscillatory external force into the system. The validity of the method was examined by evaluating the storage and loss moduli of a system composed of identical spherical particles dispersed in an incompressible Newtonian host fluid at volume fractions of , , , and . The moduli were evaluated at several frequencies of shear flow; the shear flow used here has a zigzag profile, as is consistent with the usual periodic boundary conditions. The simulation results were compared with several experiments for colloidal dispersions of spherical particles.

###### pacs:

82.70.-yDisperse systems; complex fluids and 82.20.WtComputational modeling; simulation and 83.60.BcLinear viscoelasticity^{1}

## 1 Introduction

The viscoelastic properties of the dispersion of solid particles are of particular importance in several scientific, engineering, and industrial fields. These properties are strongly affected not only by direct interactions between particles, but also by thermal fluctuations of the system and the hydrodynamic interactions (HI) acting on dispersed particles, mediated by the surrounding fluid (1); (2). The viscoelasticity of materials are most commonly characterized by the complex modulus, which represents the mechanical response of materials to temporally oscillating small shear deformations of frequency . The complex modulus consists of elastic and viscous components (i.e., the storage modulus and the loss modulus ). In general, materials tend to lose elastic energy by dissipation when the frequency of external deformation is low, but can store elastic energy when the frequency is high.

Understanding viscoelastic behavior has been the subject of both fundamental and technological interest. Especially, the linear and non-linear responses of colloidal systems, including colloidal glasses and aggregated colloidal suspensions and gels, have attracted much attention in recent year (3); (4); (28). While extensive experimental studies have been carried out to determine the viscoelastic properties of particle dispersions over a wide range of volume fractions, from dilute () to very dense (), theoretical studies encounter a fundamental problem when the density of dispersion is high, because of many-body effects and the long-range nature of the HI among dispersed particles. Developing analytical theories becomes even more difficult when the host fluids are complex, such as electrolytes or viscoelastic media. Computer simulations can thus be very powerful tools in theoretical investigations of dense dispersions in general.

The widely used simulation technique is Stokesian dynamics (SD) method (5), which is based on the Stokes approximation (Re 0), involves near-field lubrication forces and far-field many-body HIs. The SD method has succeeded in simulating the motions of particle dispersions at steady shear. However the viscoelasticity of particle dispersions at oscillatory shear has not been examined. Furthermore, the SD method is valid only for simulating the particles in a Newtonian fluid and is not applicable to the motions of particles in complex fluids, such as charged particles in electrolytes or viscoelastic media.

Various numerical methods have recently been proposed for simulating dense dispersions, including some cases where the host fluids are complex. In order to accurately track the motions of host fluids as well as the motions of dispersed particles, a variety of numerical methods have been developed. One of these methods, which we call the direct numerical simulation (DNS) method for particle dispersions, involves solving the Navier-Stokes (NS) equation for host fluids in a manner consistent with boundary conditions defined according to the particle positions and motions. Aside from several successful implementations of DNS methods for particle dispersions (6); (7); (8); (9); (10); (11); (12); (13); (14), other methods have been developed that do not rely on the Navier-Stokes equation to resolve fluid motions. The coupling methods of particle dynamics with the lattice-Boltzmann (LB) method (15) or stochastic rotational dynamics (SRD) (16); (17) are the most popular alternatives to DNS, because LB and SRD are believed to be more computationally efficient.

Problems studied here have very small Reynold number and the full NS equations for the host fluid are not necessarily solved to simulate the motions of particles in a Newtonian fluid. It is sufficient to use the methods based on the Stokes approximation. However there are several advantages for solving the NS equations for the host fluid. Those techniques mentioned above enable us to simulate the short-time motions of particle dispersions, where the coupling between the particle motion and the fluid motion remains strongly, even if Reynolds number is very small. One of the typical behavior is the power-law decay in the velocity correlations of particles, which is known as ”long-time tail”. In addition, Those methods have succeeded in simulating several complex phenomena, such as the structure formation in colloidal gels where the pressure filed of the fluid plays an important role (18) and electrophoresis of charged colloidal dispersions (19); (20).

Although extensions of those methods have already been developed for steady shear flow with DNS (26), LB (21); (22); (23), and SRD(24), there exists no successful attempt for calculating the complex moduli of the dispersions with an imposed oscillatory shear flow. The purpose of the present paper is to report our successful reproductions of experimentally observed behaviors of the storage and loss moduli in temporally oscillating shear flow.

To check the validity of the method, we first apply the method to a Newtonian host fluid of viscosity , which should exhibit a purely viscous response, and . We then apply the method to dense dispersions composed of identical spherical particles in a Newtonian host fluid. Since several experimental measurements of the frequency dependence of and for dense colloidal dispersions have already been reported (28); (29); (30); (31), we compare our numerical results with those experiments, to examine the validity of our method for dense particle dispersions.

## 2 Simulation method

Let us consider a dispersion composed of identical spherical particles of radius in a Newtonian host fluid subjected to oscillatory shear. The host fluid is described by the velocity field and the pressure field . The th dispersed particle is described by , where is the position of the particle, is the translational velocity, and is the rotational velocity. The coupling scheme between the fluid motions and the particle motions is based on the smoothed profile method (SPM), which introduces a particle density field on the entire field, for the particle domains and for the fluid domains. These domains are separated by a thin interfacial domain of thickness . The SPM is an efficient method to resolve the HIs between the fluid and particle motions; its details are given in (9); (10); (12).

The time evolution of the host fluid is governed by the Navier-Stokes equation

(1) |

with the incompressibility condition , where is the density of the fluid, is the shear viscosity, the stress tensor , and is an external force field that is introduced to enforce an oscillatory shear flow on the entire system. The body force is introdued to ensures the rigidity of particles and the appropriate boundary condition at the fluid/particle interface (9); (10); (12).

The time evolution of the th particle with mass and moment of inertia is governed by Newton’s equations of motion:

(2) | |||||

(3) |

where and are the hydrodynamic forces and torques, respectively, exerted by the host fluid on the particle. is the repulsive force that prevents the particles from overlapping, and a truncated Lennard-Jones potential, for or , is adopted in this work. Here, and . and are the random force and torque, respectively, due to thermal fluctuations. These random fluctuations are assumed to be , and , where denotes time averaging, and and are numerical parameters to control the particle temperature . The procedure for determining the temperature is described in (25).

The apparent stress of the dispersion is written as

(4) |

where , is the density of the particles, and is the volume of the system. The derivation of Eq. (4) was reported in (26).

The apparent stress consists of two terms: the first term is a stress tensor including the external force and the second term is a stress tensor arising from the acceleration of the dispersion. In experimental viscoelastic measurements of dispersions, the acceleration term can be ignored, because the relaxation time scales related to the acceleration are considerably smaller than the experimental time scales. On the other hand, in the simulations based on the DNS approach, the acceleration term strongly affects the apparent stress of the dispersion.

There are two key points for calculating the apparent stress of the dispersion: i) how to calculate the acceleration term in , and ii) how to introduce the external force . The first point is straightforward. For a simulation time step , the acceleration term can be simply calculated as . Next the external force is introduced as a body force, to enforce the following oscillatory velocity field over the entire system,

(5) | |||||

(6) |

where denotes distance in the velocity-gradient direction, and is the oscillatory shear rate, with amplitude and frequency . Here, the flow is imposed in the direction and is the length of the system in the direction. This zigzag velocity profile was first used to simulate dispersions in a steady shear (25); (26). Note that is defined to be a body force that constrains the velocity field of the dispersions, and its explicit form is written as

(7) |

If the external force is introduced as a boundary force, then the development of the velocity from the boundary to bulk (i.e., the propagation modes) are observed for short time scales, such as the kinematic time scale . In the simulation, the propagation modes influence the overall viscoelastic properties of the dispersion. On the other hand, in experimental measurements the propagation modes are ignored. Using the external forces mentioned above, we can eliminate the propagation modes numerically.

To measure the storage modulus and loss modulus , we monitor the component of the apparent shear stress, , and shear rate as a function of time. In general, the component of the stress is written as

(8) |

where is the amplitude of the stress and is the phase difference between and . In our model, the shear rate is an externally controlled parameter. By using the obtained and , we can determine the following dynamic viscoelastic moduli:

(9) |

## 3 Simulation results

Three-dimensional simulations were performed at several frequencies ranging from 0.0005 to 0.2. The amplitude of strain, , was set to be 0.2. The spatial and temporal units are expressed in terms of the lattice spacing and , respectively. The system lengths are . The parameters of the simulation are , , , , , and . The dispersed particles are initially randomly distributed.

If we consider a particle of radius m in water at room temperature, the units of space and time correspond to be m and sec, respectively. In this case the simulated range of the frequency, , is between 7.1 and 2842 kHz.

### 3.1 Test of the simulation method

In order to test the validity of our method, we applied the method to a Newtonian host fluid that does not contain dispersed particles. Figure 1 (a) shows the time evolutions of the shear rate and shear stress at . It can be seen that both curves develop in time with the same phase, , . This behavior represents the typical features of Newtonian fluids. The loss modulus was then calculated using Eq. (9) for different frequencies. Figure 1 (b) displays the frequency dependence of the loss modulus for the host fluid. The modulus and frequency are non-dimensionalized by and , respectively, where is the self-diffusion coefficient of a Brownian particle at infinite dilution, . The loss modulus increases linearly with , and its slope is equal to the viscosity of the host fluid. The solid line indicates the loss modulus of the host fluid, . These results show that the viscoelastic properties of the host fluid are correctly reproduced by the simulation.

We next investigated a concentrated dispersion composed of spherical particles fluctuating in the host fluid. The volume fraction of the dispersed particles was . Figure 2 (a) shows the time evolutions of the shear rate and shear stress at and . The behavior of the shear stress of the dispersion differs from that of the host fluid shown in Fig. 1. We can see that there are phase differences between the shear rate and shear stress (i.e., ), and the amplitude of the shear stress becomes greater than that of the shear rate.

From the obtained , and , we can calculate the storage modulus and the loss modulus at different frequencies. Figure 2 (b) shows the frequency dependence of the storage modulus and the loss modulus for the dispersion at . The modulus and frequency are scaled by and , respectively. For low frequencies, increases linearly with and increases linearly with . As increases, grows monotonically until it reaches a plateau region, while develops up to a linear region, and its slope is smaller than the slope at low frequencies. The slope at high frequencies represents the high frequency viscosity , which is given by

(10) |

From Fig. 2 the high frequency viscosity was roughly estimated as . Next we focus on the high-frequency elastic shear modulus , which is defined as

(11) |

We compared our results with a theoretical expression for of hard spheres, which has been derived by Lionberger and Russel (27). To evaluate the theoretical value, we performed a numerical integration by using the obtained and Percus-Yevick distribution function for . The calculated high-frequency modulus, , is shown in Fig. 2. We can see the values of at high frequencies approach to the theoretical value.

Furthermore, the loss modulus is larger than the storage modulus, and the dispersion behaves like a viscous fluid. Here, we define the particle relaxation time to be , where is the diffusion coefficient of a spherical particle in the dispersion at thermal equilibrium. From the equilibrium calculations (), we can estimate the diffusion coefficient at , resulting in . The characteristic frequency is also indicated as an arrow in Fig. 2. At , the onset of elasticity is clearly observed. These behaviors accurately represent the typical viscoelastic features of concentrated dispersions (28).

Finally, we examined the dispersion behavior of particles at higher concentrations, and . The frequency dependence of and for different volume fractions is shown in Fig. 3; the data at and are also plotted for comparison. For , the structures of the particles are still randomly distributed under shear, wheres for the system forms crystallize phases completely. For and , each modulus shifts to higher values than those at . At low frequencies, the dependence of both moduli is remarkable: the storage modulus rises rapidly with increasing , causing the storage modulus to become larger than the loss modulus at low frequencies and . The dispersion then becomes elastic. Furthermore, the loss modulus at is a concave curve, with a minimum at low frequencies.

We found that the viscoelastic properties of the concentrated dispersions depend strongly on the volume fraction. In addition, at , the storage modulus increases very slowly with increasing , and at , a crossover from elastic- to viscous-dominant regions is observed.

### 3.2 Comparison of simulations and experiments

The simulation results obtained in this work were compared with experimental results for the dynamic viscoelasticity of colloidal dispersions, which were measured by several group (28); (29); (30); (31). These experiments use colloidal dispersions of spherical particles dispersed in solvent. All experimental data were also plotted as a function of dimensionless variables, , , and . In Fig.4, the simulation results for elastic and loss modulus at were plotted together with those measured experimentally by Shikata and Pearson (28). No fitting parameters are used. The viscoelastic responses of the simulations agree well with those of the experiment, although the volume fraction between them is not exactly the same. This clearly confirms our simulation can provide the typical viscoelastic behavior of colloidal dispersions in fluid states.

Comparison of simulation results and several experiments, including concentrated dispersions with higher , would give us a comprehensive information for the understanding of the dynamical behavior of colloidal dispersions. Figure 5 show the experimental data on the elastic and loss modulus of concentrated dispersions over a wide range of volume fractions from 0.37 to 0.56 (28); (29); (30); (31), and our simulation results were also plotted. These experimental data were measured over a range of from to , and our simulation results lie within the range. The experimental curve shifts at lower values of as the volume fraction increases and this reflects the characteristic time scales of system become longer, and the dynamics of dispersed particles slows down.

The experimental data for high volume fraction vary widely depending on each experiment, and it is considered that the viscoelastic response becomes sensitive to the details of direct interactions between particles. The repulsive potential used in this study can produce the features of dispersions with relatively small volume fraction (See Fig. 4). The simulation at has a crystalline structure over the whole frequency range, where a colloidal cystal is formed. Thus the simulation results at are completely different from the experiments, for which the structures are in fluid and glasslike states. Recently Crassous et al (29) have examined the effect of crystallization on the viscoelasticity of the dispersions at low frequency regime in the vicinity of the glass transition, and they found the elastic modulus is greater than the loss modulus at low frequencies even in fluidlike states. This behavior is qualitatively similar to the simulation results at , and it represents solidlike respnses.. In the present work the dispersions in glassy states were not studied.

## 4 Conclusion

We have developed a DNS method for simulating dynamics of solid particles dispersed in simple and complex fluids (9); (10). This method, called SPM, has also been successfully applied to simulate properties of dispersions under several non-equilibrium conditions, such as electrophoresis of charged spherical particles under external electric fields (19). The method was then modified to introduce thermal fluctuations into the dispersions (11), so that one can also simulate situations where the thermal fluctuations and hydrodynamic interactions acting among dispersed particles are both important. We have carried out systematic simulations for dispersions composed of spherical particles with and without steady shear flow, in order to analyze the diffusion process of dispersed particles in detail (25) and to investigate the nonlinear viscosity of the system (26).

In the present paper, we report an important extension of the SPM for analyzing the viscoelastic properties of particle dispersions immersed in host fluids, by introducing a temporally oscillatory external force into the system. To be consistent with the usual periodic boundary conditions, shear flow with a zigzag profile was employed. The validity of the method was examined by evaluating the storage and loss moduli of a system comprising identical spherical particles dispersed in an incompressible Newtonian host fluid at volume fractions of , , , and , for flow rates spanning a range of frequencies, . We confirmed that the method could successfully reproduce 1) purely viscous responses for and 2) typical viscoelastic responses for and , in excellent agreement with experimental data obtained for colloidal dispersions.

To our knowledge, the present study is the first successful attempt to calculate the complex modulus of particle dispersions using DNS-type methods. Further applications of our DNS method to more complex systems, such as dispersions of non-spherical particles, particles in polymer matrices, or dispersions of aggregating particles, are promising.

### Footnotes

- offprints:

### 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)
- H. M. Wyss, K. Miyazaki,, J. Mattsson, Z. Hu, D. R. Reichman, and D. A. Weitz, Phys. Rev. Lett. 98, (2007) 238303.
- V. Carrier and G. Petekidis. J. Rheol. 53, (2009) 245.
- J. F. Brady, G. Bossis, Ann. Rev. Fluid Mech. 20, (1988) 111.
- H. H. Hu, N. A. Patankar, and M. Y. Zhu, J. Comput. Phys. 192, (2001) 427.
- H. Tanaka and T. Araki, Phys. Rev. Lett. 85, (2000) 1338.
- T. Kajishima, S. Takiguchi, H. Hamasaki, and Y. Miyake, JSME Int. J., Ser. B 44, (2001) 526.
- Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, (2005) 036707.
- Y. Nakayama, K. Kim and R. Yamamoto, Eur. Phys. J. E 26, (2008) 361.
- T. Iwashita, Y. Nakayama and R. Yamamoto, J. Phys. Soc. Jpn. 77, (2008) 074007.
- X. Lio, M. R. Maxey, and G. E. Karniadakis, J. Comput. Phys. 228, (2009) 1750.
- P. J. Atzberger, P. R. Kramer, and C. S. Peskin, J. Comput. Phys. 224, (2007) 1255.
- M. Fujita and Y. Yamaguchi, Phys. Rev. E 77, (2008) 026706.
- A. J. C. Ladd, Phys. Rev. Lett. 70, (1993) 1339.
- A. Malevanets and R. Kapral, J. Chem. Phys. 112, (2000) 7260.
- J. T. Padding and A. A. Louis, Phys. Rev. E 74, (2006) 031402.
- R. Yamamoto, K. Kim, Y. Nakayama, K. Miyazaki, and D. R. Reichman, J. Phys. Soc. Jpn. 77, (2008) 084804.
- K. Kim, Y. Nakayama, and R. Yamamoto, Phys. Rev. Lett. 96, (2006) 208302.
- T. Araki and H. Tanaka, EPL, 82 (2008) 18004.
- A. Shakib-Manesh, P. Raiskinmaki, A. Koponen, M.Kataja, and J.Timonen, J. Stat. Phys. 107, (2002) 67.
- J. Kromkamp, D. T. M. van den Ende, D. Kandhai. R. G. M. van der Sman, and R. M. Boom, J. Fluid Mech. 529, (2005) 253.
- P. M. Kulkarmi and J. F. Morris, Phys. Fluids 20, (2008) 040602.
- M. Hecht, J. Harting, and H. J. Herrmann, Phys. Rev. E 74, (2006) 021403.
- T. Iwashita and R. Yamamoto, Phys. Rev. E 79, (2009) 031401.
- T. Iwashita and R. Yamamoto, Phys. Rev. E 80 (2009) 061402
- R. A. Lionberger and W. B. Russel, J. Rheol. 38, (1994) 1885.
- T. Shikata and D. S. Pearson, J. Rheol 38, (1994) 601.
- J. J. Crassous, M. Siebenbürger, M. Ballauff, M. Drechsler, D. Hajnal, O. Henrich, and M. Fuchs, J. Chem. Phys. 128, (2008) 204902.
- T. G. Mason and D. A. Weitz, Phys. Rev. Lett. 75, (1995) 2770.
- H. Watanabe, M. L. Yao, T. Shikata, H. Niwa, and Y. Morishima, Rheol. Acta 38, (1999) 2.