Eldin Wee Chuan Lim

*National University of Singapore Singapore* 

#### **1. Introduction**

214 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

[14] M. Karimi, G. Stapay, T. Kaplan and M. Mostoller, Modelling Simul. Mater. Sci. Eng., 5

[25] Y.O. Ciftci, K.Colakoglu, S.Ozgen and S.Kazanc, J. Phys: Condens. Matter 19 (2007)

[27] W. B.Pearson, Handbook of lattice Spacing and structure of Metals and Alloys

[28] R. O. Simmons and H. Wang, Single Crystal Elastic Constants and Calculated

[31] Y. S. Touloukian and E.H. Buyco, Specific Heat: Metallic Elements and Alloys (New

Aggregate Properties,A Handbook, MIT Press, Cambridge, 1991 [29] Londolt –Bornstein New Series, vols. III-11 and III-18 ,Springer –Verlag, Berlin,1991

[32] J. H. Rose, J.R Smith., F Guinea., J. Ferrante, Phys. Rev. B. 29(6) (1984) 2963.

[37] M.H. Rice, R.G. McQueen , J.H. Walsh, Solid State Physics, 6 (1958) 1.

[36] S.K. Nayak, S.N. Khanna, B.K. Rao, P. Jena, J. Phys. Cond. Mat., 10 (1988) 10853.

[40] Y. Qi, T. Çağn, Y. Kimura, and W.A. Goddard, Phys. Rev. B, 59(5) (1999)3527.

[30] C. Kittel, Introduction to solid state physics, Wiley, New York,1986.

[33] M. Parrinello, A. Rahman, Phys. Rev. Lett. 45(11) (1980) 1196. [34] M. Parrinello, A. Rahman, J. App. Phys. 52 (12) (1981) 7182.

[12] Y.Gurler, S. Ozgen, Matt. Lett. 57 (2003) 4336.

(1997) 337.

326204.

[15] S. Erkoc, Phys. Rep. 278 (1997) 79.

[13] M. W. Finnis, J.E. Sinclair, Philos. Mag., A. 50 (1) (1984) 45.

[16] M.S. Daw, M.I. Baskes, Phys. Rev. Lett., 50 (17) (1983) 1285. [17] M.S. Daw, M.I. Baskes, Phys. Rev. B, 29(12) (1984) 6443.

[26] M.L. Verma, R.P.S. Rathore, Phys. Stat. Sol. b, 185 (1994) 93.

[20] J.S.C. Jang, C.C. Koch, J. Mater Res., 5(3)(1990) 498.

[22] J. Cai and Y.Y. Ye, Phys. Rev.B 54 (12) (1996) 8398.

,Pergamon, Oxford,1967.

York: IFI/ Plenum, 1970)

[35] L. Verlet, Phys. Rev. 159 (1967) 98.

[38] Y. Hiki, A.Granato, Phys. Rev. B, 144 (1966) 411. [39] W.B. Daniels, C.S. Smith, Phys. Rev. B, 111 (1958) 713.

[41] L.Wang, X.B. Fang, H. Li, Mater. Lett. 51(2001) 7.

[18] C. Massobric, V. Pontikis, G. Martin, Phys. Rev. B, 41(15)(1990)10486. [19] F. Cardellini, V. Contini, G.Mazzone, Scripta Metall Mater, 4(1995) 641.

[21] Y.O Ciftci. and K.Colakoğlu, Acta Physica Polanica A, (100)4 (2001) 539.

[23] S. Kazanc, Y.O. Ciftci, K.Colakoglu, and S.Ozgen,Physica B, 381(2006)96. [24] Y.O. Ciftci, K.Colakoglu, S.Ozgen and S.Kazanc, Cent. Eur. J.Phys. 4(2006)472.

