Stokesian Dynamics Simulations of Ferromagnetic Colloidal Dispersions in a Simple Shear Flow

JOURNAL OF COLLOID AND INTERFACE SCIENCE ARTICLE NO.

203, 233–248 (1998)

CS985498

Stokesian Dynamics Simulations of Ferromagnetic Colloidal Dispersions in a Simple Shear Flow Akira Satoh,* ,1 Roy W. Chantrell,† Geoff N. Coverdale,† and Shin-ichi Kamiyama‡ *Department of Mechanical Engineering, Faculty of Engineering, Chiba University, 1-33, Yayoi-cho, Inage-ku, Chiba 263, Japan; †School of Electronic Engineering and Computer Systems, University of Wales Bangor, Dean Street, Bangor, Gwynedd LL57 1UT, United Kingdom; and ‡Institute of Fluid Science, Tohoku University, 1-1, Katahira 2-chome, Aoba-ku, Sendai 980, Japan Received April 29, 1997; accepted February 26, 1998

We have investigated the behavior of clusters of ferromagnetic particles in a colloidal dispersion subjected to a simple shear flow. To do so, the Stokesian dynamics method has been used under the assumption that the effect of Brownian motion is negligible. For the case of no shear flow, the aggregate structures obtained by the Stokesian dynamics simulations agree well with Monte Carlo results qualitatively. We can, therefore, conclude that the Stokesian dynamics simulations can capture thick chainlike clusters without introducing a specific clustering algorithm, which is indispensable for Monte Carlo simulations. The behavior of the thick chainlike clusters in a simple shear flow is summarized as follows. The thick chainlike clusters decline in the shear flow direction as time advances. Since longer clusters experience larger shear forces, it is difficult for them to survive in such a situation. The thick chainlike clusters, therefore, dissociate into some short clusters. Such clusters are relatively stable in a shear flow, so that they do not decrease significantly any more. The viscosities have a strong relationship with the internal structures of the aggregates. The instantaneous viscosities, therefore, fluctuate significantly for the case of the thick chainlike clusters. q 1998 Academic Press Key Words: ferromagnetic colloidal dispersion; aggregation phenomena; Stokesian dynamics simulation; simple shear flow; thick chainlike cluster.

1. INTRODUCTION

In an ideal magnetic fluid or ferrofluid, ferromagnetic particles are uniformly dispersed in a base liquid. However, in an actual ferrofluid, the particles aggregate to form chainlike clusters along an applied magnetic field (1). Especially, water-based magnetic fluids show this tendency more clearly, and such aggregation phenomena can even be observed under an optical microscope (1), although particles dispersed in a base liquid are about 10 nm in diameter. Nakatani also has recently indicated the same result from a microscopic observation of aggregate structures (2). On the 1 To whom correspondence should be addressed. Fax: /81 043-2903039.

other hand, kerosine-based magnetic fluids have long been regarded as ideal magnetic fluids where there is no aggregate formation either in no external field or in an applied field. However, Buzmakov and Pshenichnikov (3) have recently pointed out from the measurements of the viscosity and the diffusion coefficient that even a kerosine-based magnetic fluid includes many microaggregates (or secondary particles) composed of several tens of primary particles in zero external magnetic field. Such secondary particles have a quasispherical shape and are very stable on temperature change. They interpreted that these aggregates are formed not due to magnetic forces, but molecular forces. Namely, they conclude that the secondary particles are formed at the stage of preparation due to the insufficiency of surfactant molecules on the surface of the primary particles. The magnetic moment of each particle forming a secondary particle behaves independently, and therefore the magnetic properties are described by the superparamagnetic theory. Also, they reported that chainlike clusters were not found for the case of no external field. It is noted that we arrived at the secondary particle model about ten years ago in order to explain the aggregate formation especially in a water-based magnetic fluid (4). Monte Carlo studies using the secondary particle concept have succeeded in predicting thick chainlike clusters (or microdroplets) in an applied magnetic field (5–7). This secondary particle model is conceptually the same of the microaggregates which Buzmakov and Pshenichnikov observed experimentally (3). Monte Carlo simulations (4–7) showed the process of the formation of thick chainlike clusters in an applied magnetic field in the following. Primary particles aggregate to form secondary particles in the preparation process due to an isotropic molecular attractive force distinct from the magnetic forces. Since primary particles are typically 10 nm in diameter, the magnetic moment of each particle is almost independent, which means that secondary particles have no significant magnetic moments in zero external field. An applied magnetic field induces the magnetic moment of each secondary particle, and then these secondary

233

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

0021-9797/98 $25.00 Copyright q 1998 by Academic Press All rights of reproduction in any form reserved.

coida

234

SATOH ET AL.

