- Email: [email protected]
Development of Effective Stokesian Dynamics Method for Ferromagnetic Colloidal Dispersions (Cluster-based Stokesian Dynamics Method)
Journal of Colloid and Interface Science 255, 98–106 (2002) doi:10.1006/jcis.2002.8606
Development of Effective Stokesian Dynamics Method for Ferromagnetic Colloidal Dispersions (Cluster-based Stokesian Dynamics Method) Akira Satoh1 Department of Machine Intelligence and System Engineering, Faculty of System Science and Technology, Akita Prefectural University, 84-4, Ebinokuchi, Tsuchiya-aza, Honjo 015-0055, Japan Received September 27, 2001; accepted July 15, 2002
physics engineering fields, we are required to clarify in detail the behavior of colloidal particles under conditions of external fields such as a magnetic field and a flow field. It is highly difficult to clarify the behavior of colloidal particles in an applied magnetic or electric field even under a quiescent flow field situation. If another factor of a flow field is added to this situation, this phenomenon becomes much more complicated and is extraordinary difficult to treat analytically. In this case, microsimulation methods such as Stokesian dynamics (SD) and Brownian dynamics come to play an overwhelmingly important role. Although the friction term of particles in a base liquid is sufficient to simulate the particle motion in a highly dilute dispersion, hydrodynamic interactions between particles have a significant influence on the particle motion in a nondilute dispersion. Microsimulation methods that take into account this effect more accurately are indispensable in investigating the aggregation phenomena and rheological properties under conditions of a flow field. There are two kinds of microsimulation methods for colloidal dispersions that have been developed to date: (I) Bossis and Brady’s SD method (7–12), which is based on the additivity of forces, and (II) Satoh’s method (13–16), which is based on the additivity of velocities. The former method is fully established and takes into account the lubrication effect more accurately, but a large amount of computation time is required in simulations because the inverse of the resistance matrix must be calculated for every time step. Thus, they have been attempting to develop an accelerated Stokesian dynamics method for an ordinary colloidal dispersion (17). In contrast, although the latter method is inferior to the former method in the reproduction of the lubrication effect, the calculation of the inverse matrix is unnecessary in simulations, which means that the latter method is applicable to a much larger system than the former method. It is, however, noted that the method based on the additivity of velocities is under development for usual colloidal dispersions, although it is found to be highly useful for ferromagnetic colloidal dispersions (16). The objective of the present study is to develop a new simulation method based on the additivity of forces which enables us to reduce drastically the computation time. We have conducted the three-dimensional simulations of a ferromagnetic colloidal
We have developed a new Stokesian dynamics (SD) method for nondilute colloidal dispersions, which enables us to reduce drastically the computation time. To verify the validity of the present method, which is called the “cluster-based SD method,” threedimensional simulations of a ferromagnetic colloidal dispersion have been carried out for a simple shear flow. The correlation function and viscosity have been evaluated to compare the results obtained by the present method with those obtained by the ordinary SD method and by the method of ignoring hydrodynamic interactions between particles. The results obtained here are summarized as follows. The transient properties from an initial state obtained by the present method agree well with those obtained by the ordinary method, even if a radius rclstr , which defines the cluster formation, is taken as a small value such as rclstr = 1.2d (d is the particle diameter). Also, the equilibrium properties such as the pair correlation function and viscosity obtained by the present cluster-based method are in satisfactory agreement with those obtained by the ordinary SD method. Furthermore, the cluster-based method drastically reduces the computation time to about one-fourteenth to one-seventieth that of the ordinary method. It is clear from these results that the cluster-based SD method is significantly superior to the ordinary SD method for ferromagnetic colloidal dispersions for which a large model system such as N = 1000 or 10,000 is indispensable in simulations. C 2002 Elsevier Science (USA) Key Words: ferromagnetic colloidal dispersion; aggregation phenomena; Stokesian dynamics method; additivity of forces; clusterbased method; viscosity; correlation function.
1. INTRODUCTION
The development of microcontrol techniques, with which colloidal particles dispersed in a base liquid can be microscopically controlled with regard to aggregate structure and orientational distribution, seems to be the key factor for successfully applying functional fluids such as magnetic fluids (or ferrofluids) (1), ER fluids, and MR suspensions (2) to engineering fields (3–6). This kind of technology may be classified as a modern nanotechnology since colloidal particles are generally 10 nm to micrometers in dimensions. To develop microcontrol technologies in colloid 1
Fax: +81-184-27-2188. E-mail: [email protected].
0021-9797/02 $35.00
C 2002 Elsevier Science (USA)
All rights reserved.
98
99
CLUSTER-BASED STOKESIAN DYNAMICS METHOD
dispersion under a simple shear flow situation to verify the validity of the present method. The pair correlation function and viscosity are evaluated to aid in the verification. 2. STOKESIAN DYNAMICS METHODS
2.1. Approximation of Additivity of Forces In nondilute colloidal dispersions, multibody hydrodynamic interactions are a main factor for governing the particle motion. If a simple shear flow is considered, the velocity v i and angular velocity ω i of particle i on the level of the approximation of the additivity of forces are written as (13, 18)
N
I+
(Aii∗
− I) · (v i∗ − U ∗ ) +
j=1(=i)
2 + 3
N
N
Ai∗j · (v ∗j − U ∗ )
j=1(=i) N
B˜ ii∗ · (ω i∗ − Ω∗ ) +
j=1(=i)
∗ ∗ ∗ ˜ Bi j · (ω j − Ω )
j=1(=i)
N 2 ˜ i∗ : E ∗ = F i∗ , − G 3 j=1(=i)
1 2
N
N
Bii∗ · (v i∗ − U ∗ ) +
j=1(=i)
+ I+
N
(C ii∗
[1] Bi∗j · (v ∗j − U ∗ )
j=1(=i)
− I) · (ω i∗ − Ω∗ )
j=1(=i)
+
N j=1(=i)
C i∗j · (ω ∗j − Ω∗ ) −
N
˜ i∗ : E ∗ = T i∗ , H
[2]
j=1(=i)
where the fluid is assumed to flow in the x direction with a shear rate γ˙ , F i and T i are the force and torque acting on the ambient liquid by particle i, (δ x , δ y , δ z ) are unit vectors, I is the unit tensor, Ai∗j , Bi∗j , etc., are the resistance tensors and are referred to in Refs. (18) and (19). The rate-of-strain tensor E has only E x y and E yz as non-zero components; E x y = E yx = γ˙ /2 (E x∗y = ∗ E yx = 1/2). The superscript ∗ denotes dimensionless quantities. In this study Brownian motion of particles is assumed to be negligible. Each quantity is nondimensionalized 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/γ˙ , the velocity is γ˙ a, the force is 6πηa γ˙ a, the angular velocity is γ˙ , and the torque is 8π ηa 3 γ˙ . If Fˆ ∗ describes the column vector containing the forces and torques, vˆ ∗ describes the column vector containing the translational and angular velocities of all N particles, R∗ is the resistance matrix, and Φ∗ is the column vector containing ˜ ∗1 , . . . , H ˜ ∗ including coefficients, the equations for all parG N ticles can be written in a matrix form: Fˆ ∗ = R∗ · vˆ ∗ − Φ∗ : E ∗ .
[3]
Hence, vˆ ∗ can be solved from this equation as vˆ ∗ = R∗−1 · ( Fˆ ∗ + Φ∗ : E ∗ ).
[4]
It is clear from Eq. [4] that the inverse of the resistance matrix must be calculated to obtain the translational and angular velocities. If we consider a three-dimensional system of N particles, this means that a 6N × 6N resistance matrix is treated in simulations. Hence, the simulations based on this approximation are generally limited to a small system, which has already been pointed out. 2.2. Cluster-Based Stokesian Dynamics Method Strong particle–particle interactions due to the lubrication effect arise when two particles are in a nearly touching situation. This effect, therefore, must be taken into account as exactly as possible in developing an approximate method. In a ferromagnetic colloidal dispersion, particles aggregate to form clusters along a magnetic field direction (20–22). It is presumed that the dominant factors for determining the internal structure of aggregates in equilibrium and the process of the cluster formation are mainly the magnetic forces between particles and the hydrodynamic interactions due to the lubrication effect. The drawback of the SD method above-mentioned ordinary comes from the terms Ai∗j · v ∗j and B˜ i∗j · ω ∗j in Eq. [1] and the similar terms in Eq. [2]. Due to these terms, the inverse of 6N ∗ 6N resistance matrix for velocities and angular velocities must be calculated. If it is possible to take into account only the important particles which have a more significant influence on the motion of particle i, smaller resistance matrices are sufficient for determining the particle motion. This means that a drastic reduction in the computation time is accomplished with the more reliable prediction of the particle motion as a first approximation. Inthe present study, we adopt the following approximation for Ai∗j · v ∗j , N j=1(=i)
Ai∗j · (v ∗j − U ∗ ) ≈
N
clstr
Ai∗j · (v ∗j − U ∗ ),
[5]
j=1(=i)
in which the superscript “clstr” for the summation means that only the interactions with the particles belonging to the same cluster ∗ of ∗particle i are taken into account. Similarly, the term B˜ i j · ω j and the similar terms in Eq. [2] are approximated as clstr ˜ ∗ B˜ i∗j · (ω ∗j − Ω∗ ) ≈ Bi j · (ω ∗j − Ω∗ ), j=1(=i) j=1(=i) N N ∗ ∗ ∗ ∗ ∗ clstr ∗ Bi j · (v j − U ) ≈ Bi j · (v j − U ), j=1(=i) j=1(=i) N N ∗ ∗ ∗ ∗ ∗ clstr ∗ C i j · (ω j − Ω ) ≈ C i j · (ω j − Ω ). N
j=1(=i)
N
j=1(=i)
[6]
100
AKIRA SATOH
The physical meanings of these approximations are that the particle motion is simulated by means of exactly taking account of the hydrodynamic interactions between particles in the same cluster and evaluating the inverse of such a matrix. Hence, the lubrication effect is exactly treated in each cluster on the level of the approximation of the additivity of forces. We call the simulation method based on the above-mentioned approximation the “cluster-based SD method.” 2.3. Correction Terms In the above-mentioned method, the contributions of Ai∗j · v ∗j from the particles belonging to the other clusters have been neglected. If we can take into account these contributions even approximately, the accuracy of the method may be improved. Now we treat these contributions in the following way. If the velocity and angular velocity of particle i which are obtained by the p∗ p∗ cluster-based method are denoted by v i and ω i , respectively, ∗ ∗ the correction terms, v i and ω i , are evaluated from equations
N
I+
(Aii∗
− I) · v i∗ +
j=1(=i)
2 + 3
N
N
clstr
Ai∗j · v ∗j
j=1(=i)
B˜ ii∗
·
ω i∗
+
j=1(=i)
N
clstr
B˜ i∗j
·
ω ∗j
j=1(=i)
N 2
p∗
p∗ / ˜∗ Bi j · ω j − Ω∗ , [7] Ai∗j · v j − U ∗ − 3 j=1(=i) j=1(=i) N N clstr ∗ I+ (C ii∗ − I) · ω i∗ + C i j · ω ∗j
=−
N
/
j=1(=i)
1 + 2 =−
N j=1(=i)
j=1(=i)
Bii∗
·
v i∗
+
N clstr
Bi∗j
·
v ∗j
j=1(=i)
N
N
p∗ 1
p∗ / C i∗j · ω j − Ω − Bi∗j · v j − U , [8] 2 j=1(=i) j=1(=i) /
/ in which the notation of means the summation of the contributions from the particles that do not belong to the cluster of particle i. If Eqs. [7] and [8] are solved, the final velocity v i∗ and angular velocity ω i∗ at the next time step can be obtained p∗ p∗ as v i∗ = v i + v i∗ , ω i∗ = ω i + ω i∗ . If we take into account p∗ p∗ the fact that v i and ω i satisfy Eqs. [1] and [2] with the approximations of Eqs. [5] and [6], it is easily understood that Eqs. [7] and [8] are mathematically derived from Eqs. [1] and [2] with the neglect of terms such as / Ai∗j · v ∗j . We call the simulation method with this correction procedure “Method B.” 2.4. Model of Colloidal Dispersion A ferromagnetic colloidal dispersion is considered a model dispersion for investigating the characteristics and differences of the cluster-based SD method and Method B. A particle is idealized as a spherical particle with a central point dipole coated
with a uniform surfactant (or steric) layer. A colloidal dispersion is assumed to be composed of many such particles with the same radius. The magnetic interaction energy between particles (H ) i and j, u i(m) j , the particle–field interaction energy, u i , and the interaction energy arising due to the overlapping of the steric layers, u i(Vj ) , are expressed, respectively, as (20)
3 ds (m) u i(m)∗ kT = λ = u {ni · n j − 3(ni · t i j )(n j · t i j )}, [9] j ij ri j u i(H )∗ = u i(H ) kT = −ξ ni · h, [10]
u i(Vj ) 2ri j /ds ri j /ds − 1 d u i(Vj )∗ = −2 , ln = λV 2 − kT tδ ri j tδ [11] in which ni is the unit vector denoting the direction of the magnetic moment mi (m = |mi |), t i j is the unit vector given by ri j /ri j , ri j = ri − r j , ri j = |ri j |, h = H/H, H is the applied magnetic field (H = |H|), ds is the diameter of the particle excluding the steric layer (the respective radius is expressed by as ), k is Boltzmann’s constant, T is the absolute temperature of fluid, d is the diameter of the particle including the steric layer (a is the radius), tδ is the ratio of the thickness of steric layer δ to the radius of the solid part of the particle, equal to 2δ/ds . The nondimensional parameters, λ, ξ , and λV , appearing in the above equations, are written as λ = µ0 m 2 /4π ds3 kT, ξ = µ0 m H/kT , and λV = π ds2 n s /2; λ, ξ , and λV are the dimensionless parameters representing the strengths of magnetic particle–particle, magnetic particle–field, and steric particle–particle interactions relative to the thermal energy, respectively. In these expressions, µ0 is the permeability of free space, and n s is the number of surfactant molecules per unit area on the particle surface. The force and torque of the particle acting on the ambient fluid, F i∗ and T i∗ (i = 1, 2, . . . , N ), are straightforwardly derived from the energy expressions as F i∗ =
N
N (V )∗ (H )∗ ∗ F i(m)∗ , T + F = T i(m)∗ , [12] i j ij j + Ti
j=1 ( j=i)
j=1 ( j=i)
where F i(m)∗ = −Rm j
F i(Vj )∗
8 [−(ni · n j ) t i j + 5(ni · t i j ) ri∗4j
× (n j · t i j ) t i j − {(n j · t i j )ni + (ni · t i j )n j }], [13]
2 = RV t i j ln ∗ (2/(1 + tδ ) ≤ ri∗j ≤ 2), [14] ri j
T i(m)∗ = −Rm j
2 {ni × n j − 3(n j · t i j )ni × t i j }, ri∗3j
T i(H )∗ = R H ni × h.
[15] [16]
(m) F i(m) j and T i j are the magnetic force and torque acting particle
101
CLUSTER-BASED STOKESIAN DYNAMICS METHOD
i by neighboring particles, T i(H ) is the torque acting particle i due to the deviation of the magnetic moment from an applied magnetic field in direction, and F i(Vj ) is the force due to the overlapping of the steric layers. The nondimensional numbers that appeared in Eqs. [13]–[16], Rm , R H , and RV , are due to the forces and torques being normalized by the viscous force of the shear flow and are written as Rm =
µ0 m 2 µ0 m H kT λV , RH = , RV = . 64π 2 ηs a 6 γ˙ 8πηs a 3 γ˙ 6πηs a 2 γ˙ δ [17]
Rm is the ratio of the representative magnetic particle–particle force to the representative hydrodynamic shear force, R H is the ratio of the representative torque due to an applied field to the hydrodynamic torque, and RV is the ratio of the representative repulsion due to steric layers to the hydrodynamic shear force. Viscosity data are used to discuss the characteristics and differences of the simulation methods. The viscosity increase due to the magnetic forces, ηmyx , is written as (18) ηm∗ yx =
ηmyx ηs
=−
N N N 6π 4π ∗ ∗ y F + T ∗ , [18] V ∗ i=1 j=1 i j i j x V ∗ i=1 i z ( j>i)
in which V ∗ is the dimensionless volume of the system. In this equation, the first term on the right-hand side is due to the magnetostatic forces and repulsive forces between the steric layers of each particle, and the second term is due to the torques acting on the particles including the contribution from an external field. Since we here concentrate our attention on the influence of the magnetic interactions to the viscosity, the contribution due to the stresslets is not taken into account in the present study. 3. COMPUTATION DETAILS AND SIMULATION PARAMETERS
The following conditions and values were used in conducting SD simulations in a simple shear flow situation. The results obtained by the cluster-based method and Method B are compared with those obtained by the ordinary SD method which treats the full resistance matrix. In the ordinary method, huge computation time is required because the inverse of the resistance matrix has to be calculated at every time step. We consider, therefore, a special case of a highly strong magnetic field, which leads to ω i∗ = 0 and ni = h (i = 1, . . . , N ). With this assumption, only translational velocities are evaluated and the resistance matrix becomes 3N × 3N in the ordinary SD method. The number of particles, N , is 64, the volumetric fraction is taken as ϕv = 0.15, and the number density, n ∗ (=na 3 ), is 0.0358. For this case, the length of the simulation box (cube), L ∗ (=L/a), is 12.14. The value of RV /Rm is chosen as 52.89, which corresponds to the case of λ = 9 where thick chain-like clusters are formed in a strong magnetic field (20–22). In our nondimensional sys-
tem, the strength of the shear rate is changed using the nondimensional number Rm since the dimensionless shear rate γ˙ ∗ is always unity. That is, the change of Rm = 0 → ∞ means that the shear rate changes from ∞ to zero. The values of Rm are taken as Rm = 1, 5, 20, 50, and 100 to clarify the influence of the shear rate. The cutoff radius for magnetic particle–particle inter∗ actions, rcoff (=rcoff /a), is L ∗ /2. The time interval, t ∗ , used in the present simulations is as follows: t ∗ = (0.0001, 0.00002) for Rm = (20, 100), for example. The simulations were carried out up to 2,500,000 time steps to obtain averaged values with sufficient accuracy. The resistance tensor, A∗11 , for example, is expressed using A∗ A∗ the resistance functions of X 11 and Y11A∗ as A∗11 = X 11 ee + A∗ Y11 (I − ee), in which e is the unit vector between two partiA∗ cles of interest. The value of X 11 increases rapidly (exactly, in inverse proportion to the particle–particle distance) as particles approach each other to attain a nearly touching situation (16). This type of behavior is undesirable from a simulation point of view as it may lead to the computational instability of the system, arising from the discretization of the translational motion during the simulation. We take into account the fact that the ferromagnetic particles are modeled as a solid sphere coated with a transformable material of surfactant in the present study, and, therefore, we adopted the linear approximation between r ∗ = 2.01 and 2.00 for each resistance function, as in the previous paper (16); that is, each resistance function converges to a value that is for the case of nonhydrodynamic interactions. Also, the screening effect of the third particle on the hydrodynamic interaction between two particles is taken into account in the approximate way which was adopted in the previous paper (16); the detail can be found in the paper. As already pointed out, the strength of the shear rate can be changed using the nondimensional number Rm ; that is, the change of Rm = 0 → ∞ means that the shear rate changes from ∞ to zero. 4. RESULTS
4.1. Transient Properties from an Initial State To clarify the transient properties from an initial state, the evolution of a subaveraged viscosity for 500 time steps with time is shown in Fig. 1 for Rm = 20 and 100. These results are obtained by the cluster-based method (CBM) for two cases ∗ of rclstr /2(=rclstr /a) = (1.2, 1.5), by which a cluster is defined. Also, the results for the ordinary method (OM) and for the approximation of ignoring hydrodynamic interactions (AIHI) are shown for comparison. It can be seen from Fig. 1 that the curves obtained by the CBM agree almost completely with those obtained by the OM and the AIHI in the initial stage. Since an initial configuration of particles is given using simple cubic lattice points, the curves in Fig. 1a are significantly influenced by the initial configuration and seem to oscillate with time in the beginning stage. However, after the initial stage, the subaveraged viscosity comes to change randomly around an average value. This feature of viscosity
102
AKIRA SATOH
creases rapidly once the track of the particle motion deviates from that obtained by the OM, which gives rise to the start of the deviation of the viscosity curve. Even for Rm = 20, however, if the criterion for the cluster formation is relaxed like ∗ rclstr /2 = 1.5, the curve comes to agree well with that obtained by the OM, since the resistance matrices with more particle interactions are dealt with in determining the particle motion. ∗ It is noted that, although we checked the influence of rclstr on ∗ results for five cases of rclstr = (1.05, 1.1, 1.2, 1.5, 1.7) in the present simulations, the significant difference among three cases ∗ of rclstr = (1.05, 1.1, 1.2) could not be found for both the cases of Rm = 20 and 100. Also, it is found that there is no essential difference in results obtained by the CBM and Method B. 4.2. Aggregate Structure in Equilibrium We first show a snapshot of the aggregate structure in equilibrium in Fig. 2 and the pair correlation function in Fig. 3 for
FIG. 1. Transient properties of subaveraged viscosities with time steps for two cases: (a) Rm = 20 and (b) Rm = 100. Each figure has several curves obtained by the CBM, one obtained by the OM, and another obtained by the AIHI.
curves is also valid in Fig. 1b. As the simulation proceeds, the curves for AIHI first start to deviate significantly from those for CBM and OM. This is due to the fact that the lubrication effect, together with magnetic forces, becomes a main factor for governing the particle motion after particles approach each other and aggregate to form clusters. Thus, the simulation based on AIHI cannot capture the particle motion properly in these situations. In contrast, it can be seen from Fig. 1b that the curves ∗ for rclstr /2 = 1.2 and 1.5 are in very good agreement with that obtained by the OM, since the CBM takes into account the lubrication effect exactly on the level of the additivity of forces. Such a good agreement clearly shows that the CBM can simulate the particle motion properly in the process of the cluster formation and is highly available for the simulations of transient phenomena. The agreement with the curve obtained by the OM becomes significant as the value of Rm increases, or, as the influence of magnetostatic interactions increases more significantly ∗ than that of viscous shear forces. The curve for rclstr /2 = 1.2 in Fig. 1a is in much better agreement with that obtained by the OM than that obtained by the AIHI, but it starts to deviate from that obtained by the OM at around 13,000 time steps. Since a colloidal dispersion is a multiparticle system, the difference in-
FIG. 2. Snapshot of aggregate structure obtained by the CBM for Rm = 20; the gradation technique was used to draw the snapshot so that readers can easily recognize the two wall-like clusters formed along the shear flow direction.
CLUSTER-BASED STOKESIAN DYNAMICS METHOD
103
mental result that even in a colloidal dispersion composed of rigid spherical particles without magnetic properties, particles are in a formation similar to the clusters due to the shear forces; such clusters are not due to interparticle forces such as magnetostatic interactions and are in cluster-like formation under the conditions of a shear flow. Although the curves for the different
FIG. 3. Pair correlation function obtained by the CBM for Rm = 20; θ is the angle between the magnetic field and the shear flow directions.
Rm = 20, which are obtained by CBM. Figure 2 was drawn using the gradation technique so that readers can easily recognize the two wall-like clusters formed along the shear flow direction. The definition and evaluation method of the pair correlation function are explained in a standard textbook (23), so we omit showing such a definition and evaluation method. Figure 2 clearly shows that the particles aggregate to form a wall-like structure along the shear flow direction, which has already been obtained in the previous paper (16). This aggregate formation is not essentially different from that obtained by the OM and the AIHI. Since a quantative difference in the internal structure among these simulation methods clearly appears in the curve of the pair correlation function at r ∗ = 1.99 in Fig. 3, we concentrate our attention on this curve. Figure 4 shows the pair correlation functions at r ∗ = 1.99 for the CBM, in which the results are indicated for Rm = 1, 5, and ∗ 20; Figs. 4b and 4c have the curves for four cases of rclstr to ∗ see the influence of values of rclstr . Similar results are shown in Fig. 5 for Method B. In the figures, the notation of θ is defined as the angle from the magnetic field direction toward the shear flow direction. First we consider the results shown in Fig. 4. It ∗ can be seen from Fig. 4c that the curves for rclstr /2 = 1.05, 1.1, and 1.2 almost agree well with each other and have the medium characteristics between the OM and the AIHI. The curve for ∗ rclstr /2 = 1.5 is in good agreement with that obtained by the OM, since the resistance matrices for large clusters are treated in this case. These results clearly show that the CBM can simulate properly the particle motion in equilibrium as well as in transient circumstances. It can be seen from Fig. 4a that, for the case of Rm = 1, both ∗ /2 = 1.2 and 1.5 agree with that obtained by the curves of rclstr OM almost prefectly, which is in significant contrast with the result obtained by the AIHI. The fact that the cluster-based SD method is highly useful even in the situation where the magnetic forces and viscous shear forces are of the same order suggests that the cluster-based SD method is also useful for usual colloidal dispersions. It is expected to a certain degree from the experi-
FIG. 4. Pair correlation functions at r ∗ = 1.99 for three cases: (a) Rm = 1, (b) Rm = 5, and (c) Rm = 20; the results obtained by the CBM are compared with those obtained by the OM and the AIHI.
104
AKIRA SATOH
Figure 5 shows similar results for Method B, in which the particle motion is evaluated with the correction procedure, as already explained. Characteristic features of these results are as follows. For each case of Rm = 1, 5, and 20, the correction terms do not have a significant influence on the particle motion ∗ ∗ for the case of a small value of rclstr such as rclstr /2 = 1.05 and 1.1. These curves are not so different from the corresponding curves in Fig. 4. However, the discrepancy between Method B ∗ and the OM becomes significant as the value of rclstr increases, ∗ which is clearly seen in the curve for rclstr /2 = 1.5 in Fig. 5a, ∗ ∗ for rclstr /2 = 1.2 in Fig. 5b, and for rclstr /2 = 1.5 in Fig. 5c. This incorrectness of the pair correlation function at r ∗ = 1.99 has a tendency to appear more clearly with decreasing values of Rm . Equations [7] and [8] which are satisfied by the correction terms are obtained through mathematical consideration, not through direct reflection of the physical image of the aggregation structure. Hence, it may be understandable that Method B does not give results as good as we expected. The reason why the curve shows more significant error with increasing values of ∗ rclstr is as follows. If the predictors of the velocities and angu∗ lar velocities are sufficiently accurate for small values of rclstr , the correction terms should be small and not give a significant influence on such almost exact values. However, the correction terms usually have errors which come from Eqs. [7] and [8], where terms such as / Ai∗j · v ∗j are neglected. Such errors may spread over many particles through the calculation of the inverse of the resistance matrix to give a significant influence on the particle motion. This means that more particles are influ∗ enced by the errors for the case of larger values of rclstr , since larger resistance matrices are treated in such cases. The particles are in nearly touching formation in the clusters in a strongly interacting system, so that the errors give rise to a serious disturbance on such subtly balanced particle formation. Hence, for ∗ the case of large values of rclstr , such influence are spread over more particles through the calculation of the inverse of the resistance matrices. This explains reasonably that the curves for ∗ ∗ ∗ rclstr /2 = 1.5 in Fig. 5a, rclstr /2 = 1.2 in Fig. 5b, and rclstr /2 = 1.5 in Fig. 5c, are significantly different from those obtained by the OM. 4.3. Averaged Viscosities
FIG. 5. Pair correlation functions at r ∗ = 1.99 for three cases: (a) Rm = 1, (b) Rm = 5, and (c) Rm = 20; the results obtained by Method B are compared with those obtained by the OM and the AIHI.
∗ values of rclstr are not in a good agreement with each other for the case of Rm = 5, as shown in Fig. 4b, there is a tendency that the curves approach that for the OM with increasing values ∗ ∗ of rclstr ; the value at θ = 40◦ for rclstr /2 = 1.5 seems to be too small, but the reason has not been found yet.
Figure 6 shows the influence of the shear rate on the averaged viscosities; Fig. 6a shows the results for the CBM, and Fig. 6b shows the results for Method B. Each figure has four curves ∗ for rclstr /2 = (1.2, 1.5) obtained by the OM and the AIHI. An ∗ error bar (23) is shown for rclstr /2 = 1.2, and the same order of errors is included in the averaged values for the other cases as well. For the case of CBM, shown in Fig. 6a, the curves are in good agreement with that obtained by the OM, especially under the range of Rm = 50. Also, it is not seen that there ∗ /2 = 1.2 is a significant difference between the curves for rclstr and 1.5; a good agreement is observed within the range of errors even for Rm = 100. In contrast, in the case of Method B ∗ shown in Fig. 6b, the results for rclstr /2 = 1.5 are relatively
105
CLUSTER-BASED STOKESIAN DYNAMICS METHOD
TABLE 1 Comparison of CPU times ∗ /2 = rclstr
Absolute values Relative values
Cluster-based Method B Cluster-based Method B
1.05
1.1
1.2
1.5
4 h 18 min 6 h 26 min 1 1.5
7 h 32 min 13 h 35 min 1.75 3.16
22 h 5 min 41 h 55 min 5.14 9.75
76 h 4 min 142 h 15 min 17.7 33.1
larger than those obtained by the OM, and this tendency becomes clearer with decreasing shear rates. This tendency is also ∗ observed for rclstr /2 = 1.2 at Rm = 50, and such a difference at Rm = 50 is seen to be clearly beyond errors. This is easily understood from the previous section’s conclusion that the ∗ increase in values of rclstr gives rise to incorrect aggregate structure in equilibrium. Thus, we conclude that the CBM is highly available and gives good results concerning the viscosity, but Method B gives results with larger errors with decreasing shear rates.
Additivity of forces
Without hydrodynamic interactions
12.5 days
36 min
69.8
0.14
4.4. Comparison of Computation Time Finally, we compare the computation time among the four cases, the CBM, Method B, the OM, and the AIHI; the computation time is one of the important factors for successfully applying microsimulation methods to a large colloidal dispersion. All simulations were conducted on an Alpha workstation with dual CPUs produced by Compaq. The total of about 2,500,000 time steps was needed to obtain the results with sufficient accuracy. Table 1 shows the results of the computation time for Rm = 20; ∗ the values relative to that at rclstr = 1.05 for the CBM are also shown for straightforward comparison. It is seen that a drastic reduction in the computation time can be attained for the CBM compared with the OM. For example, the simulation by CBM ∗ for rclstr /2 = 1.2 requires just one-fourteenth the computation time of the OM. Method B requires twice the computation time of the CBM since the calculation of the inverse of the resistance matrices is necessary twice. It is seen from these results that the drastic reduction in the computation time can be attained if we use CBM. Thus, we conclude that CBM is highly useful as a micro-simulation method for a non-dilute colloidal dispersion. 5. CONCLUSIONS
FIG. 6. Influence of shear rate on averaged viscosities: (a) cluster-based SD method and (b) Method B.
We have developed a new SD method for nondilute colloidal dispersions that enables us to reduce drastically the computation time. To verify the validity of the present method, which is called the “cluster-based SD method,” the three-dimensional simulations of a ferromagnetic colloidal dispersion have been carried out for a simple shear flow. The correlation function and the viscosity were evaluated to compare the results obtained by the present method with those obtained by the ordinary SD method and the method of ignoring hydrodynamic interactions between particles. The results obtained here are summarized as follows. The transient properties from an initial state obtained by the present method agree well with those obtained by the or∗ dinary method, even if the radius rclstr , which defines the cluster ∗ = 1.2d (d is formation, is taken as a small value such as rclstr the particle diameter). Also, the equilibrium properties such as the pair correlation function and the viscosity obtained by the present cluster-based method are in satisfactory agreement with those obtained by the ordinary SD method. Furthermore, the cluster-based method drastically reduces the computation time to about one-fourteenth to one-seventieth that of the ordinary
106
AKIRA SATOH
method. It is clear from these results that the cluster-based SD method is significantly superior to the ordinary SD method for ferromagnetic colloidal dispersions for which a large model system such as N = 1000 or 10000 is indispensable in simulations.
6. 7. 8. 9. 10.
ACKNOWLEDGMENTS
11. 12. 13.
The author gratefully acknowledges Proffessor R. W. Chantrell and Dr. G. N. Coverdale for their valuable advice on this work. This work was partially supported by a Grant-in-Aid for Scientific Research from the Ministry of Education and Science of Japan.
14.
REFERENCES
16. 17. 18.
1. Rosensweig, R. E., “Ferrohydrodynamics.” Cambridge Univ. Press, Cambridge, 1985. 2. Bullough, W. (Ed.), “Proceedings of the Fifth International Conference on ER Fluids, MR Suspensions, and Associated Technology.” World Scientific, Singapore, 1996. 3. Coverdale, G. N., Chantrell, R. W., Hart, A., and Parker, D., J. Magn. Magn. Mater. 120, 210 (1993). 4. Coverdale, G. N., Chantrell, R. W., Satoh, A., and Vietch, R., J. Appl. Phys. 81, 3818 (1997). 5. Coverdale, G. N., Satoh, A., Gilson, R., and Chantrell, R. W., J. Magn. Magn. Mater. 193, 322 (1999).
15.
19. 20. 21. 22. 23.
Satoh, A., J. Colloid Interface Sci. 234, 425 (2001). 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). Satoh, A., Chantrell, R. W., Coverdale, G. N., and Kamiyama, S., J. Colloid Interface Sci. 203, 233 (1998). Satoh, A., Chantrell, R. W., and Coverdale, G. N., J. Colloid Interface Sci. 209, 44 (1999). Satoh, A., Coverdale, G. N., and Chantrell, R. W., J. Colloid Interface Sci. 231, 238 (2000). Satoh, A., J. Colloid Interface Sci. 243, 342 (2001). Sierou, A., and Brady, J. F., J. Fluid Mech. 448, 115 (2001). Kamiyama, S., and Satoh, A., “Microsimulation of Colloidal Dispersions.” Asakurashoten, Tokyo, 1997. [in Japanese] Kim, S., and Karrila, S. J., “Microhydrodynamics.” Butterworth– Heinemann, Boston, 1991. Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Colloid Interface Sci. 178, 620 (1996). Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Magn. Magn. Mater. 154, 183 (1996). Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Colloid Interface Sci. 181, 422 (1996). Allen, M. P., and Tildesley, D. J., “Computer Simulation of Liquids.” Clarendon, Oxford, 1987.
Development of Effective Stokesian Dynamics Method for Ferromagnetic Colloidal Dispersions (Cluster-based Stokesian Dynamics Method) Akira Satoh1 Department of Machine Intelligence and System Engineering, Faculty of System Science and Technology, Akita Prefectural University, 84-4, Ebinokuchi, Tsuchiya-aza, Honjo 015-0055, Japan Received September 27, 2001; accepted July 15, 2002
physics engineering fields, we are required to clarify in detail the behavior of colloidal particles under conditions of external fields such as a magnetic field and a flow field. It is highly difficult to clarify the behavior of colloidal particles in an applied magnetic or electric field even under a quiescent flow field situation. If another factor of a flow field is added to this situation, this phenomenon becomes much more complicated and is extraordinary difficult to treat analytically. In this case, microsimulation methods such as Stokesian dynamics (SD) and Brownian dynamics come to play an overwhelmingly important role. Although the friction term of particles in a base liquid is sufficient to simulate the particle motion in a highly dilute dispersion, hydrodynamic interactions between particles have a significant influence on the particle motion in a nondilute dispersion. Microsimulation methods that take into account this effect more accurately are indispensable in investigating the aggregation phenomena and rheological properties under conditions of a flow field. There are two kinds of microsimulation methods for colloidal dispersions that have been developed to date: (I) Bossis and Brady’s SD method (7–12), which is based on the additivity of forces, and (II) Satoh’s method (13–16), which is based on the additivity of velocities. The former method is fully established and takes into account the lubrication effect more accurately, but a large amount of computation time is required in simulations because the inverse of the resistance matrix must be calculated for every time step. Thus, they have been attempting to develop an accelerated Stokesian dynamics method for an ordinary colloidal dispersion (17). In contrast, although the latter method is inferior to the former method in the reproduction of the lubrication effect, the calculation of the inverse matrix is unnecessary in simulations, which means that the latter method is applicable to a much larger system than the former method. It is, however, noted that the method based on the additivity of velocities is under development for usual colloidal dispersions, although it is found to be highly useful for ferromagnetic colloidal dispersions (16). The objective of the present study is to develop a new simulation method based on the additivity of forces which enables us to reduce drastically the computation time. We have conducted the three-dimensional simulations of a ferromagnetic colloidal
We have developed a new Stokesian dynamics (SD) method for nondilute colloidal dispersions, which enables us to reduce drastically the computation time. To verify the validity of the present method, which is called the “cluster-based SD method,” threedimensional simulations of a ferromagnetic colloidal dispersion have been carried out for a simple shear flow. The correlation function and viscosity have been evaluated to compare the results obtained by the present method with those obtained by the ordinary SD method and by the method of ignoring hydrodynamic interactions between particles. The results obtained here are summarized as follows. The transient properties from an initial state obtained by the present method agree well with those obtained by the ordinary method, even if a radius rclstr , which defines the cluster formation, is taken as a small value such as rclstr = 1.2d (d is the particle diameter). Also, the equilibrium properties such as the pair correlation function and viscosity obtained by the present cluster-based method are in satisfactory agreement with those obtained by the ordinary SD method. Furthermore, the cluster-based method drastically reduces the computation time to about one-fourteenth to one-seventieth that of the ordinary method. It is clear from these results that the cluster-based SD method is significantly superior to the ordinary SD method for ferromagnetic colloidal dispersions for which a large model system such as N = 1000 or 10,000 is indispensable in simulations. C 2002 Elsevier Science (USA) Key Words: ferromagnetic colloidal dispersion; aggregation phenomena; Stokesian dynamics method; additivity of forces; clusterbased method; viscosity; correlation function.
1. INTRODUCTION
The development of microcontrol techniques, with which colloidal particles dispersed in a base liquid can be microscopically controlled with regard to aggregate structure and orientational distribution, seems to be the key factor for successfully applying functional fluids such as magnetic fluids (or ferrofluids) (1), ER fluids, and MR suspensions (2) to engineering fields (3–6). This kind of technology may be classified as a modern nanotechnology since colloidal particles are generally 10 nm to micrometers in dimensions. To develop microcontrol technologies in colloid 1
Fax: +81-184-27-2188. E-mail: [email protected].
0021-9797/02 $35.00
C 2002 Elsevier Science (USA)
All rights reserved.
98
99
CLUSTER-BASED STOKESIAN DYNAMICS METHOD
dispersion under a simple shear flow situation to verify the validity of the present method. The pair correlation function and viscosity are evaluated to aid in the verification. 2. STOKESIAN DYNAMICS METHODS
2.1. Approximation of Additivity of Forces In nondilute colloidal dispersions, multibody hydrodynamic interactions are a main factor for governing the particle motion. If a simple shear flow is considered, the velocity v i and angular velocity ω i of particle i on the level of the approximation of the additivity of forces are written as (13, 18)
N
I+
(Aii∗
− I) · (v i∗ − U ∗ ) +
j=1(=i)
2 + 3
N
N
Ai∗j · (v ∗j − U ∗ )
j=1(=i) N
B˜ ii∗ · (ω i∗ − Ω∗ ) +
j=1(=i)
∗ ∗ ∗ ˜ Bi j · (ω j − Ω )
j=1(=i)
N 2 ˜ i∗ : E ∗ = F i∗ , − G 3 j=1(=i)
1 2
N
N
Bii∗ · (v i∗ − U ∗ ) +
j=1(=i)
+ I+
N
(C ii∗
[1] Bi∗j · (v ∗j − U ∗ )
j=1(=i)
− I) · (ω i∗ − Ω∗ )
j=1(=i)
+
N j=1(=i)
C i∗j · (ω ∗j − Ω∗ ) −
N
˜ i∗ : E ∗ = T i∗ , H
[2]
j=1(=i)
where the fluid is assumed to flow in the x direction with a shear rate γ˙ , F i and T i are the force and torque acting on the ambient liquid by particle i, (δ x , δ y , δ z ) are unit vectors, I is the unit tensor, Ai∗j , Bi∗j , etc., are the resistance tensors and are referred to in Refs. (18) and (19). The rate-of-strain tensor E has only E x y and E yz as non-zero components; E x y = E yx = γ˙ /2 (E x∗y = ∗ E yx = 1/2). The superscript ∗ denotes dimensionless quantities. In this study Brownian motion of particles is assumed to be negligible. Each quantity is nondimensionalized 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/γ˙ , the velocity is γ˙ a, the force is 6πηa γ˙ a, the angular velocity is γ˙ , and the torque is 8π ηa 3 γ˙ . If Fˆ ∗ describes the column vector containing the forces and torques, vˆ ∗ describes the column vector containing the translational and angular velocities of all N particles, R∗ is the resistance matrix, and Φ∗ is the column vector containing ˜ ∗1 , . . . , H ˜ ∗ including coefficients, the equations for all parG N ticles can be written in a matrix form: Fˆ ∗ = R∗ · vˆ ∗ − Φ∗ : E ∗ .
[3]
Hence, vˆ ∗ can be solved from this equation as vˆ ∗ = R∗−1 · ( Fˆ ∗ + Φ∗ : E ∗ ).
[4]
It is clear from Eq. [4] that the inverse of the resistance matrix must be calculated to obtain the translational and angular velocities. If we consider a three-dimensional system of N particles, this means that a 6N × 6N resistance matrix is treated in simulations. Hence, the simulations based on this approximation are generally limited to a small system, which has already been pointed out. 2.2. Cluster-Based Stokesian Dynamics Method Strong particle–particle interactions due to the lubrication effect arise when two particles are in a nearly touching situation. This effect, therefore, must be taken into account as exactly as possible in developing an approximate method. In a ferromagnetic colloidal dispersion, particles aggregate to form clusters along a magnetic field direction (20–22). It is presumed that the dominant factors for determining the internal structure of aggregates in equilibrium and the process of the cluster formation are mainly the magnetic forces between particles and the hydrodynamic interactions due to the lubrication effect. The drawback of the SD method above-mentioned ordinary comes from the terms Ai∗j · v ∗j and B˜ i∗j · ω ∗j in Eq. [1] and the similar terms in Eq. [2]. Due to these terms, the inverse of 6N ∗ 6N resistance matrix for velocities and angular velocities must be calculated. If it is possible to take into account only the important particles which have a more significant influence on the motion of particle i, smaller resistance matrices are sufficient for determining the particle motion. This means that a drastic reduction in the computation time is accomplished with the more reliable prediction of the particle motion as a first approximation. Inthe present study, we adopt the following approximation for Ai∗j · v ∗j , N j=1(=i)
Ai∗j · (v ∗j − U ∗ ) ≈
N
clstr
Ai∗j · (v ∗j − U ∗ ),
[5]
j=1(=i)
in which the superscript “clstr” for the summation means that only the interactions with the particles belonging to the same cluster ∗ of ∗particle i are taken into account. Similarly, the term B˜ i j · ω j and the similar terms in Eq. [2] are approximated as clstr ˜ ∗ B˜ i∗j · (ω ∗j − Ω∗ ) ≈ Bi j · (ω ∗j − Ω∗ ), j=1(=i) j=1(=i) N N ∗ ∗ ∗ ∗ ∗ clstr ∗ Bi j · (v j − U ) ≈ Bi j · (v j − U ), j=1(=i) j=1(=i) N N ∗ ∗ ∗ ∗ ∗ clstr ∗ C i j · (ω j − Ω ) ≈ C i j · (ω j − Ω ). N
j=1(=i)
N
j=1(=i)
[6]
100
AKIRA SATOH
The physical meanings of these approximations are that the particle motion is simulated by means of exactly taking account of the hydrodynamic interactions between particles in the same cluster and evaluating the inverse of such a matrix. Hence, the lubrication effect is exactly treated in each cluster on the level of the approximation of the additivity of forces. We call the simulation method based on the above-mentioned approximation the “cluster-based SD method.” 2.3. Correction Terms In the above-mentioned method, the contributions of Ai∗j · v ∗j from the particles belonging to the other clusters have been neglected. If we can take into account these contributions even approximately, the accuracy of the method may be improved. Now we treat these contributions in the following way. If the velocity and angular velocity of particle i which are obtained by the p∗ p∗ cluster-based method are denoted by v i and ω i , respectively, ∗ ∗ the correction terms, v i and ω i , are evaluated from equations
N
I+
(Aii∗
− I) · v i∗ +
j=1(=i)
2 + 3
N
N
clstr
Ai∗j · v ∗j
j=1(=i)
B˜ ii∗
·
ω i∗
+
j=1(=i)
N
clstr
B˜ i∗j
·
ω ∗j
j=1(=i)
N 2
p∗
p∗ / ˜∗ Bi j · ω j − Ω∗ , [7] Ai∗j · v j − U ∗ − 3 j=1(=i) j=1(=i) N N clstr ∗ I+ (C ii∗ − I) · ω i∗ + C i j · ω ∗j
=−
N
/
j=1(=i)
1 + 2 =−
N j=1(=i)
j=1(=i)
Bii∗
·
v i∗
+
N clstr
Bi∗j
·
v ∗j
j=1(=i)
N
N
p∗ 1
p∗ / C i∗j · ω j − Ω − Bi∗j · v j − U , [8] 2 j=1(=i) j=1(=i) /
/ in which the notation of means the summation of the contributions from the particles that do not belong to the cluster of particle i. If Eqs. [7] and [8] are solved, the final velocity v i∗ and angular velocity ω i∗ at the next time step can be obtained p∗ p∗ as v i∗ = v i + v i∗ , ω i∗ = ω i + ω i∗ . If we take into account p∗ p∗ the fact that v i and ω i satisfy Eqs. [1] and [2] with the approximations of Eqs. [5] and [6], it is easily understood that Eqs. [7] and [8] are mathematically derived from Eqs. [1] and [2] with the neglect of terms such as / Ai∗j · v ∗j . We call the simulation method with this correction procedure “Method B.” 2.4. Model of Colloidal Dispersion A ferromagnetic colloidal dispersion is considered a model dispersion for investigating the characteristics and differences of the cluster-based SD method and Method B. A particle is idealized as a spherical particle with a central point dipole coated
with a uniform surfactant (or steric) layer. A colloidal dispersion is assumed to be composed of many such particles with the same radius. The magnetic interaction energy between particles (H ) i and j, u i(m) j , the particle–field interaction energy, u i , and the interaction energy arising due to the overlapping of the steric layers, u i(Vj ) , are expressed, respectively, as (20)
3 ds (m) u i(m)∗ kT = λ = u {ni · n j − 3(ni · t i j )(n j · t i j )}, [9] j ij ri j u i(H )∗ = u i(H ) kT = −ξ ni · h, [10]
u i(Vj ) 2ri j /ds ri j /ds − 1 d u i(Vj )∗ = −2 , ln = λV 2 − kT tδ ri j tδ [11] in which ni is the unit vector denoting the direction of the magnetic moment mi (m = |mi |), t i j is the unit vector given by ri j /ri j , ri j = ri − r j , ri j = |ri j |, h = H/H, H is the applied magnetic field (H = |H|), ds is the diameter of the particle excluding the steric layer (the respective radius is expressed by as ), k is Boltzmann’s constant, T is the absolute temperature of fluid, d is the diameter of the particle including the steric layer (a is the radius), tδ is the ratio of the thickness of steric layer δ to the radius of the solid part of the particle, equal to 2δ/ds . The nondimensional parameters, λ, ξ , and λV , appearing in the above equations, are written as λ = µ0 m 2 /4π ds3 kT, ξ = µ0 m H/kT , and λV = π ds2 n s /2; λ, ξ , and λV are the dimensionless parameters representing the strengths of magnetic particle–particle, magnetic particle–field, and steric particle–particle interactions relative to the thermal energy, respectively. In these expressions, µ0 is the permeability of free space, and n s is the number of surfactant molecules per unit area on the particle surface. The force and torque of the particle acting on the ambient fluid, F i∗ and T i∗ (i = 1, 2, . . . , N ), are straightforwardly derived from the energy expressions as F i∗ =
N
N (V )∗ (H )∗ ∗ F i(m)∗ , T + F = T i(m)∗ , [12] i j ij j + Ti
j=1 ( j=i)
j=1 ( j=i)
where F i(m)∗ = −Rm j
F i(Vj )∗
8 [−(ni · n j ) t i j + 5(ni · t i j ) ri∗4j
× (n j · t i j ) t i j − {(n j · t i j )ni + (ni · t i j )n j }], [13]
2 = RV t i j ln ∗ (2/(1 + tδ ) ≤ ri∗j ≤ 2), [14] ri j
T i(m)∗ = −Rm j
2 {ni × n j − 3(n j · t i j )ni × t i j }, ri∗3j
T i(H )∗ = R H ni × h.
[15] [16]
(m) F i(m) j and T i j are the magnetic force and torque acting particle
101
CLUSTER-BASED STOKESIAN DYNAMICS METHOD
i by neighboring particles, T i(H ) is the torque acting particle i due to the deviation of the magnetic moment from an applied magnetic field in direction, and F i(Vj ) is the force due to the overlapping of the steric layers. The nondimensional numbers that appeared in Eqs. [13]–[16], Rm , R H , and RV , are due to the forces and torques being normalized by the viscous force of the shear flow and are written as Rm =
µ0 m 2 µ0 m H kT λV , RH = , RV = . 64π 2 ηs a 6 γ˙ 8πηs a 3 γ˙ 6πηs a 2 γ˙ δ [17]
Rm is the ratio of the representative magnetic particle–particle force to the representative hydrodynamic shear force, R H is the ratio of the representative torque due to an applied field to the hydrodynamic torque, and RV is the ratio of the representative repulsion due to steric layers to the hydrodynamic shear force. Viscosity data are used to discuss the characteristics and differences of the simulation methods. The viscosity increase due to the magnetic forces, ηmyx , is written as (18) ηm∗ yx =
ηmyx ηs
=−
N N N 6π 4π ∗ ∗ y F + T ∗ , [18] V ∗ i=1 j=1 i j i j x V ∗ i=1 i z ( j>i)
in which V ∗ is the dimensionless volume of the system. In this equation, the first term on the right-hand side is due to the magnetostatic forces and repulsive forces between the steric layers of each particle, and the second term is due to the torques acting on the particles including the contribution from an external field. Since we here concentrate our attention on the influence of the magnetic interactions to the viscosity, the contribution due to the stresslets is not taken into account in the present study. 3. COMPUTATION DETAILS AND SIMULATION PARAMETERS
The following conditions and values were used in conducting SD simulations in a simple shear flow situation. The results obtained by the cluster-based method and Method B are compared with those obtained by the ordinary SD method which treats the full resistance matrix. In the ordinary method, huge computation time is required because the inverse of the resistance matrix has to be calculated at every time step. We consider, therefore, a special case of a highly strong magnetic field, which leads to ω i∗ = 0 and ni = h (i = 1, . . . , N ). With this assumption, only translational velocities are evaluated and the resistance matrix becomes 3N × 3N in the ordinary SD method. The number of particles, N , is 64, the volumetric fraction is taken as ϕv = 0.15, and the number density, n ∗ (=na 3 ), is 0.0358. For this case, the length of the simulation box (cube), L ∗ (=L/a), is 12.14. The value of RV /Rm is chosen as 52.89, which corresponds to the case of λ = 9 where thick chain-like clusters are formed in a strong magnetic field (20–22). In our nondimensional sys-
tem, the strength of the shear rate is changed using the nondimensional number Rm since the dimensionless shear rate γ˙ ∗ is always unity. That is, the change of Rm = 0 → ∞ means that the shear rate changes from ∞ to zero. The values of Rm are taken as Rm = 1, 5, 20, 50, and 100 to clarify the influence of the shear rate. The cutoff radius for magnetic particle–particle inter∗ actions, rcoff (=rcoff /a), is L ∗ /2. The time interval, t ∗ , used in the present simulations is as follows: t ∗ = (0.0001, 0.00002) for Rm = (20, 100), for example. The simulations were carried out up to 2,500,000 time steps to obtain averaged values with sufficient accuracy. The resistance tensor, A∗11 , for example, is expressed using A∗ A∗ the resistance functions of X 11 and Y11A∗ as A∗11 = X 11 ee + A∗ Y11 (I − ee), in which e is the unit vector between two partiA∗ cles of interest. The value of X 11 increases rapidly (exactly, in inverse proportion to the particle–particle distance) as particles approach each other to attain a nearly touching situation (16). This type of behavior is undesirable from a simulation point of view as it may lead to the computational instability of the system, arising from the discretization of the translational motion during the simulation. We take into account the fact that the ferromagnetic particles are modeled as a solid sphere coated with a transformable material of surfactant in the present study, and, therefore, we adopted the linear approximation between r ∗ = 2.01 and 2.00 for each resistance function, as in the previous paper (16); that is, each resistance function converges to a value that is for the case of nonhydrodynamic interactions. Also, the screening effect of the third particle on the hydrodynamic interaction between two particles is taken into account in the approximate way which was adopted in the previous paper (16); the detail can be found in the paper. As already pointed out, the strength of the shear rate can be changed using the nondimensional number Rm ; that is, the change of Rm = 0 → ∞ means that the shear rate changes from ∞ to zero. 4. RESULTS
4.1. Transient Properties from an Initial State To clarify the transient properties from an initial state, the evolution of a subaveraged viscosity for 500 time steps with time is shown in Fig. 1 for Rm = 20 and 100. These results are obtained by the cluster-based method (CBM) for two cases ∗ of rclstr /2(=rclstr /a) = (1.2, 1.5), by which a cluster is defined. Also, the results for the ordinary method (OM) and for the approximation of ignoring hydrodynamic interactions (AIHI) are shown for comparison. It can be seen from Fig. 1 that the curves obtained by the CBM agree almost completely with those obtained by the OM and the AIHI in the initial stage. Since an initial configuration of particles is given using simple cubic lattice points, the curves in Fig. 1a are significantly influenced by the initial configuration and seem to oscillate with time in the beginning stage. However, after the initial stage, the subaveraged viscosity comes to change randomly around an average value. This feature of viscosity
102
AKIRA SATOH
creases rapidly once the track of the particle motion deviates from that obtained by the OM, which gives rise to the start of the deviation of the viscosity curve. Even for Rm = 20, however, if the criterion for the cluster formation is relaxed like ∗ rclstr /2 = 1.5, the curve comes to agree well with that obtained by the OM, since the resistance matrices with more particle interactions are dealt with in determining the particle motion. ∗ It is noted that, although we checked the influence of rclstr on ∗ results for five cases of rclstr = (1.05, 1.1, 1.2, 1.5, 1.7) in the present simulations, the significant difference among three cases ∗ of rclstr = (1.05, 1.1, 1.2) could not be found for both the cases of Rm = 20 and 100. Also, it is found that there is no essential difference in results obtained by the CBM and Method B. 4.2. Aggregate Structure in Equilibrium We first show a snapshot of the aggregate structure in equilibrium in Fig. 2 and the pair correlation function in Fig. 3 for
FIG. 1. Transient properties of subaveraged viscosities with time steps for two cases: (a) Rm = 20 and (b) Rm = 100. Each figure has several curves obtained by the CBM, one obtained by the OM, and another obtained by the AIHI.
curves is also valid in Fig. 1b. As the simulation proceeds, the curves for AIHI first start to deviate significantly from those for CBM and OM. This is due to the fact that the lubrication effect, together with magnetic forces, becomes a main factor for governing the particle motion after particles approach each other and aggregate to form clusters. Thus, the simulation based on AIHI cannot capture the particle motion properly in these situations. In contrast, it can be seen from Fig. 1b that the curves ∗ for rclstr /2 = 1.2 and 1.5 are in very good agreement with that obtained by the OM, since the CBM takes into account the lubrication effect exactly on the level of the additivity of forces. Such a good agreement clearly shows that the CBM can simulate the particle motion properly in the process of the cluster formation and is highly available for the simulations of transient phenomena. The agreement with the curve obtained by the OM becomes significant as the value of Rm increases, or, as the influence of magnetostatic interactions increases more significantly ∗ than that of viscous shear forces. The curve for rclstr /2 = 1.2 in Fig. 1a is in much better agreement with that obtained by the OM than that obtained by the AIHI, but it starts to deviate from that obtained by the OM at around 13,000 time steps. Since a colloidal dispersion is a multiparticle system, the difference in-
FIG. 2. Snapshot of aggregate structure obtained by the CBM for Rm = 20; the gradation technique was used to draw the snapshot so that readers can easily recognize the two wall-like clusters formed along the shear flow direction.
CLUSTER-BASED STOKESIAN DYNAMICS METHOD
103
mental result that even in a colloidal dispersion composed of rigid spherical particles without magnetic properties, particles are in a formation similar to the clusters due to the shear forces; such clusters are not due to interparticle forces such as magnetostatic interactions and are in cluster-like formation under the conditions of a shear flow. Although the curves for the different
FIG. 3. Pair correlation function obtained by the CBM for Rm = 20; θ is the angle between the magnetic field and the shear flow directions.
Rm = 20, which are obtained by CBM. Figure 2 was drawn using the gradation technique so that readers can easily recognize the two wall-like clusters formed along the shear flow direction. The definition and evaluation method of the pair correlation function are explained in a standard textbook (23), so we omit showing such a definition and evaluation method. Figure 2 clearly shows that the particles aggregate to form a wall-like structure along the shear flow direction, which has already been obtained in the previous paper (16). This aggregate formation is not essentially different from that obtained by the OM and the AIHI. Since a quantative difference in the internal structure among these simulation methods clearly appears in the curve of the pair correlation function at r ∗ = 1.99 in Fig. 3, we concentrate our attention on this curve. Figure 4 shows the pair correlation functions at r ∗ = 1.99 for the CBM, in which the results are indicated for Rm = 1, 5, and ∗ 20; Figs. 4b and 4c have the curves for four cases of rclstr to ∗ see the influence of values of rclstr . Similar results are shown in Fig. 5 for Method B. In the figures, the notation of θ is defined as the angle from the magnetic field direction toward the shear flow direction. First we consider the results shown in Fig. 4. It ∗ can be seen from Fig. 4c that the curves for rclstr /2 = 1.05, 1.1, and 1.2 almost agree well with each other and have the medium characteristics between the OM and the AIHI. The curve for ∗ rclstr /2 = 1.5 is in good agreement with that obtained by the OM, since the resistance matrices for large clusters are treated in this case. These results clearly show that the CBM can simulate properly the particle motion in equilibrium as well as in transient circumstances. It can be seen from Fig. 4a that, for the case of Rm = 1, both ∗ /2 = 1.2 and 1.5 agree with that obtained by the curves of rclstr OM almost prefectly, which is in significant contrast with the result obtained by the AIHI. The fact that the cluster-based SD method is highly useful even in the situation where the magnetic forces and viscous shear forces are of the same order suggests that the cluster-based SD method is also useful for usual colloidal dispersions. It is expected to a certain degree from the experi-
FIG. 4. Pair correlation functions at r ∗ = 1.99 for three cases: (a) Rm = 1, (b) Rm = 5, and (c) Rm = 20; the results obtained by the CBM are compared with those obtained by the OM and the AIHI.
104
AKIRA SATOH
Figure 5 shows similar results for Method B, in which the particle motion is evaluated with the correction procedure, as already explained. Characteristic features of these results are as follows. For each case of Rm = 1, 5, and 20, the correction terms do not have a significant influence on the particle motion ∗ ∗ for the case of a small value of rclstr such as rclstr /2 = 1.05 and 1.1. These curves are not so different from the corresponding curves in Fig. 4. However, the discrepancy between Method B ∗ and the OM becomes significant as the value of rclstr increases, ∗ which is clearly seen in the curve for rclstr /2 = 1.5 in Fig. 5a, ∗ ∗ for rclstr /2 = 1.2 in Fig. 5b, and for rclstr /2 = 1.5 in Fig. 5c. This incorrectness of the pair correlation function at r ∗ = 1.99 has a tendency to appear more clearly with decreasing values of Rm . Equations [7] and [8] which are satisfied by the correction terms are obtained through mathematical consideration, not through direct reflection of the physical image of the aggregation structure. Hence, it may be understandable that Method B does not give results as good as we expected. The reason why the curve shows more significant error with increasing values of ∗ rclstr is as follows. If the predictors of the velocities and angu∗ lar velocities are sufficiently accurate for small values of rclstr , the correction terms should be small and not give a significant influence on such almost exact values. However, the correction terms usually have errors which come from Eqs. [7] and [8], where terms such as / Ai∗j · v ∗j are neglected. Such errors may spread over many particles through the calculation of the inverse of the resistance matrix to give a significant influence on the particle motion. This means that more particles are influ∗ enced by the errors for the case of larger values of rclstr , since larger resistance matrices are treated in such cases. The particles are in nearly touching formation in the clusters in a strongly interacting system, so that the errors give rise to a serious disturbance on such subtly balanced particle formation. Hence, for ∗ the case of large values of rclstr , such influence are spread over more particles through the calculation of the inverse of the resistance matrices. This explains reasonably that the curves for ∗ ∗ ∗ rclstr /2 = 1.5 in Fig. 5a, rclstr /2 = 1.2 in Fig. 5b, and rclstr /2 = 1.5 in Fig. 5c, are significantly different from those obtained by the OM. 4.3. Averaged Viscosities
FIG. 5. Pair correlation functions at r ∗ = 1.99 for three cases: (a) Rm = 1, (b) Rm = 5, and (c) Rm = 20; the results obtained by Method B are compared with those obtained by the OM and the AIHI.
∗ values of rclstr are not in a good agreement with each other for the case of Rm = 5, as shown in Fig. 4b, there is a tendency that the curves approach that for the OM with increasing values ∗ ∗ of rclstr ; the value at θ = 40◦ for rclstr /2 = 1.5 seems to be too small, but the reason has not been found yet.
Figure 6 shows the influence of the shear rate on the averaged viscosities; Fig. 6a shows the results for the CBM, and Fig. 6b shows the results for Method B. Each figure has four curves ∗ for rclstr /2 = (1.2, 1.5) obtained by the OM and the AIHI. An ∗ error bar (23) is shown for rclstr /2 = 1.2, and the same order of errors is included in the averaged values for the other cases as well. For the case of CBM, shown in Fig. 6a, the curves are in good agreement with that obtained by the OM, especially under the range of Rm = 50. Also, it is not seen that there ∗ /2 = 1.2 is a significant difference between the curves for rclstr and 1.5; a good agreement is observed within the range of errors even for Rm = 100. In contrast, in the case of Method B ∗ shown in Fig. 6b, the results for rclstr /2 = 1.5 are relatively
105
CLUSTER-BASED STOKESIAN DYNAMICS METHOD
TABLE 1 Comparison of CPU times ∗ /2 = rclstr
Absolute values Relative values
Cluster-based Method B Cluster-based Method B
1.05
1.1
1.2
1.5
4 h 18 min 6 h 26 min 1 1.5
7 h 32 min 13 h 35 min 1.75 3.16
22 h 5 min 41 h 55 min 5.14 9.75
76 h 4 min 142 h 15 min 17.7 33.1
larger than those obtained by the OM, and this tendency becomes clearer with decreasing shear rates. This tendency is also ∗ observed for rclstr /2 = 1.2 at Rm = 50, and such a difference at Rm = 50 is seen to be clearly beyond errors. This is easily understood from the previous section’s conclusion that the ∗ increase in values of rclstr gives rise to incorrect aggregate structure in equilibrium. Thus, we conclude that the CBM is highly available and gives good results concerning the viscosity, but Method B gives results with larger errors with decreasing shear rates.
Additivity of forces
Without hydrodynamic interactions
12.5 days
36 min
69.8
0.14
4.4. Comparison of Computation Time Finally, we compare the computation time among the four cases, the CBM, Method B, the OM, and the AIHI; the computation time is one of the important factors for successfully applying microsimulation methods to a large colloidal dispersion. All simulations were conducted on an Alpha workstation with dual CPUs produced by Compaq. The total of about 2,500,000 time steps was needed to obtain the results with sufficient accuracy. Table 1 shows the results of the computation time for Rm = 20; ∗ the values relative to that at rclstr = 1.05 for the CBM are also shown for straightforward comparison. It is seen that a drastic reduction in the computation time can be attained for the CBM compared with the OM. For example, the simulation by CBM ∗ for rclstr /2 = 1.2 requires just one-fourteenth the computation time of the OM. Method B requires twice the computation time of the CBM since the calculation of the inverse of the resistance matrices is necessary twice. It is seen from these results that the drastic reduction in the computation time can be attained if we use CBM. Thus, we conclude that CBM is highly useful as a micro-simulation method for a non-dilute colloidal dispersion. 5. CONCLUSIONS
FIG. 6. Influence of shear rate on averaged viscosities: (a) cluster-based SD method and (b) Method B.
We have developed a new SD method for nondilute colloidal dispersions that enables us to reduce drastically the computation time. To verify the validity of the present method, which is called the “cluster-based SD method,” the three-dimensional simulations of a ferromagnetic colloidal dispersion have been carried out for a simple shear flow. The correlation function and the viscosity were evaluated to compare the results obtained by the present method with those obtained by the ordinary SD method and the method of ignoring hydrodynamic interactions between particles. The results obtained here are summarized as follows. The transient properties from an initial state obtained by the present method agree well with those obtained by the or∗ dinary method, even if the radius rclstr , which defines the cluster ∗ = 1.2d (d is formation, is taken as a small value such as rclstr the particle diameter). Also, the equilibrium properties such as the pair correlation function and the viscosity obtained by the present cluster-based method are in satisfactory agreement with those obtained by the ordinary SD method. Furthermore, the cluster-based method drastically reduces the computation time to about one-fourteenth to one-seventieth that of the ordinary
106
AKIRA SATOH
method. It is clear from these results that the cluster-based SD method is significantly superior to the ordinary SD method for ferromagnetic colloidal dispersions for which a large model system such as N = 1000 or 10000 is indispensable in simulations.
6. 7. 8. 9. 10.
ACKNOWLEDGMENTS
11. 12. 13.
The author gratefully acknowledges Proffessor R. W. Chantrell and Dr. G. N. Coverdale for their valuable advice on this work. This work was partially supported by a Grant-in-Aid for Scientific Research from the Ministry of Education and Science of Japan.
14.
REFERENCES
16. 17. 18.
1. Rosensweig, R. E., “Ferrohydrodynamics.” Cambridge Univ. Press, Cambridge, 1985. 2. Bullough, W. (Ed.), “Proceedings of the Fifth International Conference on ER Fluids, MR Suspensions, and Associated Technology.” World Scientific, Singapore, 1996. 3. Coverdale, G. N., Chantrell, R. W., Hart, A., and Parker, D., J. Magn. Magn. Mater. 120, 210 (1993). 4. Coverdale, G. N., Chantrell, R. W., Satoh, A., and Vietch, R., J. Appl. Phys. 81, 3818 (1997). 5. Coverdale, G. N., Satoh, A., Gilson, R., and Chantrell, R. W., J. Magn. Magn. Mater. 193, 322 (1999).
15.
19. 20. 21. 22. 23.
Satoh, A., J. Colloid Interface Sci. 234, 425 (2001). 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). Satoh, A., Chantrell, R. W., Coverdale, G. N., and Kamiyama, S., J. Colloid Interface Sci. 203, 233 (1998). Satoh, A., Chantrell, R. W., and Coverdale, G. N., J. Colloid Interface Sci. 209, 44 (1999). Satoh, A., Coverdale, G. N., and Chantrell, R. W., J. Colloid Interface Sci. 231, 238 (2000). Satoh, A., J. Colloid Interface Sci. 243, 342 (2001). Sierou, A., and Brady, J. F., J. Fluid Mech. 448, 115 (2001). Kamiyama, S., and Satoh, A., “Microsimulation of Colloidal Dispersions.” Asakurashoten, Tokyo, 1997. [in Japanese] Kim, S., and Karrila, S. J., “Microhydrodynamics.” Butterworth– Heinemann, Boston, 1991. Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Colloid Interface Sci. 178, 620 (1996). Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Magn. Magn. Mater. 154, 183 (1996). Satoh, A., Chantrell, R. W., Kamiyama, S., and Coverdale, G. N., J. Colloid Interface Sci. 181, 422 (1996). Allen, M. P., and Tildesley, D. J., “Computer Simulation of Liquids.” Clarendon, Oxford, 1987.