The study of magnetic nanoparticles and ferrofluids has gained considerable interests among research workers in recent years. The potential range of application that these novel magnetic nanomaterials can offer is gradually being recognized and continues to be explored. As described in a recent review article (Pamme, 2006), magnetic particles have already been used for such diverse applications as the fabrication of ferrofluidic pumps, solid supports for bioassays, fast DNA hybridization, giant magnetoresistive sensors and superconducting quantum interference devices (SQUID). At a more fundamental level, one of the most important and widely investigated aspects of magnetic nanoparticles and ferrofluids is the formation of self-organized microstructures under the influence of an externally applied magnetic field. A suspension of magnetic nanoparticles in a fluid medium can generally be considered as a single magnetic domain with macroscopic properties that are dependent on the properties of individual nanoparticles as well as the interactions between them (Rosensweig, 1985). In the presence of an external magnetic field, the magnetic domain will be oriented in the direction of the field and may approach saturation magnetization. When the external magnetic field is removed, the domain will revert to a randomly oriented state which exhibits no macroscale magnetism. Although it is well-established that the magnetization of a magnetic fluid or ferrofluid is related to the arrangement of the suspended magnetic nanoparticles, which in turn arises due to the effects of interactions between various types of forces present such as Brownian and dipoledipole interactions for example, current understanding of the kinetics, dynamics and resulting microstructure of the nanoparticle aggregation process is far from complete.

In the research literature, a variety of experimental, theoretical and computational approaches have been applied towards studies of the aggregation and microstructure formation process of magnetic nanoparticles and ferrofluids. In particular, the computational techniques that have been used for such investigations include Monte Carlo simulations (Davis et al., 1999; Richardi et al., 2008), Brownian dynamics (Meriguet et al., 2004, 2005; Yamada and Enomoto, 2008), lattice-Boltzmann method (Xuan et al., 2005), molecular dynamics simulations (Huang et al., 2005), combination of analytical density functional theory and molecular dynamics (Kantorovich et al., 2008), stochastic dynamics (Duncan and Camp, 2006) and analytical methods (Furlani, 2006; Furlani and Ng, 2008; Nandy et al., 2008). Further, several recent studies have also reported comparisons between experimental and theoretical or computational results. For example, the chain formation process of magnetic particles in an external magnetic field and under

Gelation of Magnetic Nanoparticles 217

The molecular dynamics approach to modeling of particulate systems, otherwise known as the Discrete Element Method (DEM), has been applied extensively for studies of flow behaviors in various types of granular and multiphase systems (Lim et al., 2006a, 2006b; Lim and Wang, 2006; Lim et al., 2007; Lim, 2007, 2008, 2009, 2010a, 2010b; Lim et al., 2011). For a comprehensive review, the interested reader is referred to a recent review article by Zhu et al. (2008). The methodology of DEM and its corresponding governing equations have also been presented numerous times in the research literature and only a brief description will be

The translational and rotational motions of individual solid particles are governed by

*i c ij d ij f i dd ij vdw ij e ij Bi ij*

*dv m ff ff f fff dt*

, , , , , , , lub,

1

where mi and vi are the mass and velocity of the ith particle respectively, N is the number of particles in contact with the ith particle, fc,ij and fd,ij are the contact and viscous contact damping forces respectively, ff,i is the fluid drag force that is governed by Stokes' Law, fdd,ij is the dipole-dipole interaction between particles i and j in the presence of an applied magnetic field, fvdw,ij and fe,ij are the van der Waals interaction and electrostatic repulsion between particles i and j respectively, fB,i is the random force arising due to Brownian effects, flub,ij is the lubrication force due to hydrodynamic effects, Ii is the moment of inertia of the ith particle, i is its angular velocity and Tij is the torque arising from contact forces which causes the particle to rotate. The effect of gravity is neglected in the present study.

Contact and damping forces have to be calculated using force-displacement models that relate such forces to the relative positions, velocities and angular velocities of the colliding particles. A linear spring-and-dashpot model is implemented for the calculation of these collision forces. With such a closure, interparticle collisions are modeled as compressions of a perfectly elastic spring while the inelasticities associated with such collisions are modeled by the damping of energy in the dashpot component of the model. The normal (fcn,ij, fdn,ij) and tangential (fct,ij, fdt,ij) components of the contact and damping forces are calculated

> *cn i*, ,, *<sup>j</sup> ni nij <sup>i</sup> f*

*ct i*, ,, *<sup>j</sup> ti tij <sup>i</sup> f* 

*f v nn dn i*, , *<sup>j</sup>* 

*N i i ij j <sup>d</sup> I T dt* 

(1)

(2)

*n* (3)

*t* (4)