particles aggregate to form chainlike clusters. Furthermore, under a certain condition (6), such chainlike clusters aggregate to form further thicker chainlike clusters. This concept of aggregation has been supported by the results of the Monte Carlo simulations (5–7). A water-based magnetic fluid is less stable than a kerosine-based magnetic fluid, so that much larger secondary particles are expected to be formed in zero external field. If secondary particles are large enough to overcome thermal energy, it is possible for them to aggregate to form thick chainlike clusters due to magnetostatic forces in an applied field. We propose that these thick chainlike clusters correspond to those observed with a microscope. The relationship between the dimensions of secondary particles and the strength of magnetostatic forces has already been discussed in the previous paper (6). According to the study (6), a secondary particle can be regarded approximately as a unitary particle with a single moment, so that the results of Monte Carlo simulations (5, 7) can be applied to magnetorheological fluids which are composed of micron size ferromagnetic particles (8). It is clear from many experimental results (9–11) that aggregate formation significantly modifies the magnetic and rheological properties of the fluid. In order to apply ferrofluids to fluids engineering fields, it is very important to clarify microstructures in a quiescent flow field, microstructures in a flow field, and the influence of microstructures on rheological properties. In particular, since the application of an external magnetic field is central to their physical investigation and their practical application, it is important to understand the effects of a magnetic field on the microstructure. An applied field is known to enhance the degree of aggregation and to induce an anisotropic microstructure, which is clearly important in relation to the fluid rheology. As already mentioned, we have already conducted a series of studies concerning aggregate structures of ferrofluids in a quiescent flow by means of the cluster-moving Monte Carlo method ( 5 – 7 ) . It is found from these studies that there are three typical aggregate structures which depend on the strength of magnetic particle – particle interactions: ( I ) significant aggregation does not appear for the case of weak interactions between particles; ( II ) if particle – particle interactions are relatively strong, the particles aggregate to form thick chainlike clusters along an applied magnetic field; ( III ) finally, for the case of strong interactions, thin chainlike clusters are formed, but thick chainlike clusters are not. From a simulation point of view, the distinction between types II and III aggregation arises only when the clustermoving Monte Carlo method is adopted (5–7). This method carries out single-particle Monte Carlo moves followed by a cluster analysis and a series of Monte Carlo moves of the clusters themselves. This procedure is repeated until equilibrium is achieved. Thick chains do not occur when only single-particle Monte Carlo moves are used. The reason is

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

found in the form of the potential between two parallel chains of particles. At close range this is attractive for the case of two chains offset by the particle radius. At large separations the force is repulsive, the two chains being separated by an energy barrier. In type II clusters, the energy barrier is small enough to be overcome, leading to the existence of thick chains. In type III clusters, the interactions are strong enough to form extensive linear chains but the energy barrier is too large for thick chains to form. This interpretation highlights the importance of the cluster-moving Monte Carlo method and also serves in the present study by means of the Stokesian dynamics method. Although Monte Carlo methods are highly useful for the case of thermodynamic equilibrium, they are not appropriate for investigating dynamic properties such as the change in aggregate structures in a flow field. Thus there are two main reasons for the development of dynamic models. The first is to make a comparison with previous cluster-moving Monte Carlo studies. Second, it is of importance to study the effects of flow on the aggregate structures and rheological properties. There are three main methods for investigating dynamic properties of microstructures in a flow field. First, if hydrodynamic interactions and Brownian motion are negligible, the molecular dynamics method is applicable. Second, if only Brownian motion is negligible and hydrodynamic interactions are taken into account, the Stokesian dynamics method is used. Finally, if Brownian motion is not negligible, Brownian dynamics method has to be used. Generally speaking, Brownian dynamics simulations with hydrodynamic interactions are computationally the most expensive and require huge memory area on a computer. We here focus attention on the dynamic properties of a dispersion including chainlike clusters. Brownian motion is, therefore, not important in this situation, since magnetic particle–particle interactions are much more important than the effect of Brownian motion. Hence, We have used the Stokesian dynamics method in the present study to clarify the change of aggregate structures in a simple shear flow. The additivity of velocities is adopted for multibody hydrodynamic interactions between particles. 2. STOKESIAN DYNAMICS METHOD

2.1. Particle Model A particle is idealized as a spherical particle with central point dipole coated with a uniform surfactant layer (or a steric layer). A colloidal dispersion is assumed to be composed of such particles with radius a. As already pointed out, this particle model can also be regarded as a model of secondary particles composed of many primary particles. The magnetic interaction energy between particles i and j, u (ijm) , and the particle–field interaction energy, u (i H ) , are written, respectively, as

coida

235

STOKESIAN DYNAMICS SIMULATIONS

u ij( m) * Å u (ijm) /kT

SD ds rij

Ål

u

(H) i

*Åu

(H) i

F (ijV ) Å 0

3

{ni rnj 0 3(ni rtij )(nj rtij )},

/kT Å 0 jni rh,

u (ijV ) * Å u (ijV ) /kT

H

Å lV 2 0

SD

2rij /ds d ln td rij

m0m2 , 4pd 3s kT

jÅ

02

rij /ds 0 1 td

J

, [3]

m0mH , kT

lV Å

pd 2s ns , [4] 2

in which m0 is the permeability of free space and ns is the number of surfactant molecules per unit area on the particle surface. Now it is straightforward to derive the expressions for the forces and torques acting on a particle i using the above energy equations. We show here the final results: F (ijm) Å 0

1 d Å kT lV r tij ln d rij

[2]

in which d is the diameter of the particle including the steric layer (a is the radius) and td is the ratio of the thickness of the steric layer d to the radius of the solid part of the particle, equal to 2d /ds . The nondimensional parameters appearing in the above equations, l, j, and lV, are written

lÅ

SD

[1]

in which ni is the unit vector denoting the direction of the magnetic moment mi (m Å Émi É), tij is the unit vector given by rij /rij , rij Å ri 0 rj , rij Å ÉrijÉ, h Å H/H, H is the applied magnetic field (H Å ÉHÉ), ds is the diameter of the particle excluding the steric layer (the radius is expressed by as ), k is Boltzmann’s constant, and T is the absolute temperature of fluid. The superscript ‘‘*’’ denotes dimensionless quantities. Repulsive forces arise due to the overlapping of the steric layers. The expression for the interaction energy for this case, u (ijV ) , has already been developed by Rosensweig et al. (12) as

Ì ( m) Ì ( m) 3m0m2 1 u ij Å 0 u ij Å 0 r 4 Ìri Ìrij 4pd (rij /d) 4

1 [ 0 (ni rnj )tij / 5(ni rtij )(nj rtij )tij 0 {(nj rtij )ni / (ni rtij )nj }],

[5]

Ì (V) Ì (V) u ij Å 0 u ij Ìri Ìrij