*ni r i i* (5)

**2. Computational model 2.1 Discrete element method** 

Newton's laws of motion:

presented here for sake of completeness.

according to the following equations:

1

*j*

*N i*

the effects of shear in microchannels was analyzed and the chain growth rate predicted by the Smoluchowski model was observed to be consistent with experimental observations (Brunet et al., 2005). The transport of an isolated magnetic microsphere or of very dilute suspensions where dipole-dipole interactions are negligible through a microchannel have also been investigated both experimentally and numerically (Sinha et al., 2007). The controlled aggregation of Janus magnetic nanoparticles was studied using dynamic light scattering and cryo-TEM imaging techniques and the main features of the aggregation behavior were consistent with predictions provided by a modified version of the classic Monte Carlo simulation algorithm (Lattuada and Hatton, 2007). Brownian dynamics has also been applied towards the study of motion of magnetic particles in a magnetic field gradient and shown to be in agreement with experimental measurements based on an optical detection method (Schaller et al., 2008).

While most of the investigations of magnetic nanoparticles dispersions have targeted colloidally stable systems, a few studies have also focused on suspensions undergoing aggregation. Colloidal systems undergoing aggregation exhibit complex behaviors due to several factors such as particle-particle interactions, fractal structure of individual clusters and the aggregation mechanism and kinetics. Several modeling approaches have been proposed in the literature to simulate the aggregation kinetics of colloidal systems, either based on Monte-Carlo simulations, or on population balance equations, or on a combination of the cluster mass distribution computed based on the population balance equations with the structure properties of individual clusters determined by Monte-Carlo simulations (Lattuada et al., 2004a). It was found that the average sizes and structure properties predicted in both the diffusion-limited and reaction-limited aggregation regimes were in good agreement with light scattering measurements (Lattuada et al., 2004b). In the case of magnetic nanoparticles aggregation, both experimental and computational studies have underlined substantial differences between diffusion limited aggregation in the absence and in the presence of an applied magnetic field (Tsouris and Scott, 1995; Promislow et al., 1995; Miyazima et al., 1987). In the presence of magnetic fields, clusters grow as chains aligned in the direction of the applied magnetic field. The kinetics of chain growth has been modeled using Monte-Carlo methods (Miyazima et al., 1987), Brownian Dynamics simulations (Dominguez-Garcia et al., 2007), and population balance equations (Martinez-Pedrero et al., 2008). All of these studies have demonstrated how the average size of chains grows as a power law of time, with an exponent that depends upon the particle volume fraction and the strength of the dipolar interactions (Climent et al., 2004). However, no studies to the best of our knowledge have focused on gelation and percolation of magnetic dispersions at high particles volume fractions.

In this work, we report the first application of a modified version of the Discrete Element Method (DEM) towards the simulation of magnetic nanoparticle aggregation with and without an external magnetic field. The various types of interparticle forces that are important in the dynamics of the aggregation process, such as dipole-dipole interactions, van der Waals forces, electrostatic forces and Brownian effects, are taken into account through the incorporation of the respective force models into the classical DEM model. The effects of overall solid fraction and the presence or absence of an external magnetic field on the propensity of such magnetic nanoparticles to aggregate and the microstructure of the resulting clusters or chain-like assemblies formed are investigated computationally.

#### **2. Computational model**

216 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

the effects of shear in microchannels was analyzed and the chain growth rate predicted by the Smoluchowski model was observed to be consistent with experimental observations (Brunet et al., 2005). The transport of an isolated magnetic microsphere or of very dilute suspensions where dipole-dipole interactions are negligible through a microchannel have also been investigated both experimentally and numerically (Sinha et al., 2007). The controlled aggregation of Janus magnetic nanoparticles was studied using dynamic light scattering and cryo-TEM imaging techniques and the main features of the aggregation behavior were consistent with predictions provided by a modified version of the classic Monte Carlo simulation algorithm (Lattuada and Hatton, 2007). Brownian dynamics has also been applied towards the study of motion of magnetic particles in a magnetic field gradient and shown to be in agreement with experimental measurements based on an

While most of the investigations of magnetic nanoparticles dispersions have targeted colloidally stable systems, a few studies have also focused on suspensions undergoing aggregation. Colloidal systems undergoing aggregation exhibit complex behaviors due to several factors such as particle-particle interactions, fractal structure of individual clusters and the aggregation mechanism and kinetics. Several modeling approaches have been proposed in the literature to simulate the aggregation kinetics of colloidal systems, either based on Monte-Carlo simulations, or on population balance equations, or on a combination of the cluster mass distribution computed based on the population balance equations with the structure properties of individual clusters determined by Monte-Carlo simulations (Lattuada et al., 2004a). It was found that the average sizes and structure properties predicted in both the diffusion-limited and reaction-limited aggregation regimes were in good agreement with light scattering measurements (Lattuada et al., 2004b). In the case of magnetic nanoparticles aggregation, both experimental and computational studies have underlined substantial differences between diffusion limited aggregation in the absence and in the presence of an applied magnetic field (Tsouris and Scott, 1995; Promislow et al., 1995; Miyazima et al., 1987). In the presence of magnetic fields, clusters grow as chains aligned in the direction of the applied magnetic field. The kinetics of chain growth has been modeled using Monte-Carlo methods (Miyazima et al., 1987), Brownian Dynamics simulations (Dominguez-Garcia et al., 2007), and population balance equations (Martinez-Pedrero et al., 2008). All of these studies have demonstrated how the average size of chains grows as a power law of time, with an exponent that depends upon the particle volume fraction and the strength of the dipolar interactions (Climent et al., 2004). However, no studies to the best of our knowledge have focused on gelation and percolation of magnetic dispersions at high

In this work, we report the first application of a modified version of the Discrete Element Method (DEM) towards the simulation of magnetic nanoparticle aggregation with and without an external magnetic field. The various types of interparticle forces that are important in the dynamics of the aggregation process, such as dipole-dipole interactions, van der Waals forces, electrostatic forces and Brownian effects, are taken into account through the incorporation of the respective force models into the classical DEM model. The effects of overall solid fraction and the presence or absence of an external magnetic field on the propensity of such magnetic nanoparticles to aggregate and the microstructure of the

resulting clusters or chain-like assemblies formed are investigated computationally.

optical detection method (Schaller et al., 2008).

particles volume fractions.

#### **2.1 Discrete element method**

The molecular dynamics approach to modeling of particulate systems, otherwise known as the Discrete Element Method (DEM), has been applied extensively for studies of flow behaviors in various types of granular and multiphase systems (Lim et al., 2006a, 2006b; Lim and Wang, 2006; Lim et al., 2007; Lim, 2007, 2008, 2009, 2010a, 2010b; Lim et al., 2011). For a comprehensive review, the interested reader is referred to a recent review article by Zhu et al. (2008). The methodology of DEM and its corresponding governing equations have also been presented numerous times in the research literature and only a brief description will be presented here for sake of completeness.

The translational and rotational motions of individual solid particles are governed by Newton's laws of motion:

$$m\_i \frac{dv\_i}{dt} = \sum\_{j=1}^{N} \left( f\_{c, \text{ij}} + f\_{d, \text{ij}} \right) + f\_{f,i} + f\_{dd, \text{ij}} + f\_{vdw, \text{ij}} + f\_{e, \text{ij}} + f\_{\text{B},i} + f\_{\text{turb,ij}} \tag{1}$$

$$I\_i \frac{d o o\_i}{d t} = \sum\_{j=1}^{N} T\_{ij} \tag{2}$$

where mi and vi are the mass and velocity of the ith particle respectively, N is the number of particles in contact with the ith particle, fc,ij and fd,ij are the contact and viscous contact damping forces respectively, ff,i is the fluid drag force that is governed by Stokes' Law, fdd,ij is the dipole-dipole interaction between particles i and j in the presence of an applied magnetic field, fvdw,ij and fe,ij are the van der Waals interaction and electrostatic repulsion between particles i and j respectively, fB,i is the random force arising due to Brownian effects, flub,ij is the lubrication force due to hydrodynamic effects, Ii is the moment of inertia of the ith particle, i is its angular velocity and Tij is the torque arising from contact forces which causes the particle to rotate. The effect of gravity is neglected in the present study.