, (ds ° rij ° d),

[8]

in which ni and tij have already been defined after Eq. [2], F (ijm) and T (ijm) are the magnetic force and torque acting on particle i from neighboring particles, T (i H ) is the torque acting on particle i due to the deviation of the magnetic moment from an applied magnetic field in direction, and F (ijV ) is the force due to the overlapping of the steric layers. Hij is the magnetic field at the position of the particle i made by the dipole moment of the particle j. 2.2. Hydrodynamic Interactions between Two Spheres Since the relative velocity between a particle and the fluid is generally much less than the fluid velocity (namely, the Reynolds number, Re, satisfies Re ! 1), the flow is governed by the Stokes equations (13). This means that one can obtain analytical solutions of the flow field, since the Stokes equations are linear partial differential equations. The solutions for a two-solid-particle system in a simple shear flow have already been obtained by some authors (14–16). If the velocities and angular velocities of particles 1 and 2 are denoted by v1 and v2 , and v1 and v2 , these quantities are expressed in terms of the forces, F1 and F2 , and torques, T 1 and T 2 , of the particles acting on the ambient fluid (16), n1 0 U(r1 ) n2 0 U(r2 ) v1 0 V v2 0 V

Å

a11 a12 bH 11 bH 12 gI 1 a21 a22 bH 21 bH 22 gI 2 b11 b12 c11 c12 hH 1 b21 b22 c21 c22 hH 2

F1 / h F2 / h T1 / h T2 / h E

,

[9]

where r1 and r2 are the position vectors of the particles, h is the viscosity of the base liquid, and U and V are the velocity and angular velocity of the fluid in the absence of the particles. We here consider a simple shear flow, so the velocity vector U can be written as U(r) Å V 1 r / Err, in which E is the rate-of-strain tensor. Also, a, b, c, g, and h are called mobility tensors. For example, the mobility tensor a12 is written in the form

2

T (ijm) Å m0mi 1 Hij Å 0

m0m 1 r 3 4pd (rij /d) 3

1 {ni 1 nj 0 3(nj rtij )ni 1 tij },

T

AID

(H) i

Å m0mHni 1 h,

JCIS 5498

/

a12 Å x a12 t21t21 / y a12 (I 0 t21t21 ), [6] [7]

6g44$$$461

06-10-98 08:54:41

[10]

in which x a12 and y a12 are called mobility functions and are dependent only on the distance between particles. I is the

coida

236

SATOH ET AL.

unit tensor. The other mobility tensors can be written in similar forms (16). In the present study, the values of the mobility functions have been evaluated according to the numerical method shown by Kim and Mifflin (16). The tabulated data are referred to and interpolated to evaluate the values necessary in the run of the simulations. 2.3. Additivity of Velocities for Multibody Hydrodynamic Interactions Multibody hydrodynamic interactions can be obtained by extending the theory of a two-particle system in order to deal with a dense colloidal dispersion. The additivity of forces and the additivity of velocities are mainly used for taking into account the multibody hydrodynamic interactions. Since the former approximation can more accurately reproduce the lubrication effect which arises when particles nearly touch, Brady, Bossis, and their colleagues have therefore been vigorously developing the Stokesian dynamics method based on this approximation (17–20). However, since the inverse of the resistance matrix has to be calculated according to their method, the simulations are therefore limited to a small system; for example, the number of particles, N, is 25 and 27 in the papers (21, 22), respectively (several cases are for N Å 123 in (22)). Although the additivity of velocities is an inferior approach as concerns the accuracy of the lubrication effect, one can expand the simulations to a larger system (for example, N Å 1000). Since we are interested in aggregation phenomena in a flow field, a larger system is desirable. Hence the additivity of velocities is adopted in the present study. In this approximation, the velocity vi and angular velocity vi of particle i are written (16,23) ni 0 U(ri ) Å

1 h

F

1 Fi / 6pa

N

∑ j Å1 ( xi )

N

∑

/

G

N

D

1 I r Fi 6p a

j Å1 ( xi )

F

∑ j Å1 ( xi ) N

/

N

N

j Å1 ( xi )

j Å1 ( xi )

∑ j Å1 ( xi ) N

/

S

cii 0

G

JCIS 5498

/

D

1 Ti 8pa 3

Each quantity can be written in dimensionless forms in the usual way for a simple shear flow with the following representative values: the length is a (the particle radius), the time is 1/ gh , the velocity is gh a, the force is 6phagh a, the angular velocity is gh , the torque is 8pha 3gh , and so forth. Then the dimensionless velocity and angular velocity are written N

n* i Å y*dx / F* i /

N

06-10-98 08:54:41

∑ (a*ii 0 I)r F*i j Å1 ( xi )

[12]

N

j Å1 ( xi )

6g44$$$461

[13]

j Å1 ( j xi )

Once the velocities and angular velocities at time t are known, one can evaluate these values at the next time step; for example, ri (t / Dt) Å ri (t) / Dtvi (t). Our particle model has a steric layer, not a solid one. The solutions of the mobility tensors for this kind of particle have not been obtained yet, so the expressions for the solid particle of diameter d ( Å 2a) are used in the present study. However, they are used together with the procedure that the hydrodynamic interactions between particles are regarded as zero if the steric layers overlap. This procedure does not give a significant influence on aggregate structures, since the repulsive interactions due to the overlapping of the steric layers have more dominant factor than the lubrication effect in the overlapping situation. Finally, we have to point out that the approximation of the additivity of velocities does not necessarily satisfy the positive definiteness of the mobility matrix. This drawback may cause a serious error in determining the particle positions and velocities for the case of a usual colloidal system; this may be why there are few Stokesian dynamics simulations based on the additivity of velocities for nonmagnetic colloidal dispersions. However, for the case of a strongly interacting system such as ferrofluids and electrorheological fluids, the positive definiteness of the mobility matrix is essentially satisfied if the cutoff radius for hydrodynamic interactions between particles is introduced. We will mention this matter in more detail in the latter part of Section 3.