Contact and damping forces have to be calculated using force-displacement models that relate such forces to the relative positions, velocities and angular velocities of the colliding particles. A linear spring-and-dashpot model is implemented for the calculation of these collision forces. With such a closure, interparticle collisions are modeled as compressions of a perfectly elastic spring while the inelasticities associated with such collisions are modeled by the damping of energy in the dashpot component of the model. The normal (fcn,ij, fdn,ij) and tangential (fct,ij, fdt,ij) components of the contact and damping forces are calculated according to the following equations:

$$f\_{cn,i\dot{j}} = -\left(\kappa\_{n,i}\mathcal{S}\_{n,i\dot{j}}\right)n\_{\dot{i}}\tag{3}$$

$$(f\_{ct,ij} = - (\kappa\_{t,i} \delta\_{t,ij})\mathbf{t}\_i \tag{4}$$

$$f\_{dn,i\dagger} = -\eta\_{n,i} (\upsilon\_r \cdot n\_i) n\_i \tag{5}$$

$$f\_{\text{dft,ij}} = -\eta\_{t,i} \left\{ (\upsilon\_r \cdot t\_i) \mathbf{t}\_i + \left( a \mathbf{o}\_i \times \mathbf{R}\_i - a \mathbf{o}\_j \times \mathbf{R}\_j \right) \right\} \tag{6}$$

Gelation of Magnetic Nanoparticles 219

The electrostatic repulsion between particles due to the so called double layer forces is

2

, nb is the concentration of ions.

*ze h*

1 exp

2

*e ij o i*

exp <sup>2</sup>

 

where is relative permittivity, o is the absolute permittivity of free space, z is the valency of ions, q is the surface charge, e is the fundamental electronic charge, κ-1 is the Debye decay

For the size of nanoparticles simulated, it is also pertinent to consider the random forces arising due to Brownian effects. The algorithm for simulating Brownian forces is similar to

> 2 *i B i f i*

where t is the time step used in the simulation and ni' is a unit vector with a random direction. It is well-established that Brownian effects become less significant at small separation distances between particles due to the presence of hydrodynamic lubrication effects. As such, Brownian forces were set to zero for surface-to-surface distances less than 1

Hydrodynamic interactions due to lubrication effects become important at small surface-tosurface separation distances between particles. The lubrication force between two spheres is described by the lubrication theory and may be calculated as follows (Russel et al., 1989):

Here, as with the calculation of van der Waals forces described earlier, when the surface-tosurface separation distance between two particles is less than 1 nm, hij is fixed at 1 nm to

The simulation conditions applied were based as much as possible on the materials and methods used for experimental studies reported in the literature so that a meaningful comparison between the simulations and experiments can be made. Spherical nanoparticles of diameter 70 nm and density 1000 kg m-3 were simulated within a pseudo-threedimensional computational domain. The dimensions of the computational domain were 5 m 5 m with thickness in the spanwise direction equivalent to one particle diameter. The

6

<sup>2</sup>

16 *r i i ij i ij v n <sup>d</sup> <sup>f</sup> <sup>n</sup> h* 

*<sup>h</sup> kT <sup>f</sup> R q*

 

'

*kT t n* (11)

(12)

*ij*

(10)

*ij*

described by the DLVO theory (Russel et al., 1989):

2 2 2 *o*

 

1

length given by

nm in all simulations.

**2.3 Hydrodynamic interactions** 

avoid the singularity in the above equation.

**3. Simulation conditions** 

,

1 2

that for generating a Gaussian white noise process (Russel et al., 1989):

, 12

lub,

*<sup>d</sup> <sup>f</sup>* 

*b kT ezn* 

where n,i, n,ij, ni, n,i and t,i, t,ij, ti, t,i are the spring constants, displacements between particles, unit vectors and viscous contact damping coefficients in the normal and tangential directions respectively, vr is the relative velocity between particles and Ri and Rj are the radii of particles i and j respectively. If , , *f f ct ij cn ij* tan , then 'slippage' between two contacting surfaces is simulated based on Coulomb-type friction law, i.e. , , *f f ct ij cn ij* tan , where tan is analogous to the coefficient of friction.

#### **2.2 Long range interactions**

In the presence of an externally applied magnetic field, paramagnetic nanoparticles can be considered as magnetic single-domains with a permanent magnetic moment, i proportional to their volume, 3 6 *i i i <sup>d</sup> <sup>M</sup>* , where Mi is the intensity of magnetization. By the superparamagnetic magnetization law for a monodisperse, colloidal ferrofluid (Rosensweig, 1985),

$$\frac{M\_i}{\phi\_s M\_d} = \coth \alpha\_i - \frac{1}{\alpha\_i} \tag{7}$$

where 3 6 *od i i M Hd kT* , o is the magnetic permeability of free space, Md is the saturation magnetization of the bulk magnetic solid, H is the magnetic field strength, k is the Boltzmann's constant, T is thermodynamic temperature and s is the volume fraction of solid present. The anisotropic dipole-dipole interaction energy Edd,ij is then given by (Rosensweig, 1985)

$$E\_{dd,ij} = \frac{1}{4\pi\mu\_o} \left[ \frac{\mu\_i \cdot \mu\_j}{r\_{ij}^3} - \frac{3}{r\_{ij}^5} (\mu\_i \cdot r\_{ij}) (\mu\_j \cdot r\_{ij}) \right] \tag{8}$$

where i and j are the magnetic moments of particles i and j respectively and rij is the displacement vector between the two particles. The dipole-dipole force of interaction acting on particle i is then derived from fdd,ij = -Edd,ij.

The van der Waals forces of interaction and electrostatic repulsion between particles are calculated as follows (Russel et al., 1989):

$$f\_{vdw,ij} = \frac{H\_a}{6} \frac{d\_i^6}{\left(h\_{ij}^2 + 2d\_ih\_{ij}\right)^2 \left(h\_{ij} + d\_i\right)^3} m\_i \tag{9}$$

where Ha is the Hamaker constant and hij is the surface-to-surface separation distance along the line of centers of particles i and j. When the actual surface-to-surface separation distance between two particles is less than 1 nm, hij is fixed at 1 nm to avoid the singularity in the above equation.

, where Mi is the intensity of magnetization. By the

*ti r i i i i j j* (6)

, then 'slippage' between two contacting

(7)

, where tan

(8)

(9)

*f v tt R R dt i*, , *<sup>j</sup>*

where n,i, n,ij, ni, n,i and t,i, t,ij, ti, t,i are the spring constants, displacements between particles, unit vectors and viscous contact damping coefficients in the normal and tangential directions respectively, vr is the relative velocity between particles and Ri and Rj are the radii

In the presence of an externally applied magnetic field, paramagnetic nanoparticles can be considered as magnetic single-domains with a permanent magnetic moment, i proportional

superparamagnetic magnetization law for a monodisperse, colloidal ferrofluid (Rosensweig,

<sup>1</sup> coth *<sup>i</sup> i s d i*

magnetization of the bulk magnetic solid, H is the magnetic field strength, k is the Boltzmann's constant, T is thermodynamic temperature and s is the volume fraction of solid present. The anisotropic dipole-dipole interaction energy Edd,ij is then given by

, 3 5

where i and j are the magnetic moments of particles i and j respectively and rij is the displacement vector between the two particles. The dipole-dipole force of interaction acting

The van der Waals forces of interaction and electrostatic repulsion between particles are

, <sup>2</sup> <sup>3</sup> <sup>6</sup> <sup>2</sup> <sup>2</sup> *a i vdw ij i*

where Ha is the Hamaker constant and hij is the surface-to-surface separation distance along the line of centers of particles i and j. When the actual surface-to-surface separation distance between two particles is less than 1 nm, hij is fixed at 1 nm to avoid the singularity in the

*H d <sup>f</sup> <sup>n</sup>*

*h dh h d*

*ij i ij ij i*

6

 

1 3

 

*i j dd ij i ij j ij o ij ij E r r r r*

4

, o is the magnetic permeability of free space, Md is the saturation

surfaces is simulated based on Coulomb-type friction law, i.e. , , *f f ct ij cn ij* tan

*M M*

of particles i and j respectively. If , , *f f ct ij cn ij* tan

3 6 *i i i <sup>d</sup> <sup>M</sup>* 

is analogous to the coefficient of friction.

3

on particle i is then derived from fdd,ij = -Edd,ij.

calculated as follows (Russel et al., 1989):

*od i*

*M Hd kT*

**2.2 Long range interactions** 

to their volume,

6

 

*i*

(Rosensweig, 1985)

above equation.

1985),

where

The electrostatic repulsion between particles due to the so called double layer forces is described by the DLVO theory (Russel et al., 1989):

$$f\_{e, \dot{\imath}\dot{\jmath}} = 2\pi\varepsilon\varepsilon\_o \left(\frac{kT}{ze}\right)^2 \kappa R\_i q^2 \frac{\exp\left(-\kappa h\_{\dot{\imath}\dot{\jmath}}\right)}{1 - \exp\left(-\kappa h\_{\dot{\imath}\dot{\jmath}}\right)}\tag{10}$$

where is relative permittivity, o is the absolute permittivity of free space, z is the valency of ions, q is the surface charge, e is the fundamental electronic charge, κ-1 is the Debye decay

$$\text{Length given by } \kappa^{-1} = \left(\frac{\varepsilon x\_o kT}{2e^2 z^2 n\_b}\right)^{\frac{1}{2}}, \text{n.s. is the concentration of ions.}$$

For the size of nanoparticles simulated, it is also pertinent to consider the random forces arising due to Brownian effects. The algorithm for simulating Brownian forces is similar to that for generating a Gaussian white noise process (Russel et al., 1989):

$$f\_{B,i} = \sqrt{12\pi\mu\_f \left(\frac{d\_i}{2}\right) kT \Delta t} \stackrel{\cdot}{n\_i} \tag{11}$$

where t is the time step used in the simulation and ni' is a unit vector with a random direction. It is well-established that Brownian effects become less significant at small separation distances between particles due to the presence of hydrodynamic lubrication effects. As such, Brownian forces were set to zero for surface-to-surface distances less than 1 nm in all simulations.

#### **2.3 Hydrodynamic interactions**

Hydrodynamic interactions due to lubrication effects become important at small surface-tosurface separation distances between particles. The lubrication force between two spheres is described by the lubrication theory and may be calculated as follows (Russel et al., 1989):

$$f\_{\rm hub,ij} = \frac{6\pi\mu\left(v\_r \cdot n\_i\right)}{h\_{\rm ij}} \frac{d\_i^2}{16} n\_i \tag{12}$$

Here, as with the calculation of van der Waals forces described earlier, when the surface-tosurface separation distance between two particles is less than 1 nm, hij is fixed at 1 nm to avoid the singularity in the above equation.

#### **3. Simulation conditions**

The simulation conditions applied were based as much as possible on the materials and methods used for experimental studies reported in the literature so that a meaningful comparison between the simulations and experiments can be made. Spherical nanoparticles of diameter 70 nm and density 1000 kg m-3 were simulated within a pseudo-threedimensional computational domain. The dimensions of the computational domain were 5 m 5 m with thickness in the spanwise direction equivalent to one particle diameter. The

Gelation of Magnetic Nanoparticles 221

Fig. 1. Aggregation patterns of magnetic nanoparticles in the absence of an external magnetic field at 10-3 s physical time obtained from the modified DEM simulations. The

solid volume fractions applied were (a) 0.05, (b) 0.10, (c) 0.15 and (d) 0.20.

(a) (b)

(c) (d)

numbers of nanoparticles simulated were 485, 975, 1460 and 1950 which correspond to solid volume fractions of 0.05, 0.10, 0.15 and 0.20 respectively. Here, solid volume fraction is defined to be the ratio of the total volume of all nanoparticles present to the volume of the pseudo-three-dimensional domain. To ensure numerical stability and accuracy, a relatively small time step of 10 ps was applied for all simulations carried out in this study. Table 1 summarizes the values of pertinent material properties and system parameters applied in the simulations. At the start of each simulation, the positions of all nanoparticles were assigned randomly within the computational domain such that no overlap between any two nanoparticles occurred. Periodic boundary conditions were applied on all four sides of the computational domain so as to eliminate any possible effects that may arise due to the presence of boundaries. The application of such boundary conditions also allowed the possibility of simulating a large system using a significantly smaller computational domain which leads to more efficient utilization of computing resources.


Table 1. Material properties and system parameters for DEM simulations