/

∑ a*ij r F*j / 2 j Å1 ( xi )

where the fluid flows in the x direction with shear rate gh , U(ri ) Å gh yi dx , V Å 0 ( gh /2) dz , and ( dx , dy , dz ) are unit vectors. The force and torque on the ambient liquid from

AID

N

Ti Å ∑ T (ijm) / T (i H ) .

j Å1 ( j xi )

1 I r Ti 8pa 3

∑ cijrTj / ∑ hH i :E,

j Å1 ( xi )

[11]

j Å1 ( xi )

∑ biirFi / ∑ bijrFj / N

N

Fi Å ∑ (F (ijm) / F (ijV ) ),

2.4. Nondimensional Numbers Characterizing Phenomena bH iirTi

∑ bH ijrTj / ∑ gI i :E,

/

1 h

aii 0

N

aijrFj /

j Å1 ( xi )

vi 0 V Å

S

the particle i, namely, the force Fi and the torque Ti acting on the particle i from the ambient particles, are as follows:

N

/

N

∑ bH *ii r T*i

j Å1 ( xi ) N

∑ bH *ij r T*j / 2 ∑ gI *i :E*,

j Å1 ( xi )

coida

D

S

j Å1 ( xi )

[14]

237

STOKESIAN DYNAMICS SIMULATIONS

1 3 v* i Å 0 dz / 2 2

S

N

N

j Å1 ( xi )

j Å1 ( xi )

∑ b*ii r F*i / ∑ b*ij r F*j

D

N

∑ (c*ii 0 I)r T*i

/ T* i /

j Å1 ( xi ) N

N

j Å1 ( xi )

j Å1 ( xi )

/

∑ c*ij r T*j / ∑ hH *i :E*,

[15]

where the dimensionless mobility tensors are a* ij Å 6paaij ,

2 b* ij Å 4pa bij ,

g* ij Å gij /2a,

h* ij Å hij .

3 c* ij Å 8pa cij ,

[16]

The dimensionless force F* i and torque T* i are written N

( m) (V) F* i Å ∑ (F ij * / F ij * ),

N

( m) (H) T* *, i Å ∑ T ij * / T i

j Å1 ( j xi )

j Å1 ( j xi )

[17]

FIG. 1. Potential curves for magnetic fluids.

where F ij( m) * Å 0 Rm

8 4 r* ij

values of the magnetic moments, the applied magnetic field, and the properties of the steric layers. The expressions for RH /Rm and RV /Rm are written

[ 0 (ni rnj )tij / 5(ni rtij )

1 (nj rtij )tij 0 {(nj rtij )ni / (ni rtij )nj }],

F T

(V) ij

( m) ij

* Å RV tij ln * Å 0 Rm

S D

2 3 r* ij

2

,

r* ij

(2/(1 / td ) ° r* ij ° 2),

{ni 1 nj 0 3(nj rtij )ni 1 tij },

T (i H ) * Å RHni 1 h.

[18] [19] [20] [21]

The nondimensional numbers, Rm , RH , and RV, appear due to the forces and torques being normalized by the viscous force of the shear flow, and are written Rm Å

m0m2 , 64p 2ha 6gh

RH Å

m0mH , 8pha 3gh

RV Å

kTlV . 6pha 2gh d

[22]

Rm means the ratio of the representative magnetic particle– particle force to the representative hydrodynamic shear force; RH means the ratio of the representative torque due to an applied field to the hydrodynamic torque, and RV means the ratio of the representative repulsion due to steric layers to the hydrodynamic shear force. Since the nondimensional shear rate gh * is always unity, we cannot change the flow rate using this quantity. However, by changing the value of Rm with constant values of RH /Rm and RV /Rm , one can change the shear rate with constant

AID

JCIS 5498

/

6g44$$$461

RH /Rm Å

06-10-98 08:54:41

(1 / td ) 3 j r , l 4

RV /Rm Å

(1 / td ) 4 lV r , l 3td

[23]

where l, j, and lV are shown in Eq. [4]. It is clear from Eq. [23] that RH /Rm and RV /Rm are dependent only on particle properties, not on flow properties. Also, it is noted that the change of Rm Å 0 r ` means that the shear rate changes from ` to zero. 3. PARAMETERS FOR SIMULATIONS

The following values were used in conducting Stokesian dynamics simulations. Three values of RV /Rm are chosen (95.2, 52.89, 23.8), which are obtained for l Å (5, 9, 20) to compare the present results with the previous ones obtained by the Monte Carlo simulations (7). We here set the steric properties as lV Å 150 and td Å 0.3. The potential curves for this case are shown in Fig. 1; s is the surface-tosurface distance. The minimum values of the potential curves are about (5, 9, 20)kT for l Å (5, 9, 20), respectively. It is therefore expected from our previous Monte Carlo results (7) that thick chainlike clusters are formed along an applied magnetic field for RV /Rm Å 52.89. We consider a strong magnetic field as RH /Rm Å 3.3, which corresponds to j Å 30 for l Å 5. This means that particle–field interactions are much larger than particle–particle ones. The number of

coida

238

SATOH ET AL.

particles, N, is 1000, and the number density, n* ( Å na 3 ), is 0.01. For this case, the length of the simulation box (cube), L* ( Å L/a), is 46.42. The cutoff radius for particle– particle interactions, r * cof f ( Å rcof f /a), is 16. The values of the time interval Dt* have to be selected carefully. If the time interval is not much shorter than the characteristic times for shear flow and for the particle motion due to magnetic interactions, the particles may overlap to an unreasonably large extent thereby causing large repulsive forces and divergence of the system. Since the force in the nondimensional system is characterized by Rm , the characteristic time Dt* for the case of the nondimensional force Rm and the length DL* is expressed as Dt* Å DL*/Rm . If we take DL* as DL* Å 0.1td , the characteristic time can be obtained as Dt* Å 0.1td /Rm . Also, the time interval Dt* has to be much shorter than the characteristic time for a shear flow, so we adopted the time interval as Dt* Å min(0.001, Dt*). A shear flow is generated by Lees–Edwards periodic boundary condition (24), namely, the upper and lower periodic boxes are made to slide in the shear direction. Next, we mention about the necessity of the cutoff radius ( hydro ) . The for the hydrodynamic interactions, denoted by r cof f additivity of velocities is just an approximation, so that some restrictions may be introduced together with this approximation. In the level of the additivity of velocities, the hydrodynamic interactions between particles i and j are developed neglecting the influence of the third particle. If the particles i and j are separated enough for a third particle to intervene between them, the influence of particle j on particle i is screened by the existence of the third particle. Taking into account this kind of configuration, therefore, may cause unrealistic aggregate structures. This is why we introduced the ( hydro ) cutoff radius for the hydrodynamic interactions, r cof . Acf cording to this approximation, the hydrodynamic interactions are ignored if the particle distance rij becomes larger ( hydro ) ( hydro ) ( hydro ) than r cof . We here took r cof as r cof Å 1.9d ( Å 2 f f f 1 1.9a). Although the introduction of the cutoff radius eventually means the negligence of the far-field effect, we presume that this influence on the aggregation formation is not so significant, because we here consider a strongly interacting system due to magnetostatic forces. In this case, the aggregation formation is mainly governed by magnetostatic interactions and the lubrication effect. We mentioned about the possibility of the loss of the positive definiteness of the mobility matrix for the case of the additivity of velocities in the last part of Section 2.3. The loss of the positive definiteness strongly depends on the particle arrangement in the system; a group of second or third nearest particles gathering around a particle of interest causes a loss of positive definiteness. Thus the introduction of the cutoff radius for the hydrodynamic interactions significantly improves this situation, and the positive definite-

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

ness of the mobility matrix, therefore, is substantially satisfied. Finally, we mention about the computation time for some typical cases. Many simulations were conducted on a standard workstation (RS6000, produced by IBM). A total of about 20,000 time steps were needed to obtain the aggregate structure shown in Fig. 2a, and, for this case, the CPU time was about 52 h. Also, the CPU time was about 240 h to obtain the instantaneous viscosity curve for Rm Å 100, which is shown in Fig. 7a, and the number of total time steps was about 100,000. 4. RESULTS

4.1. Aggregation Structures in a Quiescent Flow Figure 2 shows the results for Rm Å 1 1 10 7 and Dt* Å 3 1 10 09 , which corresponds to zero shear rate. Figure 2a and 2b are for RV /Rm Å 52.89 and 23.8, respectively. The figure on the left-hand side is viewed from an oblique angle, and the figure on the right-hand side from the magnetic field direction. The magnetic field is applied in the y direction (in the upper direction on paper). The value of 52.89 is chosen to compare the present results with the previous ones obtained by the Monte Carlo simulations (7), in which thick chainlike clusters were observed. If the aggregate structures shown in Fig. 2a are viewed from various directions, thick chainlike formation can be recognized. Hence, it is seen that the Stokesian dynamics simulations can also capture thick chainlike clusters. They are, however, somewhat less pronounced than our previous Monte Carlo results (7), which we ascribe to the fact that this previous study does not take into account the repulsive forces due to the steric layers. For the case of very strong magnetic particle–particle interactions, relatively thin chainlike clusters are formed along the field direction, as shown in Fig. 2b. These clusters are not so thick compared with those for the previous case. We may conclude from these results that Stokesian dynamics results agree well with Monte Carlo results qualitatively. It is interesting to note, as described earlier, that our previous study found it necessary to invoke cooperative motion of several particles using a cluster-moving Monte Carlo approach. This initial study of quiescent flow properties confirms that the Stokesian dynamics approach successfully and naturally takes account of the cooperative motion in predicting thick chainlike structures. This is of course an important feature in studying the flow properties of the system, which forms the major aim of the current work. 4.2. Aggregation Structures in a Shear Flow Figure 3 and 4 show the change of aggregate structures subjected to a simple shear flow, from the initial configuration, given in Fig. 2a. Figure 3 is obtained for Rm Å 100, and Fig. 4 is for Rm Å 1. The values of D X/L mean the

coida

STOKESIAN DYNAMICS SIMULATIONS

239

FIG. 2. Aggregation structures in a quiescent flow: (a) RV /Rm Å 52.89; (b) RV /Rm Å 23.8.

travel length of the upper and lower periodic boxes in the shear direction. It is seen from Fig. 3 that the thick chainlike clusters decline in the shear direction as time advances. At

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

longer times, these thick chainlike clusters dissociate into short clusters, which are relatively stable. Since the applied magnetic field is very strong, the directions of the magnetic

coida

240

SATOH ET AL.

FIG. 3. Aggregation structures in a shear flow for Rm Å 100: (a) D X/L Å 0.2; (b) D X/L Å 0.6; (c) D X/L Å 2.0.

AID

JCIS 5498

/

6g44$$5498

06-10-98 08:54:41

coida

STOKESIAN DYNAMICS SIMULATIONS

241

FIG. 3—Continued

moments almost coincide with the field direction. This behavior of the thick chainlike clusters is very physically reasonable. Namely, longer clusters experience larger shear forces, so that it is difficult for long clusters to survive in such a situation. Once long clusters dissociate into short clusters of appropriate length, such clusters are relatively stable. Hence, after that, the aggregate structures do not change significantly. Figure 4 shows the results for the case of a strong shear flow. In this case, the characteristic time for magnetic interactions is much longer than that for the shear flow, so that shear forces dominate and all particles move in the shear direction. This means that the clusters disappear as time increases. This fact will become much more clear later when we consider the results of correlation functions. It is noted that the magnetic moments remain pointing to the field direction even at long times. For the case of extremely strong shear flow field (not shown here), the spatial configuration is almost the same as the previous one. But, for this case, the shear flow is very strong, so that the particles themselves rotate with the macroscopic angular velocity. 4.3. Pair Correlation Functions We have investigated the above-mentioned aggregate structures in terms of pair correlation functions. Figures 5 and 6 show the change in the pair correlation functions with time for

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

RV /Rm Å 52.89. Figure 5 is obtained for Rm Å 100, and Fig. 6 is for Rm Å 1. In these figures, the u axis is taken from the field direction to the shear flow direction. The sharp peaks at r/a á 2, 4, 6, . . . , and u á 0 mean that significant chainlike clusters are formed along the field direction. It is seen from Fig. 5 that these correlation functions do not change significantly after D X/L Å 0.6. This suggests that the dissociation process has almost finished at this time. So, after that, the aggregate structures do not essentially change. This is in good agreement with the snapshot results shown in Fig. 3. Namely, after the long chainlike clusters dissociate into some short clusters, the aggregate structures do not change significantly. So we may say that the cores of the thick chainlike clusters do survive under the influence of the shear flow. Next, we consider the results shown in Fig. 6. Each peak moves in the u direction, decreases in height, and finally disappears. This means that the clusters decline in the shear flow direction as time advances, and the aggregate structures disappear finally. These properties of the pair correlation functions agree well with those shown in Fig. 4. Figure 6 is extremely informative as to the detail of the breakup of clusters in the shear flow. At D X/L Å 0.6, a second nearest-neighbor peak is still apparent, but this has disappeared by D X/L Å 1, leaving clusters in the form of dimers. Rather than disappearing immediately the relatively stable dimer rotates into the flow direction and presumably

coida

242

SATOH ET AL.

FIG. 4. Aggregation structures in a shear flow for Rm Å 1: (a) D X/L Å 0.6; (b) D X/L Å 1.0; (c) D X/L Å 2.0.

AID

JCIS 5498

/

6g44$$5498

06-10-98 08:54:41

coida

STOKESIAN DYNAMICS SIMULATIONS

243

FIG. 4—Continued

break up because this is an energetically unfavorable orientation from the point of view of the magnetostatic energy. We have the similar results for the case of the thin chainlike formation, but the conclusions are similar to those for the previous case. So we omit these results in this paper. 4.4. Viscosities Finally, we discuss the results of viscosity calculations. Figure 7 shows the evolution of the instantaneous viscosities with time. Figure 7a is for RV /Rm Å 52.89, and Fig. 7b is for RV /Rm Å 23.8. The instantaneous shear viscosity hyx is generally defined using the off-diagonal element tyx of the apparent stress tensor t and the shear rate gh as hyx Å tyx / gh . We here concentrate our attention on the instantaneous viscosity due to the forces and the torques between m m . The expression of h yx can be written as (23) particles, h yx N

m h yx *Å

m h yx 6p Å0 ∑ h V *iÅ1

N

∑ y*ij F*ijx / j Å1 ( j úi )

N

4p ∑ T*iz , [24] V *iÅ1

where V * is the system volume. The first term on the righthand side in Eq. [24] is due to the particle–particle forces, and the second term is due to the torques. It is noted that the contribution from the stresslet is not considered in the present study. Each figure has two curves: one is for a very small shear

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

rate, and the other is for a small shear rate. It is seen from these curves that, at the initial stage, the viscosities increase, and then converge to a certain value. The increase in the viscosities at the initial stage results from the deformation of the internal structures of the chainlike clusters. After the long chainlike clusters dissociate into short clusters, the resistance of such clusters becomes almost constant. Therefore the viscosity values themselves also converge to a specific (constant) values. We can also see from these figures that there exist fluctuations of the instantaneous viscosity which is very large for the thick chainlike cluster formation in comparison with those for thin chainlike formation. This is very physically reasonable, because thick chainlike clusters have a tighter structure, so they have stronger resistance to the deformation due to the shear forces. Also, it is seen that the fluctuation in the instantaneous viscosities is larger for larger values of Rm . This means that, in order to obtain viscosity values with sufficient accuracy, an ensemble average procedure with various initial configurations is indispensable. Figure 8 shows the evolution of the ensemble averaged m viscosities, » h yx … , with time. Such ensemble averaged viscosities were obtained, first from the ensemble average of the results obtained by the simulations with six different initial configurations and then the subaverage of those data every D X Å 0.1L. Figure 8a is for RV /Rm Å 52.89, and Fig. 8b is for RV /Rm Å 23.8. Although the data for Rm Å 500 do

coida

244

SATOH ET AL.

FIG. 5. Pair correlation functions for Rm Å 100; (a) D X/L Å 0.2; (b) D X/L Å 0.6; (c) D X/L Å 2.0.

not have sufficient accuracy yet, we can see the following m qualitative properties of the viscosities » h yx … . The property that the viscosities increase at the initial stage of time and then converge to a certain value, becomes more clear in Fig.8 than in Fig.7. For the both cases, Fig. 8a and b, the values at D X/L Å 1 depend on the values of Rm : specifically, the values of the viscosities increase as the shear rates decrease. Hence, the present colloidal dispersion has non-Newtonian properties, which agrees with the theoretical and experimental results (9-11). For the case of the thick chainlike m cluster formation, Fig. 8a, a large value of » h yx … persists, up to large strains of D X/L á 0.8, for Rm Å 500 compared with the case of the thin chainlike cluster formation, Fig. 8b. As pointed out previously, this is due to the property of

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

the thick chainlike formation that the thick chainlike clusters have stronger resistance to the deformation owing to the shear forces. Figure 9 shows the ensemble averaged viscosities in equilibrium for RV /Rm Å 52.89 as a function of the shear rate. Each data is shown with an error bar, which was estimated according to the method described in the textbook (25). It is clearly shown in Fig. 9 that the colloidal dispersion indicates nonNewtonian properties, as already pointed out, and the viscosity increases with the values of Rm , or the inverse of the shear rate. The viscosity begins to increase at Rm á 1, where the magnetostatic forces start to become a main factor for determining the cluster formation compared with the shear forces. We took a relatively dilute dispersion (n* Å 0.01) in the present

coida

STOKESIAN DYNAMICS SIMULATIONS

245

FIG. 6. Pair correlation functions for Rm Å 1: (a) D X/L Å 0.6; (b) D X/L Å 1.0; (c) D X/L Å 2.0.

study, so that a significant increase in the viscosity is not obtained, only about three times the viscosity of the base liquid. If we consider a dense colloidal dispersion, a large viscosity increment is expected to be observed. Now we compare the present result with experimental data. The viscosity of ferrofluids is strongly dependent on the volumetric fraction, the strength of the external magnetic field, the temperature, etc., so that experimental data for a similar condition has to be chosen to compare. Figure 10 shows the experimental data of the reduced rotational viscosity as a function of the strength of the magnetic field, which was obtained by Odenbach and Sto¨rk (26) using their origi-

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

nal magnetic fluid rheometer for the volumetric fraction of 0.07. The volumetric fraction of our colloidal dispersion is 0.042, and therefore not significantly different from that of Odenbach and Sto¨rk’s dispersion. Since we concentrate our attention on a very strong magnetic field case, we compare the viscosity data for H á 35 [mT] in Fig. 10 with our result in Fig.9. The reduced viscosity in Fig. 10 is read to be about (2.7, 2.3, 1.2) for gh Å (1.05, 2.62, 5.23) [s 01 ], respectively. We are using nondimensional quantities, so that a criterion point of Rm has to be chosen before comparing the viscosity results. We here regard the case of Rm Å 40 as a criterion for comparison, for which the viscosity data

coida

246

SATOH ET AL.

FIG. 7. Instantaneous viscosities: (a) RV /Rm Å 52.89; (b) RV /Rm Å 23.8.

of the simulation is almost equal to the experimental data, 1.2. Figure 9 shows the viscosity data, about (2.8, 1.9, 1.2), for Rm Å (200, 100, 40), which correspond to gh Å (1.05, 2.62, 5.23) [s 01 ] in the experiment. It is seen that the simulation data agree well with the experimental ones qualitatively and also quantatively. Our results were obtained under the condition that the magnetic field is extremely strong, so that

the experimental data are expected to be larger for such a strong magnetic field as giving a saturation magnetization. This means that the experimental data are supposed to show larger values than the simulation values, if we compare these results for the same case of an extremely strong field. This is quite understandable because the experimental dispersion is more dense than the simulation one. The simulations for electrorheological fluids are similar to those for ferrofluids, although the dispersed particles are significantly different in dimensions. Hence, if the nondimensional quantities are used in such simulations, we may compare the present results with those for electrorheological fluids. Unfortunately, there are few studies of computer simulations, in which hydrodynamic interactions are taken into consideration and a large system such as N Å 500–1000 is adopted. The results obtained by Melrose (27) are for N Å 108 and 256 and for the volumetric fraction of 0.31, but the hydrodynamic interaction among particles are not considered. Since their dispersion of interest, therefore, is much more dense than ours, we here attempt to compare the results qualitatively. First, we have to note that the nondimensional number used in his paper, Mn , corresponds to ( 18 Rm ) in our paper. It is seen from his paper that the reduced viscosity begins to increase monotonically at Rm á 1 from zero and attains 10 at Rm á 50, which is very similar to our result shown in Fig. 9. The significant increase in the viscosity seems to be physically reasonable because he concentrated his attention on a much more dense system than ours. We can, therefore, conclude that our result shown in Fig. 9 agrees well with Melrose’s curve qualitatively. Finally, we briefly discuss the dependence of the particle number of the system, N, the cutoff radius for magnetostatic interactions, rcof f , and the cutoff distance for hydrodynamic ( hydro ) interactions, r cof , on the ensemble averaged viscosity. f The viscosity data for various cases of such parameters are given in Table 1 for Rm Å 100 and RV /Rm Å 52.89. It is seen from Table 1 that the particle number and the cutoff

FIG. 8. Ensemble averaged viscosities: (a) RV /Rm Å 52.89; (b) RV /Rm Å 23.8.

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

coida

247

STOKESIAN DYNAMICS SIMULATIONS

TABLE 1 Dependence of the Number of Particles, N, the Cutoff Radius for Magnetic Forces, rcoff, and the Cutoff Distance for Hydrody(hydro) namic Interactions, r coff , on the Ensemble Averaged Viscosity, m »hyx…/h, for Rm Å 100 and RV/Rm Å 52.89

FIG. 9. Ensemble averaged viscosities in equilibrium for RV /Rm Å 52.89 as a function of shear rate.

radius for magnetic forces do not give significant influence on the viscosity data. It has already been pointed out in the previous paper (6) that the cutoff radius of rcof f Å 16 is sufficiently large for obtaining physically reasonable cluster formation. This is supported by the data shown in Table 1 again. Also, we see that the viscosity data are significantly dependent on the cutoff radius for hydrodynamic interactions. As pointed out in Section 3, it is not appropriate to use a value of the cutoff radius larger than 2d for the present system. Table 1 shows that a small cutoff radius leads to a large value of the viscosity. We plan to do further studies about the dependence of the cutoff radius on the cluster formation and the viscosity. 5. CONCLUSIONS

We have investigated the behavior of clusters of ferromagnetic particles in a colloidal dispersion subjected to a simple

N

r*coff

(hydro)* rcoff

»hmyx…/h

1000 729 512

16 16 16

3.8 3.8 3.8

1.86 { 0.07 1.53 { 0.16 1.88 { 0.09

1000 1000

12 10

3.8 3.8

1.98 { 0.09 2.12 { 0.08

1000 1000 1000

16 16 16

3.4 3.0 2.6

2.72 { 0.08 4.56 { 0.17 5.00 { 0.14

shear flow. To do so, the Stokesian dynamics method has been used under the assumption that the effect of Brownian motion is negligible. For the case of no shear flow, the aggregate structures obtained by the Stokesian dynamics simulations agree well with Monte Carlo results qualitatively. We can, therefore, conclude that the Stokesian dynamics simulations can capture thick chainlike clusters without introducing a specific clustering algorithm, which is indispensable for Monte Carlo simulations. The behavior of the thick chainlike clusters in a simple shear flow is summarized as follows. The thick chainlike clusters decline in the shear flow direction as time advances. Since longer clusters experience larger shear forces, it is difficult for long clusters to survive in such a situation. The thick chainlike clusters, therefore, dissociate into some short clusters. Such clusters are relatively stable in a shear flow, so that they do not decrease significantly any more. The viscosities have a strong relationship with the internal structures of the aggregates. The instantaneous viscosities, therefore, fluctuate significantly for the case of the thick chainlike clusters. ACKNOWLEDGMENTS The present work was done while the first author (A.S.) was staying at the University of Wales Bangor and Keele University, and A.S. gratefully acknowledges the financial support of the British Council for the one-year stay in UK. We also acknowledge the financial support of EPSRC.

REFERENCES

FIG. 10. Reduced rotational viscosity obtained experimentally by Odenbach and Sto¨rk (26).

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

1. Hayes, C. F., J. Colloid Interface Sci. 52, 239 (1975). 2. Nakatani, I., in ‘‘Proceedings of the International Symposium on Hydrodynamics of Magnetic Fluids and its Applications’’ (S. Kamiyama and R. Massart, Eds.), p. 9. Tohoku University, Japan, 1997. 3. Buzmakov, V. M., and Pshenichnikov, A. F., J. Colloid Interface Sci. 182, 63 (1996). 4. Satoh, A., and Kamiyama, S., in ‘‘Continuum Mechanics and Its Applications’’ (G. A. C. Graham and S. K. Malik, Eds.), p. 731. Hemisphere Publishing, New York, 1989.

coida

248

SATOH ET AL.

5. Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Colloid Interface Sci. 178, 620 (1996). 6. Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Magn. Magn. Mater. 154, 183 (1996). 7. Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Colloid Interface Sci. 181, 422 (1996). 8. Bullough, W. (Ed.), ‘‘Proceedings of the Fifth International Conference on ER Fluids, MR Suspensions, and Associated Technology.’’ World Scientific, Singapore, 1996. 9. Kamiyama, S., Koike, K., and Oyama, T., J. Magn. Magn. Mater. 39, 23 (1983). 10. Kamiyama, S., Koike, K., and Wang, Z., JSME Int. J. 30, 761 (1987). 11. Sudou, K., Tomita, Y., Yamane, R., Ishibasihi, Y., and Otowa, H., Bull. JSME 26, 20 (1983). 12. Rosensweig, R. E., ‘‘Ferrohydrodynamics,’’ p. 48. Cambridge Univ. Press, Cambridge, 1985. 13. Russel, W. B., Saville, D. A., and Schowalter, W. R., ‘‘Colloidal Dispersions,’’ p. 30. Cambridge Univ. Press, Cambridge, 1989.

AID

JCIS 5498

/

6g44$$$461

06-10-98 08:54:41

14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27.

Arp, P. A., and Mason, S. G., J. Colloid Interface Sci. 61, 21 (1977). Jeffrey, D. J., and Ohnishi, Y., J. Fluid Mech. 139, 261 (1984). Kim, S., and Mifflin, R. T., Phys. Fluids 28, 2033 (1985). Bossis, G., and Brady J. F., J. Chem. Phys. 80, 5141 (1984). Brady J. F., and Bossis, G., J. Fluid Mech. 155, 105 (1985). Durlofsky, L., Brady J. F., and Bossis, G., J. Fluid Mech. 180, 21 (1987). Brady J. F., Phillips, R. J., Lester, J. C., and Bossis, G., J. Fluid Mech. 195, 257 (1988). Bossis, G., and Brady J. F., J. Chem. Phys. 91, 1866 (1989). Phung, T. N., Brady J. F., and Bossis, G., J. Fluid Mech. 313, 181 (1996). Kamiyama, S., and Satoh, A., ‘‘Microsimulation of Colloidal Dispersions.’’ Asakurashoten, Tokyo, 1997. [In Japanese] Lees, A. W., and Edwards, S. F., J. Phys. C 5, 1921 (1972). Allen, M. P., and Tildesley, D. J., ‘‘Computer Simulation of Liquids,’’ p. 192. Clarendon Press, Oxford, 1987. Odenbach, S., and Sto¨rk, H., J. Magn. Magn. Mater. 183, 188 (1998). Melrose, J. R., Mol. Phys. 76, 635 (1992).

coida