Dynamic numerical simulations of magnetically interacting suspensions in creeping flow

Powder Technology 279 (2015) 146–165

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

Dynamic numerical simulations of magnetically interacting suspensions in creeping flow R.G. Gontijo, F.R. Cunha ⁎ University of Brasilia, Department of Mechanical Engineering, Fluid Mechanics of Complex Flows Lab — VORTEX, Brasilia, DF 70910 900, Brazil

a r t i c l e

i n f o

Article history: Received 5 February 2015 Received in revised form 11 March 2015 Accepted 21 March 2015 Available online 14 April 2015 Keywords: Magnetic suspensions Sedimentation Magnetization Particle interactions Numerical simulation

a b s t r a c t The equations governing the motion of N magnetic particles suspended in a viscous fluid at low Reynolds and finite Stokes numbers are solved by direct numerical simulations for different Péclet numbers. The Langevin dynamics simulations include all dipole–dipole magnetic interactions for force and torque. An external applied magnetic field and near field interactions represented by contact and repulsion forces are also considered. Repulsive forces are modeled using a variation of the screened-Coulomb type potential. The initial particle distribution is an ergodic ensemble in which each member consists of N mutually impenetrable spheres whose centers are randomly distributed in a prismatic cell of volume V with wall boundaries. The stability of the proposed numerical method and its convergence in calculating some relevant macroscopic properties of the magnetic suspension are carefully examined. The simulations are used to investigate structure transition from an isotropic random distribution of particles to other structures in the presence of an external magnetic field and magnetic particle–particle interactions. The simulations show dimmers and short chain formation in the suspension space at low volume fraction. When the volume fraction is increased long chains and thin anisotropic structures may be observed along the magnetic field direction. The numerical method is also used to calculate the steady state equilibrium magnetization, and accurate results are obtained for different particle volume fractions ϕ in
 agreement with O ϕ3 asymptotic theories. The method presented is able to consider up to 3000 particles with accuracy and computational efficiency. Typical configurations and particle trajectories are also shown and discussed from a physical point of view. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Particulate systems formed by solid particles immersed in Newtonian liquids appear in many practical situations in civil, chemical and oil engineering. Many theoretical, experimental and numerical works have been done in order to investigate the behavior of such systems. Usually the authors explore the mean sedimentation velocity and velocity fluctuations for particulate systems with hydrodynamic interactions. In this context Davis and Acrivos [1] have provided a deep discussion on the sedimentation motion of noncolloidal monodisperse and polydisperse spherical particles at low Reynolds number. The authors also discussed ways to enhance the sedimentation velocity by adopting inclined channels known as the Boycott effect [2]. In several cases sedimentation is used to separate the particulate phase and may be extremely important to accelerate this process. Another important result obtained in the study of particle sedimentation at low Reynolds number is the empirical relation of Richardson–Zaki [3] which states that the ⁎ Corresponding author. Tel.: +55 6133072314; fax: +55 6133073259. E-mail addresses: [email protected] (R.G. Gontijo), [email protected] (F.R. Cunha).

http://dx.doi.org/10.1016/j.powtec.2015.03.033 0032-5910/© 2015 Elsevier B.V. All rights reserved.

mean sedimentation velocity 〈U〉 of a particulate system with hydrodynamically interacting particles is given by 〈U〉 = Us(1 − ϕ)n, where U s ¼ ð2=9ηÞa2 Δρg denotes the Stokes velocity of an isolated particle, a is the particle radius, Δρ represents the density difference between the particle and the fluid, μ is the fluid viscosity, g is the intensity of gravity acceleration, ϕ represents the particle volume fraction and n = 5.1 for rigid spherical particles sedimenting at low Reynolds number. Regarding the physics of hydrodynamically interacting suspensions, we still have some interesting open questions. For instance the computation of hydrodynamic interactions by using Green's functions of incompressible creeping flows [4] leads to a well known divergence problem where some important statistical properties such as the mean velocity and the mean velocity variance depends on the size of the system. This problem was first investigated theoretically by [5–7]. While the determination of the sedimentation velocity was renormalized by Batchelor's calculations, the theoretical calculation of the velocity variance in sedimentation seems to be still an open question in suspension mechanics [8]. This divergence, however, is not observed in experiments carried out by [9,10]. This problem is directly related with the long range nature of the fluid velocity disturbances due to particle–particle hydrodynamic interactions that one sphere

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

produces on the surrounding ones by viscous diffusion in creeping flows. The velocity disturbance produced by a particle in sedimentation decays slowly like (1/r), where r denotes the distance from the particle and an arbitrary point within the fluid. An excellent review on this problem covering recent theoretical, numerical and experimental works was done by Guazzelli and Hinch [8]. We will investigate the velocity fluctuations in a magnetic suspension produced by dipole–dipole interactions. The disturbances caused by a particle on which a net external magnetic torque and force act decay as 1/r3 and 1/r 4, respectively. Regarding practical applications of magnetic suspensions, we may say that they started to be explored by physicists and engineers due to its wide range of applications. These suspensions are basically made by adding small particles of a ferromagnetic material into a carrier fluid, usually a Newtonian liquid, that can be water, ester or oil based. Magnetic suspensions consisting of nanometric particles are denominated in the current literature as ferrofluids whereas suspensions formed by micrometric sized particles are called magnetorheological suspensions [11]. These suspensions are used in several areas. Bashtovoi et al. [12] have studied the use of magnetic fluids for capillary ascension. They found that the fluid meniscus deformation when subjected to an applied magnetic field is responsible for pumping the fluid inside a capillary tube without the need of an external device. Larachi and Munteanu [13] investigated different ways of changing the effect of gravitational force by a controlled magnetic field for liquid and gas flows of magnetic fluids in porous media. In the biomedicine area, Brusentsov et al. [14] have explored the use of a bio-compatible magnetic fluid for the treatment of tumors by a technique known as magnetohydrodynamic thermochemotherapy. Recently, Gontijo and Cunha [15] have carried out experiments and performed theoretical calculations on thermomagnetic convection for increasing the efficiency of cooling systems. Recent advances in magnetic fluids rheology and hydrodynamic, including a theoretical part on the governing equation for magnetic fluids and experimental observations of these complex fluids have been reviewed by several authors [16,17]. The primary focus on early development stages of ferrohydrodynamics was to control the motion of magnetic fluids by applying an external field. However, as more studies were done on this topic, interesting features naturally arose. For instance, the magnetoviscous effect has been discovered, this mechanism corresponds to an increase in viscosity of ferrofluids under application of an external magnetic field. Theoretical studies on this problem were done by [25–27] to describe the flow of a magnetic fluid with viscosity variation during the flow due to an applied external magnetic field. Some controversies between the closure problems of the governing equations for modeling the flow of non-symmetric magnetic fluids are noticed in these works and in the current literature [11,16,17]. The first experimental observations in the field of magnetic field assisted fluidization, a general dimensional analysis and scaling study have been successfully conducted by [18,19]. A theoretical work on the stability of magnetic fluidized bed suspensions was developed recently by [20]. Some authors [22–24] have also used interesting computational techniques, such as Lattice–Boltzmann and control volume based finite element method, to explore the flow field of magnetic fluids under the presence of external fields. Since these suspensions have been extensively explored in several applications it is important to understand their complex physical behavior by direct simulations in order to propose constitutive equations for a continuum description of these fluids as the stress tensor and the magnetization evolution equation. This is a controversial point and is still in progress when we compare, for instance, the formulation presented in a paper by [28] on the stress tensor for a magnetic fluid with what is currently proposed by [11,17,16]. Direct dynamic simulations of a collection of magnetically interacting particles throughout the suspension volume with periodic boundary conditions certainly is a promising way for understanding the physics of magnetic suspensions, in particular their microstructural behavior

147

and how it affects the stability of such a suspension and its mean transport properties. In previous works we have examined the behavior of dilute magnetic suspensions with hydrodynamic and magnetic interactions based on the relative trajectory analysis of two interacting particles. The effect of particle rotation motion and particle inertia on these interactions was explored. Particle trajectories, particle aggregation and particle self-diffusion induced by particle–particle hydrodynamic and magnetic interactions have been also investigated [20,21,29–32]. Numerical codes performing many interacting bodies simulation are usually based on Langevin dynamics (LD) or non-equilibrium molecular dynamics (NEMD), static Monte-Carlo (MC) simulations and on a more theoretical approach describing a probability density function of a suspension quantity (e.g. dipole orientation) which requires the solution of a Fokker–Planck nonlinear differential equation as part of the problem [17]. The work is focused on using Langevin dynamics simulations to model the dynamic behavior of magnetically interacting suspensions. This method basically solves for each particle suspended in a medium its momentum equation based on the Langevin stochastic equation that comes from Newton's second law of motion applied to colloidal particles, subjected to Brownian motion. The forces and torques are computed in each time-step of the simulation. Several works have used LD or similar approaches such as lattice Boltzmann simulations successfully to simulate sedimenting particles interacting hydrodynamically in monodisperse and polydisperse suspensions [17,33–40]. However, just few works in the current literature have investigated by direct dynamic simulations velocity fluctuations, transport properties and calculation of bulk magnetization in Brownian and non-Brownian magnetic suspensions of magnetically interacting particles [17]. Nowadays, there is a significant lack of numerical studies concerning important computational details on the simulation of magnetic suspensions with particle inertia and rotation in systems with periodic and non-periodic boundary conditions. With the use of periodic boundary conditions, the summations of particle–particle interactions, O(N 2) operations, become lattices sums and the accurate, but expensive technique proposed by Ewald [41] should be applied to accelerate convergence. We show in this work that it is possible to save a significant computational time by avoiding a full calculation of particle interactions with periodic boundary conditions for suspensions with volume fraction around 5%. In this limit we use lattice sums in a hybrid way, calculating these periodic sums only for magnetic particle interactions induced by torque-dipole–dipole interactions which decay like 1/r 3. As the dipole–dipole interactions due to forces decay faster like 1/r 4 this hybrid approximation gives accurate results and the method becomes very attractive computationally to simulate magnetic suspensions that in general have an average particle volume fraction below 5%. This paper also aims to examine the behavior of particle velocity fluctuations and structure formation in magnetic suspensions. The accuracy of our LD simulations is illustrated by comparing the numerical prediction of the equilibrium magnetization with known theoretical results for the behavior of this property as a function of particle volume fraction proposed by [28,42].

2. Formulation of the problem Consider a suspension of N magnetic rigid, spherical and monodisperse particles with diameter 2a and density ρs immersed in a fluid with an absolute temperature T with viscosity η and density ρf. The location of each sphere's center is described by a set of vectors rN and their angular displacement by the vector di. The spheres have a permanent magnetic dipole moment mi = mdi, (i = 1,…,N) where m is the dipole intensity and di is a unit vector that gives the particle dipole orientation. The linear and angular velocities of all spheres is defined by a set of vectors uN and ωN respectively. Thus, a configuration of the magnetic suspension is completely specified by the object C = {r1, …, rn, u1, … un; d1, …, dn, ω1, … ωn} (see Fig. 1).

148

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

interactions promote particle rotation and consequently the rotation of their fixed dipole moments. We also consider rotational inertia and extra torques imposed by an external magnetic field. Under these conditions, the dimensional governing equation for the translational motion of a particle i inside a monodisperse suspension of N particles is given by

Ms

dui 4 3 i i i i i ¼ −6πηau þ Δρg πa þ f B þ f m þ f r þ f c ; dt 3

ð1Þ

where Ms is the mass of an arbitrary particle i, ui denotes the particle velocity, t is the time, a is the particle radius, Δρ = ρs − ρf represents the density difference between the carrier fluid and the particles and g denotes the gravitational acceleration vector. Here, fBi, fmi, fri and fci are the non-hydrodynamic forces exerted on the particle i, namely Brownian, magnetic, repulsive and contact forces, respectively. Now, the rotational motion of the particles is described by Newton's second law for the balance of torques on the particle i. We have

Fig. 1. A sketch of the problem, illustrating a typical magnetic suspension simulated in this work. In this figure we may notice that the dipole moments of the particles are initially distributed randomly. The initial position of the particles is also set randomly. This picture illustrates only one configuration. In our notation gravity will always be parallel to the z axis as shown in this sketch and the applied field, when applied, will also be in the same direction of gravity.

The particles interact magnetically through their permanent magnetic dipole moments and also with an external applied magnetic field. The particulate system is subjected to a net gravitational force, due to the difference between the particle and the fluid densities, ρs and ρf respectively. Hydrodynamic drag is also considered in the regime of Re ≪ 1, where Re = ρfaUs/η denotes the particle Reynolds number. It is important to notice that in this work the only hydrodynamic force acting on the particles is the viscous drag. This drag is purely produced by the motion of the particles that occurs in a low Reynolds number regime. Since we consider particle inertia, the influence of hydrodynamic interactions is smaller than it would be for the approaching limit of nonmassive particles. It is also worthy to notice that the parameter responsible for measuring the influence of hydrodynamic interactions is the volume fraction of particles. In this work we will mostly work in diluted regimes so this approximation of null hydrodynamic interactions will be shown to be in good agreement with our theoretical validation base. Other hydrodynamic forces such as virtual mass, Oseen [43] and Basset [44] drags are not considered. Since we are considering Re ≪ 1 and non-null Stokes number, as shown by Sobral et al. [45] the dominant hydrodynamic force on the particles is the Stokes hydrodynamic drag. The size of the particles can be nanometric or micrometric. Therefore, both Brownian and non-Brownian magnetic suspensions may be simulated in this work. Particle inertia and rotation are also considered. It is assumed that the dipole moment m is fixed to the particle, and so it rotates with the particle angular velocity with no delay. Physically, we are assuming the condition that the Néel relaxation time τN for the magnetic moment inside the particle to align with H without particle rotation, is much larger than the Brownian rotational relaxation time τB. So, on this time scale a particle needs to rotate in the surrounding fluid in order to align m with the field H. In addition, under this condition the vorticity of the flow can be much larger than the intrinsic angular velocity of the dipole moment. A typical time scale of the flow like τflow = a/Us it is also very small compared with τN, but it could be of the same order or much smaller than τB. The ratio between τB/τflow is defined as the Péclet number, Pe. Thus, when Pe ≪ 1 the condition of a Brownian magnetic suspension is assumed, whereas for Pe ≫ 1 Brownian motion becomes unimportant on the particle motion. Brownian and magnetic torques due to the magnetic dipole moment

i

Js

dω 3 i i i ¼ −8πηa ω þ T B þ T m ; dt

ð2Þ

where Js is the polar moment of inertia of the particle, ωi is the particle i angular velocity, TBi denotes Brownian torques and Tmi represents the magnetic torque acting on the particle i. Solving Eq. (2) is very important to properly represent the dynamics of magnetic suspensions as shown in a recent work by [20]. Relevant physical properties, such as the magnetization of a magnetic suspension are computed based on the time evolution process of the particle dipole orientation. In the next section, we present the models for all non-hydrodynamic forces fmi, fBi, fri, and fci and also non-hydrodynamic torques TBi and Tmi. 2.1. Magnetic particle interactions: fmi and Tmi The potentials representing magnetic interactions in the dipolar matter, considering both dipole–dipole interactions and dipole– external field interactions are mathematically expressed respectively as [17] ψi j ¼

X μ 0 mi m j h i≠ j

4πr 3i j



i di  d j −3 di  ^ri j d j  ^ri j ;


 ^ ; ψi ¼ −μ 0 mi H di  h

ð3Þ

ð4Þ

where the lower indexes i and j denote particles i and j respectively, μ 0 is the magnetic permeability of the free space (4 π × 10−7 H ⋅ m−1), mi and mj are the magnitude of the magnetic dipole moments of particles i and j with directions di and d j respectively, rij is the distance between the centers of the particles i and j, ^r i j is the unitary vector in the direction that links particle i to particle j, H is the intensity of an external magnetic field applied on the particles with direction ĥ. The magnetic force acting on a particle i of the suspension is given by
 i f m ¼ − ∇ψi j þ ∇ψi ;

ð5Þ

and the magnetic torque on a particle i Tmi is calculated as
 i T m ¼ −di  ∇di ψi j þ ∇di ψi ;

ð6Þ

where ∇di is a vector operator that accounts the derivatives with respect to the orientation of the magnetic dipole moment of the ith particle. The

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

calculation of the terms defined in Eqs. (5) and (6) results in the following expressions for fmi and Tmi, respectively [17] i

fm ¼

X 3μ 0 mi m j h
i≠ j

4πr4i j





 i di  d j r^i j þ di  ^ ri j d j þ d j  ^ ri j di −5 di  ^ ri j d j  ^ri j ^ri j

þ μ 0 mi di  ∇H;

i

Tm



 X 3μ 0 mi m j  1
¼ − di  d j þ d j  ^ri j di  ^r i j 3 3 4πr i j i≠ j h
i ^ þ μ 0 mi H di  h :

ð7Þ

ð8Þ

will be a particle relaxation time defined by the instant between the time a particle receives the influence of an external force and the time for a particle responses to this perturbation. In particular, when a magnetic suspension is simulated, attractive forces will appear along the time evolution of the suspension. If particle inertia is considered, the repulsion force given by Eq. (9) is not strong enough to prevent particle overlap. This indicates that a particle can collide with another and also with the solid boundaries. These contact forces are modeled here using a Hertz tension that depends on the mechanical properties of the surfaces in contact. The contact force acting on the particle i is modeled by i

2.1.1. Near field contact forces: fri and fc Since hydrodynamic interactions are not considered in the present work, there are no lubrication forces acting on near field-contact particles. In order to simulate a repulsive force that emulates lubrication forces on near contact particles and in those particles in the neighborhoods of a solid wall, a numerical fictitious force is defined as [46]
i fr

¼ C 1 ð6πηaÞui e



−Ci j

2

 ^er :

ð9Þ

A scheme illustrating the gap for near-field interaction between two particles and between a particle and the wall with the quantities used in Eq. (9) is shown in Fig. 2. Here C1 and C2 are calibration constants related to the intensity and range of the near field interaction, respectively. These constants are examined during the numerical simulations of two approaching spheres taking into account magnetic interactions between the particles [56]. The variable ui denotes the magnitude of the ith particle velocity, ij is the distance between the surface of two near-contact particles or differently between a solid wall and the surface of a particle in its neighborhoods, given by ij = rij − 2a (particle–particle) or ij = riw − a, rij denotes the distance between the centers of the particles, riw represents the distance between the center of a closing particle to a wall and êr is the unit vector in the direction of the repulsion and depends on the configuration of the particulate system only. This approach was proposed by [46] and successfully adopted by [57] when investigating the motion of a concentrated sedimenting blob of rigid particles at low Reynolds number. A similar repulsive force known as the screened-Coulomb type was used in the works developed by [47]. The use of this repulsive force approximation is much cheaper computationally than the numerical implementation of lubrication forces by means of considering the resistance matrix of hydrodynamic interactions. In this sense, this approach does not require resistance matrix inversion at each small time-steps in order to capture the near field particle motion as commonly treated in Stokesian dynamics simulation [48]. In Eq. (1) a contact force fc is also required by the simulations with particle inertia. Since this formulation considers particle inertia, there

Repulsion direction Motion direction

(a)

(b)

Fig. 2. A simple sketch of the near field repulsion forces. (a) The gap used to account for repulsive forces between two particles. b corresponds to the gap used for modeling the wall effect as a repulsive force between a particle and wall. Both situations consider a variation of the screened-Coulomb type potential that represents these near touching solid boundaries as an artificial lubrication force.

1=2 3=2 i j ^er ;

f c ¼ C 3 εb

ð10Þ

where C3 is a calibration constant of the model and  is a material constant defined in terms of the particle elastic modulus; the Young modulus E and the Poison coefficient ν. The quantity  is given by 1 ð1−ν Þ ¼2 : ε E

ð11Þ

In Eq. (11) ν denotes the Poisson modulus of the spheres and E is the Young modulus of the particles. The parameter b that appears in Eq. (10), for a monodisperse suspension is the half of the particle radius. The calculation of the parameters ε and b changes for polydisperse suspensions having particles with different size and density. Fig. 3 presents a schematic showing a typical configuration of particles overlapping. The distance ij is also shown in the figure. 2.1.2. Brownian force and torque: fBi and TBi The Brownian force is the last non-hydrodynamic force from Eq. (1) to be modeled. The model for this stochastic force is well known from literature [49] and it is given by the solution of the stochastic Langevin's differential equation, assuming that Brownian fluctuations are isotropic, have no memory and using the dissipation–fluctuation theorem. Under this condition a random force fb(t) is considered a stationary white noise with the random fluctuations non-correlated on the time scale of the particle motion. Mathematically, this condition is expressed in the form   0   0 h f b ðt Þi ¼ 0; f b ðt Þf b t ¼ F B δ t−t ;

ð12Þ

where δ is the Dirac delta distribution and FB = 12πηakBTδ relates the intensity of the Brownian force with viscous force which dissipates the random thermal fluctuations of the particles [49]. Here, kB denotes the Boltzmann constant. Now taking the trace of Eq. (12) it leads to the fluctuation and dissipation theorem for force 

Repulsion direction

149

 0   0 f b ðt Þ  f b t ¼ ð6πηaÞð6kB T Þδ t−t ;

ð13Þ

Repulsion direction

Motion direction

Fig. 3. A schematic of the contact force under the condition of particle overlap. In order to generate a repulsive force we introduce a non-linear (scales with 3/2 ij ) Hertz contact force.

150

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

by the same way for the Brownian torque we found 

  0 
 3 0 T b ðt Þ  T b t ¼ 8π ηa ð6kB T Þδ t−t :

i

Tm ¼

 6Dt 1=2 ξ; δτ

TB ¼ 8π ηa

3

 6Dr 1=2 ξ: δτ

ð15Þ

ð16Þ

Here Dt = kBT/(6π ηa) and Dr = kBT/(8π ηa3) are, respectively the translational and rotational Brownian diffusion coefficients of Stokes– Einstein [50]. In the present context, δτ is a typical time-step related to Brownian thermal fluctuations, ξ represents a random unitary vector associated with the modeling of Brownian forces and torques given by ξ¼

 1
ξx ^ex þ ξy ^ey þ ξz ^ez jξ j

ð17Þ

where ξx, ξy and ξz are three different random numbers in the interval qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi [−1,+1] and jξ j ¼ ξ2x þ ξ2y þ ξ2z . 2.1.3. Dimensionless equations and physical parameters In order to make the governing equations non-dimensionalized, we define the following dimensionless quantities by using a, Us as the typical length and velocity scales u  tU s ωa  and ω ¼ ; t ¼ ; a Us Us



u ¼

ð18Þ

where U s ¼ ð2=9ηÞΔρga2 is the Stokes velocity of the particle. From this point the upper index * will be suppressed to make the notation as simpler as possible. Now governing Eqs. (1) and (2) describing the motion of a particle i in the suspension are written in terms of these dimensionless quantities as follows i

St

du i ^ þ F mi þ f Bi þ f ri þ f ci ¼ −u þ g dt

ð19Þ

and St r

dωi i i i i ¼ −ω þ T m þ T B þ T m : dt

ð20Þ

The non-hydrodynamic forces fmi, fBi, fri, and fci and torques Tmi and TBi acting on the particle i are given as follows: i

Fm ¼





 i 24λ 1 X h
di  d j ^r i j þ di  ^r i j d j þ d j  ^r i j di −5 di  ^r i j d j  ^r i j ^ri j Pe r 4i j i≠ j hα i d  ∇h ; þ Pe i

ð21Þ i

fB ¼



6 Peδτ

1=2


i

f r ¼ C 1 ui e

i j

−C

i 3=2 f c ¼ P c i j ^er

2

ξ

ð22Þ

 ^ er

ð25Þ

ð14Þ

Using Eqs. (13) and (14) we can model the stochastic force fB(t) and stochastic torque TB, respectively, as being f B ¼ 6πηa

 

  α
 24λ 1 X 1
^ ; − di  d j þ d j  ^r i j di  ^r i j þ di  h 3 Per ri j i≠ j 3 Per

i

6 Per δτ

1=2

ξ:

ð26Þ

Here St = MsUs/(6πηa2) is the Stokes number related with the particle inertia, that represents the ratio between the particle relaxation time and the convective time scale of the particle. Pe = Usa/Dt is the Péclet number, it expresses the ratio between the Brownian diffusive time scale and the convective time scale. Controlling we can simulate Brownian and non-Brownian magnetic suspensions. In addition, Str = JsUs/(8πμ a4) is the equivalent rotational Stokes number and Per = Us/Dr a is the rotational Péclet number. Pc = Cεa2/(6πμaUs) is the contact parameter, it measures the relative importance between the contact force and the hydrodynamic drag force. As the particle has a higher stiffness, the contact time will decrease and the parameter Pc increases. In all the simulations of this work we have used Pc = 100. Another important physical parameter that appears in our set of governing equations is α = μ0mdH0/(kBT). This is a classical physical parameter used in the ferrohydrodynamics context and it expresses a relation between the intensity of forces due to an external magnetic field of intensity H0 exerted on the particles with dipole moments of magnitude md and Brownian forces. Again, T is the absolute temperature of the fluid and kB is the Boltzmann constant. λ = μ 0m2d/(4πkBTd3) denotes the ratio between the induced force by magnetic dipole–dipole interaction and Brownian forces, with d = 2a. Another important parameter that does not appear explicitly in the governing equations for the motion of the particles is the volume fraction of the particles ϕ = nvp, where vp is the volume of a particle and n is the density number, n = N/V, with V being the total volume of the suspension. 2.1.4. Computation of periodic magnetic particle–particle interactions induced by torque In this work we also compute particle–particle interactions induced by magnetic torques with the use of periodic boundary conditions by replicating a finite number of particles periodically throughout the suspension volume. The sum of particle interactions via torque becomes a so-called lattice sum and Ewald summation technique in the physical and reciprocal spaces can be used to accelerate convergence of the resulting lattice sums [51]. As the magnetic particle interactions for the force has a faster decay like (1/r4) the sums of magnetic interactions via induced forces will not be calculated by lattice sums here. That is the mean reason why we have called our numerical method hybrid. On the other hand, works developed by [52–54] and [55] suggest that the slower decay of magnetic torque interactions, (1/r3), can lead to a divergence of the suspension transport properties that depend strictly on rotational dynamics. One property that fits into this description is the suspension magnetization. In order to compare the range of validity in which we may use a hybrid approach for computing particle–particle interactions by induced torque, fluid magnetization are computed and compared for both cases: simulations with periodic and also nonperiodic particle interactions in the torque sums. The dimensionless magnetic torque exerted on a particle i counting the interactions with the neighboring particles j and particle images are computed by the lattice sums technique in the physical, x, and wavelength, k, spaces as follows [51] 2

ð23Þ ð24Þ



TB ¼

Tm

8λ X 1 ¼ 4− T ðxÞ þ 3 Per x∈L 1 L

X ^ k∈ L;k≠0

3 T 2 ðkÞ5 þ


 ^ α di  h Per

ð27Þ

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

where T 1 ðxÞ and T 2 ðkÞ are given by




 T 1 ðxÞ ¼ di  d j B ri j − di  ^ri j d j  ^r i j C r i j

ð28Þ


2


 πk ^ d  ^k e ξ cos 2π ^k  ^r : T 2 ðkÞ ¼ 4π di  k j ij

ð29Þ

  Here B(rij) and C r i j are scalar functions given by 
 
 2ξr  2 2 −ξ r i j Þ −3 ij ri j B r i j ¼ erf c ξr i j þ pffiffiffi eð π

ð30Þ


 C ri j ¼

number of boxes in the physical and reciprocal spaces. For a given number of boxes (in this work we used 125 lattices in our periodic simulations), this condition is controlled by the parameter ξ. For ξ = π1/2V−1/3, the method is both accurate and computationally efficient. The computations of the magnetic torque dipole–dipole interactions between all particles in the physical box and the image particles on the lattice are accounted in each time-step. When a particle in the physical box reaches its bottom, it returns to the top. The same occurs for particles that cross any of the faces of the central lattice. A mathematical definition of the periodicity condition is given by P i j ¼ ð1 þ φi Þδi j

and
 !  3erf c ξr i j 2ξr i j
2 2 −ξ2 r 2i j Þ ð pffiffiffi 3 þ 2ξ ri j e ; þ r 4i j r 4i j π

Z

x

−t 2

e

dt;

0

ð33Þ

where ð31Þ

where erfc is the complementary error function, defined as 2 erfcðxÞ ¼ 1− pffiffiffi π

151

ð32Þ

and ξ is a parameter, responsible for controlling the convergence of the method. In this paper we used ξ = π1/2V−1/3, where V is the volume of the central lattice [51]. 3. Numerical method and preliminary results We have developed a non-commercial Fortran based research code. This code was fully programmed by the authors and is responsible for accounting for each individual particle's translational and rotational motion. It uses a non-commercial Intel Fortran Compiler to produce an executable file that reads an entrance text file with the characteristics of the suspension that will be simulated. It also makes all the required statistics of the suspension and provides the transport properties that we will present in the next sections based on an ensemble average through several simultaneous numerical experiments. For the postprocessing we use TecPlot360, which enables us to produce geometrical representations of the suspension's microstructure in different instants of time by reading the output files generated through our numerical code. In this section we present in details our numerical method for constructing dipole–dipole interactions among N particles suspended in a volume V under condition of an applied magnetic field. The numerical method use a fourth order Runge Kutta scheme to solve the time developing of the suspension dynamics governing by the dimensionless differential Eqs. (18) and (19). As mentioned before, we use a hybrid method for accounting magnetic interactions. Particle interactions induced by force-dipole (which decays like 1/r4) are computed without periodic sums, i.e. particle images throughout the periodic space of the replicated suspension do not entry in the interaction sums. However, depending on the particle volume fraction ϕ and the parameter λ the interactions induced by torque-dipole are computed by the same non-periodic approach used for force-dipole interactions or by using a periodic boundary conditions approach at higher particle volume fraction, beyond 5%. Consequently, the latter case requires the use of Ewald summation technique to accelerate the convergence of resulting sums [41,51]. For simulations using non-periodic interactions for both force and torque, the presence of physical walls on the simulated box is considered. On the other hand, under condition of periodic torquedipole interactions a number of boxes surrounding the main cell defines a lattice system where the number of particle in the physical box is replicated periodically in the wholly space. These periodic boxes are incorporated for accounting the image particles contribution to the particle interactions sums. The convergence of the method depends on the

8 < −li ; φi ¼ 0; : li ;

if xi N li if 0 b xi b li ; i ¼ 1; 2; 3 : if xi b 0

ð34Þ

Here Pij(l1, l2, l3) is a second rank tensor operator applied to the vector xi, that represents the position of each i particle of the suspension in a given configuration. In this notation δij is the Kronecker's delta operator and l1, l2 and l3 denote the lengths of each edge of the prismatic lattice. The periodic simulations here have used 125 boxes in order to replicate a finite number of particle (e.g. N = 1000) periodically throughout the suspension volume. This procedure enhances the computational time substantially. An analysis of the computational costs when using such a procedure is also presented this work. Details on the code structure, convergence, numerical stability analysis and validation are given in the next section. 3.1. Generation of the initial particle configuration An important step required for the time developing process of the suspension is the generation of an initial particle distribution with random and independent positions of the particles inside the simulated box. This initial condition is performed at each realization of the dynamical simulations. Fig. 4 depicts a typical initial configuration of the physical box. If any particle overlap occurs during the generation of the initial condition a random displacement with a Brownian routine is imposed on the collided particles in order to generate initial conditions as close as possible to a spatially statistical homogeneous particle distribution in the absence of particle clusters or any heterogeneity. The initial configurations of the particle dipole moments are also generated with a procedure of random and independent orientations. 3.2. Definition of a numerical time-step We define the time-step of our dynamic numerical simulations in terms of the suspension's physical parameters. The time step is taken as the following relation

 St Pe ; ; 0:001 : δt ¼ min 10:0 10:0

ð35Þ

Here min represents minimum value. Since Stokes and Péclet numbers represent ratios of time scales related to the dynamics of the particle motion, they are used in order to define a numerical time-step for the simulations. The value 0.001 works as a numerical filter to avoid high time-steps that may cause particle overlap when two particles are sufficiently near. It should be important to note that a static timestep is being used in this formulation. Actually, when two particles have a tendency of overlapping by the influence of magnetic attractive

152

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165 Z

Z

X

X

Y

(a)

Y

(b) Z

X

(c)

Y

(d)

Fig. 4. Several initial configurations for 1000 particles with volume fraction ϕ = 10 %. (a) A three-dimensional view of an arbitrary realization. (b) A side view of the same initial configuration. (c) Details of particle microstructure. (d) Side view without the initial dipole moments of the magnetic particles.

forces, an interesting, but undesired numerical phenomenon occurs. Expression (21) shows that the magnetic force due to dipole–dipole interaction scales with (1/r4), where r is the distance between the centers of the particles. If two particles approach each other by attractive magnetic forces, the intensity of the numerical value of this attractive force increases very fast. In order to avoid particle overlap it is important to turn off this interaction for very close particles. For example, for St = 0.1 and λ = 0.5, the distance between the surface of the particles in which we must turn off magnetic attractive forces to avoid particle overlap is about 0.1a. Since the particle is massive there will be a time-delay between the instant in which this force is turned off and the instant that the particle responds to the action of its neighboring.

Particle 1 a

Even using a fictitious repulsive force it is difficult to suppress the effect of attractive magnetic forces when the particles touch each other in the presence of particle inertia. The numerical gap between the particles was calibrated based on the minimum distance in which the magnetic attractive forces could be suppressed to prevent particle overlap. We consider a simple experiment of two neutrally buoyant particles with initial positions fixed in a distance of ten particle radius from each other. Gravity forces and external magnetic fields were turned off, being the particle motion induced only by attractive magnetic forces between particle dipoles (see Fig. 5 for details of particle configuration). Several values of Stokes number St and the magnetic parameter defined as α* = 24λ were tested for constants C1 = 10 and C12 =

Particle 2 a

Di = Di 1

Di = Di 2

10 a

y

x Fig. 5. Particle configuration used to calibrate the numerical gap in order to avoid particle overlap induced by magnetic interaction. The dipole moments of the spheres have the same magnitude and sense in order to induce particle attraction. The dipole orientation may be interpreted as a positive monopole at the end of the arrow while the begging is a negative monopole. Poles with opposite signs attract each other.

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

0.01 of the repulsive force giving in Eq. (22). The values of these constants were settled up in order to produce a smooth motion during particles approaching each other in very close distances. Under these conditions, the behavior of the numerical gap as a function of the magnetic interacting parameter is shown in Fig. 6. With the gap characterization for different Stokes number and dipole–dipole magnetic interacting parameter, particle overlapping is avoided during the simulation even under condition of massive particles. Based on the results presented in Fig. 6, we have proposed the following fitting curves to be used in the numerical simulations in order to control particle overlapping at several conditions of St and α *:   0:0975 : For St ¼ 0:1→ Gap α ¼ 1:6853α

ð36Þ

  0:1399 : For St ¼ 1→ Gap α ¼ 1:8401α

ð37Þ

  0:1628 : For St ¼ 10 → Gap α ¼ 2:3409α

ð38Þ

With this gap calibration it is possible to run dynamic numerical simulations for magnetic suspensions with particle inertia with a fixed time-step. This implies a lower computational cost in our simulations, since the use of prohibitive small time steps is avoided for approaching particles. In this work, we will consider at least a small inertial effect (i.e. St ~0.1) in order to compute and explore some new physics from the dynamic behavior of a magnetic suspension of slightly massive particles. 3.3. Numerical evolution scheme The forces acting on all the particles at different realizations are calculated simultaneously. The velocity of each particle is computed by using a fourth order Runge–Kutta scheme and the position of the particles at the subsequent time-step is also calculated by numerical integration of dxi i ¼u; dt

ð39Þ

where the upper index i denote the ith particle in a particular realization. The same numerical procedure is also used for computing the 8

di

nþ1

 n n n ¼ di þ ω  di δt;

ð40Þ

where the subscript i denotes the ith particle and the superscript n represents the current step of the simulation. 3.4. Stability and convergence of the method It is well-known [37,39,40] that the calculations of particle hydrodynamic interactions results in a non-convergent problem on the statistical properties of the suspension due to slow decay of these interactions like (1/r), where r represents the distance between the centers of the particles. In order to avoid such a problem a periodic system with periodic boundary conditions should be used. The disadvantage in using periodic boundary conditions is related to a significant increase of the computational cost in dynamic simulations. In the context of magnetic interactions, however, dipole–dipole interactions decay considerable faster then hydrodynamic interactions, like 1/r4 for the force dipole interactions and 1/r3 for torque dipole interactions. We first performed dynamic simulations using non-periodic magnetic many-body interactions for both force and torque dipoles. Next a hybrid method is proposed with periodicity for particle interaction in a periodic box only for the torque-dipole interactions. We shall see that this hybrid procedure demands much lower computational cost when compared with simulations based on the calculations of full periodic many-body interactions for torque and force dipoles. In the present section it is also shown that the proposed hybrid approach is both accurate and computationally efficient, and the statistical properties of the suspension such as the velocity variance provide convergent values. 3.4.1. Statistical quantities The particle average velocity of the suspension at an arbitrary timestep and a specific realization is defined as follows uð j; t Þ ¼

N 1X u ð j; t Þ: N i¼1 i

ð41Þ

Here j denotes the realization and the subscript i represents a given particle of the suspension. Now, the square of particle velocity fluctuations in a typical j realization at the same time-step is calculated as being 02

2

ð42Þ

We assume a statistical process with 〈u′〉 = 0. In this case, the mean velocity variance σ ′ = 〈u′2〉 − 〈u′〉2 reduces to 〈u′2〉. So, the mean velocity variance of the magnetic suspension is given by the following expression

6

Gap

angular velocity of each particle and to evolve the dipole moment of the particles with the following discrete Eq. (39).

u i ð j; t Þ ¼ ½ui ð j; t Þ−uð j; t Þ :

St = 0.1 St = 1.0 St = 10.0

7

153

D

5

N rea X N E 1 1X 02 02 u ð j; t Þ; u ðt Þ ¼ Nrea N j¼1 i¼1 i

ð43Þ

4

where Nrea represents the number of realizations. Note that 〈u′2(t)〉 is a vector quantity that measures velocity fluctuations induced by particle– particle magnetic interactions. For account those fluctuations we define the square root of the mean velocity variance as

3

2 200

400

α*

600

800

1000

Fig. 6. Numerical gaps in which the magnetic attractive force must be suppressed in order to prevent particle overlap when particle inertia is considered. The thin continuous line denotes St = 10, the dashed line denotes St = 1 and the thicker continuous line St = 0.1.

ζ ðt Þ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 2 0 2 0 2 〈 ux ðtÞ 〉 2 þ 〈 uy ðtÞ 〉 2 þ 〈 uz ðtÞ 〉 2 :

ð44Þ

In some cases we will analyze the mean velocity variance for a fixed time in which we have no aggregate formation. In these cases we will refer to this quantity as the equilibrium mean velocity variance ζ0.

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

2

1.5

0.1

ζ0

3.4.2. Active range of the suspension Since this formulation considers the presence of physical walls, we have defined an active range of the suspension in which the statistical treatment of the suspension is always representative. Particles that are close to the wall have their velocities strongly influenced by near solid boundaries and are consequently excluded from the statistical treatment of the suspension. We define a suitable distance D from the walls in order to always have a numerical box with a sufficient number of active particles where the statistical analysis of the data should be representative (Fig. 7). In the absence of cluster formation a numerical experiment has been done for computing the mean value of the velocity variance as a function of D for the following combinations of the physical parameters Pe = 1000, St = 0.1, α⁎ = 10, α = 10 and ϕ = 10.0%. These simulations were performed with 500 particles and 3 simultaneous realizations. The results of the velocity variance with D ranging from zero to five particle radius can be seen in Fig. 8. The numerical test indicates that for a distance D from the wall equal to 4 particle radius the mean value of the equilibrium velocity variance is already convergent. By this result we define an active range of the suspension as being the internal box with dimensions: (‘ − 4a) × (‘ − 4a) × (h − 4a), where ‘ denotes the length and h is the height of the numerical. The typical behavior of the square root of the mean velocity variance of a non-Brownian magnetic suspension varying with time is presented in Fig. 9. The statistics was calculated over ten independent realizations with 300 particles. We examine the behavior of this statistical quantity for the particle velocity fluctuations being produced by dipole–dipole magnetic interactions only. We can see in Fig. 9 the presence of a distribution of discrete impulses of the mean velocity variance of the suspension. However, between the impulses the variance assumes a constant value. Usually, in a non-magnetic suspension with particle–particle hydrodynamic interactions the mean velocity variance is a statistically stationary quantity [37,39]. In contrast, this behavior has not been observed in the current simulation of a magnetic suspension in the absence of hydrodynamic interactions. We argue that the presence of particle clusters in the suspension induced by dipole–dipole magnetic interactions leads to this time dependence behavior of the variance

0.05

ζ0

154

1 0 4

D

4.5

5

0.5

0 0

1

2

3

4

5

D Fig. 8. Particle velocity variance as a function of the wall distance D. The points represent numerical values and the straight line denotes the average of the last three points.

seen in Fig. 9. Fig. 10a shows for a very short time interval the details of the mean velocity variance behavior. In addition, Fig. 10b gives the distribution at the same time interval of the number of clusters with exactly three particles, denoted by N3c . In these findings we can see a direct relation between the formations of typical 03 particle-clusters induced by dipole–dipole interactions as the time evolves with the discrete distribution of the variance impulses related to the particle velocity fluctuations. This variance time dependence occurs because, during the formation of a cluster some very close particles may have velocity fluctuations much greater than the system average velocity as a consequence of the strong dipole–dipole interactions which decays like 1/r4. In other words, in this near field configuration, a particle can experience a high velocity fluctuation in a very short time, that justifies the discrete impulses distribution shown in Fig. 9. Now, when considering just the values of the mean velocity variance that are not influenced by the variance time dependence during a cluster formation, we can present a study of the variance dependence on the system size. Fig. 11 shows a plot of the mean variance ζ0 in equilibrium

0.12

ζ(t)

0.08

0.04

0 5

10

15

20

25

t Fig. 7. A numerical box defined as being an active range of the suspension in order to compute the average macroscopic quantities of the suspension with a sufficient number of mobile particles and representative statistics.

Fig. 9. Typical mean velocity variance behavior in the dynamic simulation of this work for a magnetic suspension with St = 0.1, Pe = 1000, α* = 20 and ϕ = 5.0 %.

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

(a)

(b)

0.1

30

0.08

28

0.06 3

26

Nc

ζ(t)

155

0.04

24 0.02 22 0 20 10.7

10.8

10.9

11

10.7

10.8

t

t

Fig. 10. (a) Dimensionless mean velocity variance as a function of time for a non-Brownian magnetic suspension; (b) number of clusters with three particles as the time evolves for the same interval of the plot (a). In these plots St = 0.1, Pe = 1000, α* = 20 and ϕ = 5.0 %.

(taken in a given time without aggregate formation) as a function of the number of particles for a fixed particle volume fraction, ϕ. It is seen that for a number of particles N = 333 the mean particle velocity variance converges to a constant value equal to 6.1 × 10−3 for ϕ = 0.1 % and 6.3 × 10−3 for ϕ = 1.0 % with Pe = 1000 and λ = 0.1. This value remains constant as the number of particle in the numerical box increases. This finding indicates that the faster decay of the dipole–dipole interactions (i.e. 1/r3, 1/r4) was sufficient to produce a variance associated with velocity fluctuations independent of the system size. This is in contrast with the problem of variance divergence in the context of particle hydrodynamic interactions of sedimenting particles, where the interactions decay slowly, like 1/r. Other important property to be examined with respect to numerical convergence is the equilibrium magnetization of the magnetic suspension. This macroscopic physical property is defined as the average orientation of the particles dipole moment with respect to the direction of an applied external magnetic field. The equilibrium magnetization M0 in a given time-step can be calculated as M0 1 1 j ¼ hcos θii ; M s N Nrea

ð45Þ

where Ms is the saturation magnetization, θ is the angle between the direction of the dipole moment of a particle i in a realization j at a given time-step and an applied magnetic field direction. The results of the numerical simulations for the equilibrium magnetization shown in Fig. 12 with a particle volume fraction of ϕ = 0.05

(a)

reveal again a convergence of the suspension magnetization as the system size increases. We can see that the convergence is reached for a number of particles N = 500 and the value of the equilibrium magnetization M0 = 0.9Ms. 4. Results and discussions 4.1. Magnetization One of the most important properties of a magnetic suspension is its equilibrium magnetization. This macroscopic property is closely related to the particle's rotational motion under the action of an applied magnetic field. The equilibrium magnetization is defined as the average alignment of the particle's magnetic dipoles in the direction of the applied field. We can see this rotational mechanism from the numerical simulations developed here. Fig. 13 shows a typical response of the dipole moments orientations in the presence of an external field. Fig. 13a corresponds to the initial condition (i.e. t = 0) in the absence of a magnetic field. Each particle has a dipole moment aligned independently in a random direction in the numerical box. In this configuration the average orientation of dipole moments over all particles is zero, which implies a null magnetization in the considered fluid volume. Fig. 13b reveals how the particles dipoles respond to an applied magnetic field. This configuration was taken from the dynamic simulation in t = 0.2 Stokes time. We can expect that with the application of an external field a magnetic torque, Tm = mp × H, is exerted on each particle due to the interaction between the external field vector and the dipole moment of the individual particle. In all of our simulations a

(b) 0.008

ζ0

ζ0

0.008

0.004

0.004

10-3

10-2

1/N

10-1

10-3

10-2

10-1

1/N

Fig. 11. Variance convergence with respect to the system size (a) ϕ = 0.1 % and (b) ϕ = 1.0 %. The circles represent the numerical values, while the straight line denotes the average of the last three points (for the smallest values of 1/N).

156

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

1

4.2. Theoretical models and code validation The simplest theoretical model used to calculate the equilibrium magnetization of a very dilute magnetic suspension of non-interacting dipoles is the well-known Langevin's model [11]. This is therefore an Oϕ model, given by Eq. (45).

M0 / Ms

0.95

 1 ¼ Md Lðα Þϕ; M0 ¼ Ms coth α− α

0.9

0.85

0.8

10-3

10-2

10-1

1/N Fig. 12. Convergence of the equilibrium magnetization with respect to the system size for ϕ = 0.05. The circles represent the numerical values, while the dashed line denotes the average of the last three points (for the smallest values of 1/N).

typical scale of the magnetic relaxation time and the particle rotation is much smaller than the Neel time scale so that the case of fixed dipole moments on the particles has been considered. In the presence of an external field the magnetic torque produces particle rotational motion and consequently the particle's fixed dipole moments tend to align in the direction of the applied field. In a condition in which the ratio of magnetic and Brownian torques α is sufficiently large (i.e. α ~ 1) the net alignment of the particle dipole moments is seen to be in the field direction as revealed in Fig. 13b. When all particles dipole moments are aligned in the field direction we say that the equilibrium magnetization of the suspension reaches its saturation value Ms. In this condition the magnetization of the suspension is simply the volume average ϕMd, where Md denotes the magnetization of the solid material of the particles.

ð46Þ

where L is the well known Langevin function [11]. In addition, there are also more robust theoretical models considering O(ϕ2) and O(ϕ3) corrections of dipole–dipole interactions for two and three particles. In these models the equilibrium magnetization can be predicted as a function of α, λ and ϕ [28,42]. The accuracy of the numerical simulations developed here is verified by comparing the results of our numerical experiments with these models. Fig. 14 shows a plot of the equilibrium magnetization as a function of the parameter α for a volume fraction ϕ = 0.05, a Péclet number of Pe = 0.1 and λ = 1.0. We can see a close agreement between the first (Langevin [11]) and second order [42] models with the present numerical values. For conditions of stronger influence of the magnetic field. i.e. α ≫ 1, all models collapse within the numerical simulations. This occurs because in the presence of more intense magnetic fields Brownian thermal fluctuations are unimportant and the influence of dipole–dipole interactions become a small effect compared with the effects of an external field acting on the particles. So, the parameter ϕ and λ produce a weak interaction effect for λ ≫ 1. The biggest differences between the models shown in Fig. 14 focuses on regions where α ~ 1. When α ~ 1 the energy of the external field is compared with those from Brownian fluctuations. In this limit particle's interactions are not overwhelmed by the external field contribution. For α ≫ 1 the external field dominates particle interactions and all the models tend to the same value. We can also see that the numerical values tend to underestimate the higher order asymptotic solutions. This is an indication that non-periodic dipole–dipole magnetic interactions computations are valid only for small or moderate values of ϕ and λ.

Z Z

X X

(a)

Y

(b) Fig. 13. Typical time evolution of the equilibrium magnetization of the suspension — ϕ = 0.05.

Y

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

157

0.4 1

0.6

M0 / Ms

1

0.8

M0 / Ms

M0 / Ms

0.8

0.4

0.6

0.8

0.2

M0 / Ms

0.4

0.2

0

0

0

1

2

3

4

α

5

6

0

0 0

5

10

15

α

0

0.4

0.8

λ

Fig. 14. Equilibrium magnetization as a function of α for ϕ = 0.05 and λ = 1.0. The simulations were performed for 500 particles in the numerical box and the statistics was taken over 50 independent realizations of the suspension. The solid line denotes the Langevin
 model, the dashed line represents the O ϕ3 model of [42], the dash–dot line denotes Jansons' model [28] and the filled black circles denote the present numerical results with non-periodic torque computations. We also show in the insert a plot of the magnetization at lower values of α. As α → 0 we can identify a paramagnetic regime with the linear dependence being fitted by the known asymptotic solution for α ≪ 1: M0/ϕMd = α/3.

Since the biggest discrepancy between the theoretical models and numerical simulations is seen for α ~ 1, we present a plot in Fig. 15 that shows the behavior of the suspension magnetization as a function of the volume fraction of particles for α = 1. Now we shall present results that show that the numerical procedure developed here is both accurate and computationally efficient. We shall see that, in certain limits the numerical results are not just accurate, but also represent practically an exact solution.

0.06

M0/Md

1

λ

0.2

20

0.04

Fig. 16. Equilibrium magnetization as a function of λ for ϕ = 0.05, α = 1.0. The simulations were performed for 500 particles with statistics taken over 50 suspension realizations. The
 solid line denotes the Langevin model, the dashed line represents the O ϕ3 model of [42], the filled black circles denote the numerical results for a non-periodic computation of magnetic dipole–dipole torque interactions, while the blank circles represent numerical results for periodic computations (i.e. using lattice sums for particle interactions) of magnetic dipole–dipole interactions for torque. We also show in detail the behavior obtained with error bars in a larger range of the y-axis.

For α ~ 1 the intensity of the forces and torques induced by the applied field is in the same order of magnitude as translational and rotational fluctuations induced by Brownian motion. In this regimes particle motion is not dominated by dipole interaction with the applied field. So, depending on the value of λ dipole–dipole magnetic interaction become also important, being possible to investigate by our computer simulations how the steady state equilibrium magnetization depends on particle volume fraction. Fig. 15 shows that for small values of ϕ ≪ 1 (i.e. around 0.01) all theoretical models and the numerical solutions give very similar results. In this dilute regime, for λ ~ 1, dipole–dipole magnetic interactions are not comparable with other forces and torques that govern the suspension dynamics. However, when increasing particle volume fraction, the dipole–dipole magnetic interactions become important on the suspension evolution. In this regime, it is seen that the numerical predictions computing particle interactions without periodic sums (i.e. using lattice
 sums with particle images) tend to underestimate the O ϕ3 asymptotic theories. So, these results indicate that dipole–dipole magnetic interactions when computed in a non-periodic way tend to underestimate the equilibrium magnetization values for more concentrated regimes (i.e. ϕ ~ 0.05 for λ ~ 1). Even in dilute regimes we already observe small anisotropic chains of two or three particle dipoles that are close enough to interact magnetically. The numerical code developed here is

0.02

0 0

0.05

φ

0.1

0.15

Fig. 15. Dimensionless equilibrium magnetization as a function of ϕ for α = 1.0 and λ = 1.0. The simulations were performed for 500 particles inside the periodic box with statistics taken over 50 different realizations. The solid line denotes the Langevin's model [11],
 the dashed line represents the O ϕ3 model of [42], the filled black circles denote the numerical results for a non-periodic computation of magnetic dipole–dipole torque interactions, while the blank circles represent numerical results for periodic computations of magnetic dipole–dipole interactions for torque (i.e. using lattice sums for particle interactions).

Table 1 Computation time × number of particles for computation of periodic torques and nonperiodic computations. N

Periodic torques (s)

Non-periodic (s)

50 100 200 500 700 1000 1500

2.46 9.43 38.18 333.42 653.21 1315.30 3017.09

0.04 0.11 0.41 2.93 5.98 13.60 31.55

158

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

Table 2 Relative error × volume fraction of particles for computation of periodic torques and nonperiodic computations. ϕ

Periodic torques (%)

Non-periodic (%)

0.01 0.03 0.05 0.07 0.09 0.11

1.54 2.03 1.08 0.74 1.45 2.42

3.10 7.04 9.53 12.94 14.23 15.68

able of capturing structure transition and chain formation. We shall see this dynamics in the next sections. Our numerical procedure performs the magnetic interactions between all the particles within the suspension and, it is not limited by the boundaries imposed in the theories. We may then speculate, that the use of non-periodic (i.e. minimum image procedure) dipole– dipole magnetic interactions to represent an infinite suspension may lead to a convergence problem in the transport properties calculations of the suspension in more concentrated regimes. On the other hand, the use of periodic computations using lattice sums in the reciprocal and physical spaces for magnetic torques due to dipole–dipole interactions has produced results in excellent agreement with the theoretical O
 ϕ3 model [42] for values of ϕ up to 0.1 and λ ~ 1. For still higher values of ϕ we can see an appreciable deviation between the numerical simulations results and the higher order asymptotic theory. This discrepancy can be explained by the account of all dipole–dipole magnetic interactions by the numerical method while just three interacting particles


 are considered by the theory O ϕ3 . So, the numerical simulations are quite important to simulate more concentrated suspension. In addition, Fig. 16 shows a plot of the dimensionless equilibrium magnetization as a function of the dipole–dipole interaction parameter, λ. It is seen that numerical simulations just considering non-periodic magnetic particle interactions produce an average value of the magnetization pretty close to the Langevin model applied for very diluted magnetic suspensions. However, when we compute the magnetic particle interactions in a periodic way by using lattice sums for the magnetic torque dipoles (not for magnetic forces which decay like 1/r4) the results are seen to be in a very good agreement with the asymptotic
 solution O ϕ3 [42]. The insert in Fig. 16 shows the error bars related to the statistics of the numerical simulation over 50 realizations of 500 particles. We can see that both simulations results using periodic (i.e. lattice sums or multiple images) and non-periodic (minimum image) magnetic dipole–dipole particle interactions are inside the error bars. Consequently, at this low particle volume fraction, ϕ = 0.05, and λ ~ 1 even the computations without periodic interactions produce acceptable results inside the error bars interval. This scenario may change significantly as ϕ and λ increase since these parameters increase the effect of dipole–dipole magnetic interactions on the particle motion. This is stronger for dipole–dipole torque interaction that decays like 1/r3 instead of dipole–dipole force interaction which goes like 1/r4. For dipole–dipole force interactions we use the procedure of minimum image (non-periodic procedure) for all simulations. We have verified through several tests that computing just dipole–dipole torque interactions periodically (i.e. with lattice sums) is sufficient to predict correct values of the equilibrium magnetization even at moderate particle volume fraction and λ, i.e. ϕ ≈ 0.1 and λ ~ 1.

Z

X

Z

Y

(a)

X

(b) Z

X

(c)

Y

Z

Y

X

Y

(d)

Fig. 17. Configurations of the particle trajectories in a magnetic suspension for different values of the physical parameters. (a) Trajectory of a test particle in a magnetic suspension in the absence of magnetic and hydrodynamic particle interactions, φm = 0, and Brownian motion, Pe → ∞. (b) This particle trajectory considers magnetic interactions only for Pe → ∞ and φm = 1. (c) Represents a particle trajectory for a suspension without dipole–dipole magnetic interactions, i.e. λ = 0, but undergoing Brownian motion with Pe = 1. (d) Particle trajectory in the presence of dipole–dipole magnetic interactions and Brownian motion with λ = 1 and Pe = 1.

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

It should be important to note that the results presented in Fig. 16 were obtained for α ~ 1. In this regime the particle suffers random fluctuations due to Brownian motion which produces an intrinsic tendency of numerical results with large error bars. When we show the results without the error bars as depicted in the plot of Fig. 16, the difference between the numerical results using the procedure of non-periodic and periodic computations of magnetic torque interaction increase as λ increases. The same behavior is seen when comparing the numerical
 results with the asymptotic solution O ϕ3 for both procedures of computing dipole–dipole magnetic interactions. In addition to these results, we present in Tables 1 and 2 a summary of our numerical tests, showing the computational cost and the relative error when using the numerical procedure of minimum image for computing particle interactions (i.e. non-periodic procedure) and when performing the numerical simulation with lattice sums or multiple particle-images as a mean of representing an infinite suspension. We
 use the asymptotic solution O ϕ3 developed by [42] as the exact reference to quantify the relative errors. The simulations were performed in a serial processing scheme with an i7 processor of 2.00 GHz and 8 Gb of RAM. The numerical research code, developed by the authors uses Linux Platform and runs in a serial processing scheme, but performing different numerical realizations simultaneously. The great increase in the computational cost due to periodic computations comes from the fact that in this case we must account for particle interactions not only in the physical lattice, but also in 125 imaginary lattices spread in the physical and reciprocal space. From the analysis done so far we can conclude that the minimum image procedure for computing dipole–dipole magnetic interactions in a non-periodic way can be applied under two conditions. Firstly, for a dilute suspension (i.e. ϕ ~ 0.01) we could use non-periodic computations

159

of magnetic particle interactions to estimate the ferrofluid's magnetization. Even for arbitrary values of λ and α the predictions are accurate and the method seems to be computationally efficient. As we can see in Tables 1 and 2, a considerable amount of computational time can be saved when using this non-periodic procedure for particle interactions. Secondly, for more concentrated suspensions, i.e. ϕ N 0.01, this methodology would produce accurate results only for λ ≪ 1 at arbitrary values of α, or for the extreme case when a strong external field is applied, α ≫ 1 with an arbitrary values of λ . Table 2 indicates that in non-dilute regimes of the suspension with α ~ 1 and λ ≥ 1 we need to use the hybrid method with periodic dipole–dipole interactions for torque dipoles and non-periodic interactions for force dipoles. This guarantees relative errors bellow 3% for ϕ around 0.1. In the next sections we will show the potential of the developed code for exploring structure transition, typical configurations of Brownian and non-Brownian magnetic suspensions and chain formation in the presence and absence of an external applied magnetic field. 4.3. Typical configurations and structures In this section we examine typical configurations of the suspension's microstructure for Brownian and non-Brownian magnetic suspensions. When exploring Brownian regimes α and λ are used as the magnetic parameters. For non-Brownian suspensions we introduce two new dimensionless parameters ψm and φm. These parameters are defined respectively as follows:

φm ¼

μ 0 m2d ; 8π2 ηa51 U s

and ψm ¼

md H 0 : 6πηa21 U s

ð47Þ

Fig. 18. Time evolution a Brownian magnetic suspension during sedimentation. In these simulations ϕ = 0.05, Pe → ∞, ψm = 20, φm = 1.6 and St = 0.1. (a) Initial particle configuration in the suspension with the dipole moments of each particle distributed randomly and independently; (b) denotes an amplification of the suspension core after 33 a/Us; (c) suspension microstructure at t = 66a/Us; and (d) final particle distribution at t = 100a/Us.

160

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

The choice of the parameters ψm and φm is more appropriate for the context of non-Brownian suspensions since Pe cannot be made infinity and it still accounts for magnetic effects by using Eqs. (18) and (20). So when dealing with simulations of non-Brownian suspensions the parameter φm accounts for the contributions of force/torque dipole–dipole interactions whereas ψm for the intensity of the applied field. The effect of dipole–dipole magnetic interactions can be seen by tracking the trajectory of a test particle obtained for a condition of St = 0.1 and ϕ = 0.05. Fig. 17a, b shows the action of the dipole–dipole magnetic interactions on the trajectory of a test particle in the suspension. For φm = 1 as shown in Fig. 17b particle trajectory presents fluctuations compared with the vertical straight-line (i.e. for φm = 0) of Fig. 17a as a direct consequence of the particle magnetic interactions. It would be important to note that the trajectory developed by a test particle in a non-Brownian suspension has a more deterministic characteristic than the motion of a test particle in a Brownian suspension even in the presence of particle inertia (Fig. 17c). The presence of particle inertia implies a nonzero particle relaxation time and the effect of the Brownian random walk is attenuated on the particle trajectory since the particle does not respond instantly to the random kicks imposed by the carrier fluid molecules. Finally, Fig. 17d shows a result due to

the combined effect of Brownian motion and magnetic interactions on the motion of a test particle. It is observed in this condition the addition of the random character on the particle trajectory by the presence of the Brownian motion. So, under this condition the particle motion presents trajectory deviation due to both dipole–dipole magnetic interaction and Brownian motion. The particle trajectories illustrated in Fig. 17d present a scenario where the effect of magnetic interactions for λ ~ 1 is in the same order of Brownian motion fluctuations when Pe ~ 1. In this sense particle fluctuations induced by magnetic interactions and Brownian motion can produce an important mixing in a magnetic suspension which is different of the standard ordinary thermal diffusion due to pure random velocity fluctuations. Fig. 18 shows a typical time evolution of a nonBrownian magnetic suspension during sedimentation. This numerical simulation was performed for 1000 particles, St = 0.1, ψm = 20, φm = 1.6, ϕ = 5.0 % and Pe → ∞. It is seen that in a short time interval, small clusters of few particles (i.e. two or three particles) are formed and as the time evolves long chains aligned with the field direction can be observed. It is important to observe that the initial condition is generated in the absence of any particle cluster. Therefore, the observed structures are all formed during

Fig. 19. Typical time evolution process of a magnetic suspension. For images a and b we have ϕ = 0.05, Pe → ∞, ψm = 0, φm = 0.1 and St = 0.1 (no applied field). For images c and d: ϕ = 0.05, Pe → ∞, ψm = 10, φm = 0.1 and St = 0.1 (applied field). Configurations a and c correspond to the initial condition with the particles randomly and independently distributed. Pictures b and d show the suspension configuration after 100 a/Us of simulation.

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

the dynamic simulation as a direct consequence of the fluctuations produced in particle trajectories by force–torque dipole–dipole magnetic interactions and dipole–magnetic field interaction. Actually, a rapid orientation of the particle dipoles in the direction of the applied external field is closely related with torque dipole–field interactions. In addition, we show in Fig. 19 the differences in the cluster structure induced by particle magnetic interactions. These pictures of the suspension time developing were obtained after two test cases. In Fig. 19a and b we consider only the magnetic effect of dipole–dipole interactions with no influence of an external field. On the other hand in Fig. 19c and d an external field is incorporated to the simulation. The question is to verify the difference between the structure formations in terms of cluster anisotropy in both cases. Again, cluster formation is seen even in a non-Brownian magnetic regimes induced by fluctuations in the motion of the particles due to force–torque dipole–dipole magnetic interactions. In the absence of an external magnetic field as considered in the simulation results of Fig. 19a and b the shape of the clusters formed resembles a more isotropic structure without a preferential direction. In contrast, in the presence of an applied field as shown in Fig. 19c and d the suspension evolves to a more heterogeneous anisotropic structure. In presence of a magnetic field, we observe the formation of small clusters of two or three particles that evolve to longer anisotropic chains aligned to the field direction. This behavior can be related to two physical mechanisms. First the magnetic torque acting on each particle due to the presence of an external field tends to rotate the dipoles in its direction and then, attractive forces due to dipole–dipole interactions start to induce the formation of long chains [56]. Fig. 20 shows a close detail of the microstructures presented in Fig. 19. The amplified image of the suspension shown in Fig. 20b indicates an aggregate structure of eleven particles composed of two dimmers, one micro-aggregate with three particles and another one of four particles. This behavior has motivated us to explore by numerical simulation the topology of an initially magnetic spherical blob sedimenting under action of gravity and with particles inside the blob interacting magnetically. The interesting question to investigate should be on the final structure of a spherical magnetic blob at a sufficiently long time. 4.4. Evolution of a sedimenting spherical blob of magnetic particles Now, we examine how a spherical blob of magnetic particles evolve in time. In Fig. 21 we present a typical time evolution process of a

161

spherical blob of polarized magnetic particles during sedimentation. The hydrodynamic interactions between the polarized particles are not considered in these simulations. Computer simulations investigating the influence of particle hydrodynamic interactions on the blob motion and topology were carried out by [38,57]. In the present simulations we focus on the effects of dipole–dipole interactions on the blob only. In Fig. 21 the volume fraction of particles inside the blob is defined as ϕB, given by

ϕB ¼ N 

a RB

3

;

ð48Þ

where N is the number of particles inside the blob and RB represents the radius of gyration of the equivalent spherical blob. We define RB = 0.5 × max(rij), where max(rij) is the maximum distance between two arbitrary particles inside the blob. The simulations of [38,57] show that individual particles scape from the aggregate varying its topology due to velocity fluctuations produced by hydrodynamic interactions. In contrast, it is seen in Fig. 21 that a magnetic spherical blob tends to preserve its original shape during all sedimentation process as a consequence of attractive forces produced by dipole–dipole interactions. However a slight top and bottom flattening deformation in the blob is observed. An interesting finding resulting from the effect of dipole–dipole interactions is a tendency to form a heterogeneous size distribution of smaller cluster structures inside the blob. So, even the aggregate preserving its topology the particle distribution inside the aggregate changes drastically as the time evolves. Certainly, this internal mechanism could be used to speed up solid–liquid separation process in magnetic suspensions. An interesting question motivated by the results shown in Fig. 21 would be as follows: are dipole–dipole magnetic interactions strong enough to prevent particle scape from the blob even in the presence of fluctuations induced by hydrodynamic interactions? In order to investigate the qualitative behavior of a magnetic spherical blob in sedimentation, we carried out a simple experimental observation in the Laboratory of Complex Flows at the University of Brasilia. For this end, a spherical blob of approximately 1244 iron particles was suspended in a rectangular box of 0.01875 m3 containing pure silicone oil. We considered two situations of the particles inside the blob. First, the blob time evolution process was performed as the

Fig. 20. Comparison of a typical microstructure selected arbitrarily from the middle of a magnetic suspension for ϕ = 0.05, Pe → ∞, φm = 0.1 and St = 0.1. (a) Isotropic cluster structure in a simulation under condition of ψm = 0 (i.e. in the absence of magnetic field) and (b) anisotropic chain structure for ψm = 10.

162

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

Fig. 21. Time evolution of a magnetic spherical blob during sedimentation, for ϕB = 0.03, Pe → ∞, ψm = 0, φm = 0.1 and St = 0.1. (a) Initial state of the blob with particles randomly and independently distributed inside it. b and c correspond to the time evolution scheme of the blob at 50 and 100 Stokes time (a/Us) of simulation, respectively.

16000

980

(a)

(b)

975

(m2)

(m)

970 965 960 955

8000

950 2

4

6

8

2

10

4

Number of samples

6

8

10

Number of Samples 0.09

(d)

f

0.06

0.03

0

(c)

800

1000

1200

d

Fig. 22. Statistical analysis of the particles used in the experiment of a blob during sedimentation. (a) Accumulated average diameter of the particles as a function of the number of samples. (b) Accumulated variance as a function of the number of samples. (c) Typical sample of particles observed in the microscope. (d) Frequency histogram of the particle size distribution in comparison with a Gaussian distribution.

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

particles were not magnetized. However in the second test, the particles inside the blob were magnetized for approximately seven days by using a neodymium permanent magnet that produces an induced magnetic filed of 3000 G. The silicone oil used had a viscosity η = 95 Pa∙s and density ρf = 0.97 g/cm3. The mean diameter of the particles was 972.79 ± 108.01 μm and their mean density ρs = 7.92 g/cm3. The average volume of the blob was vB = 10−6 m3 and its volume fraction about ϕB ≈ 0.5. A rigorous statistical analysis on the particle size distribution was done by using an Olympus BX51 optical microscope with an available data analysis software. The statistic was performed over 10 samples containing approximately 25 particles each. We obtained the particle mean diameter and its associated variance. The results are plotted as a function of the number of samples. Fig. 22a and b indicates that a number of 10 samples is quite enough to provide a convergent statistics with the number of samples. Fig. 22c shows a typical sample observed in the microscope, whereas Fig. 22d depicts the frequency histogram of the particle diameter and a comparison with a Gaussian probability distribution function. Fig. 23 shows a snapshot of the blobs obtained from experimental observations for the cases of non-magnetized particles inside the blob (Fig. 23a) and a magnetic blob (Fig. 23b). In the pictures shown in Fig. 23 we have a blob of non-magnetized particles (a), whereas in b the blob is composed of magnetized particles as described above. In both pictures the blobs have traveled a distance equivalent to 20RB. We can see that in Fig. 23a the blob has a varying topology due to the scape of several non-magnetized particles to the wake region of the blob

163

as have been already observed by [38,57]. This typical wake is formed due to the backflow effect that is closely related to hydrodynamic interactions. These particle interactions have a slow decay like (1/r), so the velocity field induced at the bottom of the box due to the motion of the blob is felt by all the particles within it. A quite different behavior is observed in Fig. 23b. The blob preserves its topology without a systematic particle scape and sedimenting as a rigid body motion. It is seen that the magnetic effect generated by interparticle attractive forces between the magnetized particles are sufficiently strong to suppress velocity fluctuations induced by hydrodynamic interaction across the blob domain. A typical wake does not form when the particles are magnetized, even in the absence of an applied field. This means that the short range nature of magnetic effects (1/r4) dominates the long range of hydrodynamic interactions (1/r) inside the blob domain for strong intensities of the magnetic parameter λ. This justifies why we have obtained through our computer simulations a very similar qualitative picture of a sedimenting blob even in the absence of the hydrodynamic interactions. A clear illustration of this scenario can be seen by a comparison between Figs. 21 and 23b that shows the blob time developing process obtained by numerical computations and experimental observations, respectively. During our experiments to produce Fig. 23 we have seen that in some cases even in the presence of strong magnetic interactions above a critical gyration radius the systematic particle scape on the wake may be observed. As it has been shown [38] through numerical simulations with hydrodynamic interactions, the mechanism of wake

Fig. 23. Snapshots obtained from experimental observations of a spherical blob during sedimentation. The experiments were carried out in the Fluid Mechanics Laboratory of Complex Flows — VORTEX, at the University of Brasilia. (a) Blob of non-magnetized particles with varying topology due to trajectory fluctuation produced by inter-particle hydrodynamic interactions. (b) Preservation of the blob structure when the particles are magnetized. The dipole–dipole magnetic attractive forces are strong enough to suppress the fluctuations induced by interparticle hydrodynamic interactions.

164

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

formation is particle–particle long range hydrodynamic interactions. In the physical reality we will observe both magnetic and hydrodynamic interactions and cannot dissociate these mechanisms or investigate magnetic interactions separately. We may then propose the existence of a critical dimensionless aggregate radius RB/a that will be a function of the following parameters: RB ¼ f ðϕB ; φm Þ: a

ð49Þ

In Eq. (49) ϕB is the hydrodynamic interaction parameter while φm is responsible for accounting the intensity of particle–particle magnetic interactions that would try to maintain the shape of the initially spherical blob. In addition, we show in Fig. 24 a qualitative comparison between the microstructures of a dilute non-Brownian suspension obtained from dynamic numerical simulations (a) and experimental observations (b). The suspensions were studied in the absence of an applied magnetic field for a particle volume fraction ϕ = 0.01, Pe = 103, λ = 1, St = 0.1. The experimental picture was obtained using maghemite particles with average diameter of 100 μm immersed in mineral oil with viscosity ρf = 0.085 Pa ⋅ s (see [30], for details). Both, numerical and experimental results show presence of small clusters throughout the suspension with a clear predominance of dimmers. The simulations have shown a typical behavior of dilute magnetic suspensions in a good agreement with the experimental observations. 5. Final remarks In this work we have presented a numerical method for simulating Brownian and non-Brownian suspensions of magnetically interacting particles. Some features usually neglected in standard numerical simulations, such as particle inertia and particle rotation, were considered in the present work. This work has examined the magnetic effects on the suspension's magnetization due to torque–force dipole–dipole interparticle interactions and also the influence of an applied magnetic field. The calibration of the distance in which magnetic interactions should be considered in order to avoid numerical instabilities allowed the use of a fixed time-step through the suspension time developing.

(a)

In this way we were able to perform low cost computations of dilute magnetic suspensions (i.e. ϕ b 0.05) without periodic dipole–dipole magnetic interactions. For more concentrated suspensions we have used a hybrid method computing torque dipole–dipole interactions, which decays like 1/r3, with periodic Ewald sums and force dipole– dipole interactions in a non-periodic way using a minimum image method. This was a much more economical method than directly using the full periodic sums for force and torque dipole–dipole interactions and surprisingly producing very accurate values of macroscopic properties such as the equilibrium magnetization. The present numerical simulations were validated by comparing our numerical predictions of the suspension magnetization with the well-established asymptotic solutions O(ϕ3) developed by [42]. The tests have shown that the numerical values of the suspension's magnetization were always in a very good agreement with the theoretical models. However, we have observed that the numerical values of the steady state magnetization were slightly higher than the theoretical asymptotic solutions at particle volume fractions around 0.1. This discrepancy has been physically interpreted based on the cluster structures formed during the suspension time evolution process. We have seen that even in dilute regimes the suspension's microstructure is dominated by small clusters, typically with two or three particles. These structures tend to interact with neighboring particles and with closer clusters, increasing substantially the suspension magnetization. This complex dynamics is certainly not captured by the asymptotic theories even at higher orders. We have also shown the potential and application of our dynamic simulations of magnetic suspensions for predicting structure formation and possible breakup of these structures. The anisotropic long chain formation in the direction of the applied field was well captured by our simulations. The formation of smaller isotropic blobs was observed as being more related with dipole–dipole interparticle interactions in the absence of a magnetic field. Finally, when investigating the time evolution process of a sedimenting blob of magnetic particles we have concluded that the short range nature of magnetic forces was sufficiently strong to attenuate any fluctuations due particle hydrodynamic interactions preserving in this way the magnetic blob topology. We have not observed systemic particle scape from the blob domain as the particles are magnetized. This behavior was observed experimentally and confirmed by our numerical simulations.

(b)

Fig. 24. Aggregate formation on a non-Brownian dilute magnetic suspension. (a) Particle distribution given by numerical simulations (a) and (b) experimental observations of a suspension time developing of magnetic composites with an average diameter of 100 μm [30]. The particle volume fraction was about 0.01 and Pe ≫ 1.

R.G. Gontijo, F.R. Cunha / Powder Technology 279 (2015) 146–165

Acknowledgments The authors would like to thank the Brazilian agencies CAPES and CNPq for their partial support of this work. References [1] R.H. Davis, A. Acrivos, Sedimentation of noncolloidal particles at low Reynolds number, Annu. Rev. Fluid Mech. 17 (91) (1985) 91–118. [2] A.E. Boycott, Sedimentation of blood corpuscles, Nature 104 (1920) 532. [3] J.F. Richardson, W.N. Zaki, Sedimentation and fluidization: I, Trans. Inst. Chem. Eng. 32 (1954) 35–52. [4] O.A. Ladyzhenskaya, The Mathematical Theory of Viscous Incompressible Flow, Gordon & Breach, New York, 1969. [5] G.K. Batchelor, Sedimentation in a dilute dispersion of spheres, J. Fluid Mech. 52 (1972) 245–268. [6] G.K. Batchelor, Sedimentation in a dilute polydisperse system of interacting spheres. Part 1. General theory, J. Fluid Mech. 119 (1982) 379–408. [7] R.E. Caflisch, J.H.C. Luke, Variance in the sedimentation speed of a suspension, Phys. Fluids 283 (1985) 759–760. [8] E. Guazzelli, E.J. Hinch, Fluctuations and instability in sedimentation, Annu. Rev. Fluid Mech. 43 (2010) 97–116. [9] H. Nicolai, E. Guazzelli, Effect of the vessel size on the hydrodynamic diffusion of sedimenting spheres, Phys. Fluids 7 (1) (1995) 3–5. [10] P.N. Segre, E. Herbolzeimer, P.M. Chaikin, Long-range correlations in sedimentation, Phys. Rev. Lett. 79 (13) (1997) 2574–2577. [11] R.E. Rosensweig, Directions in ferrohydronynamics, J. Appl. Phys. 57 (1) (1985) 4259–4264. [12] V. Bashtovoi, P. Kuzhir, A. Reks, Capillary ascension of magnetic fluids, J. Magn. Magn. Mater. 252 (2002) 265–267. [13] F. Larachi, M.C. Munteanu, Varying gravity force using magnetic-field emulated artificial gravity: application to cocurrent gas–liquid flows in porous media, Ind. Eng. Chem. Res. 2010 (49) (2010) 3623–3633. [14] N.A. Brusentsov, T.N. Brusentsova, E. Yu Filinova, Yu.A. Pirogov, D.A. Kupriyanov, A.I. Dubina, M.N. Shumskikh, L.I. Shumakov, N. Ya Yurchenko, E.A. Anashkina, A.A. Shevelev, Magnetohydrodynamic thermochemoterapy of malignant tumors with nanopreparations with magnetic resonance monitoring, Pharm. Chem. J. 42 (No. 4) (2008). [15] R.G. Gontijo,, F.R. Cunha, Experimental investigation on thermo-magnetic convection inside cavities, J. Nanosci. Nanotechnol. 12 (2012) 9198–9207. [16] C. Rinaldi, A. Chaves, S. Elborai, X. He, M. Zhan, Magnetic fluid rheology and flows, Curr. Opin. Colloid Interface Sci. 10 (2005) 141–157. [17] S. Odenbach, Colloidal magnetic fluids: basics, development and application of ferrofluids, Lect. Notes Phys. 763 (2009) (Springer, Berlin Heidelberg). [18] J.Y. Hristov, Fluidization of ferromagnetic particles in a magnetic field. Part 1: the effect of field orientation on bed stability, Powder Technol. 87 (1996) 59–66. [19] J.Y. Hristov, Magnetically assisted gas solid fluidization in a tapered vessel: first report with observations and dimensional analysis, Can. J. Chem. Eng. 86 (2008) 470–492. [20] F.R. Cunha, R.G. Gontijo, Y.D. Sobral, Symmetry breaking of particle trajectories due to magnetic interactions in a dilute suspension, J. Magn. Magn. Mater. 326 (2013) 240–250. [21] F.R. Cunha, Y.D. Sobral, R.G. Gontijo, Stabilization of concentration waves in fluidized beds of magnetic particles, Powder Technol. 241 (2013) 219–229. [22] M. Sheikholeslami, M. Gorji-Bandpyb, Free convection of ferrofluid in a cavity heated from below in the presence of an external magnetic field, Powder Technol. 256 (2014) 490–498. [23] M. Sheikholeslami, M. Gorji-Bandpyb, D.D. Ganji, Lattice Boltzmann method for MHD natural convection heat transfer using nanofluid, Powder Technol. 254 (2014) 82–93. [24] M. Sheikholeslami, D.D. Ganji, Ferrohydrodynamic and magnetohydrodynamic effects on ferrofluid flow and convective heat transfer, Energy 75 (2014) 400–410. [25] M. Zahn, D.R. Greer, Ferrohydrodynamic pumping in spatially uniform sinusoidally time-varying magnetic fields, J. Magn. Magn. Mater. 149 (1995) 165–173. [26] M.I. Shiliomis, K.I. Morozov, Negative viscosity of ferrofluid under alternating magnetic fluid, Phys. Fluids 6 (1994) 2855–2861. [27] B.U. Felderhof, Flow of a ferrofluid down a tube in an oscillating magnetic field, Phys. Rev. E 64 (2001) 021508.1–021508.7.

165

[28] K.M. Jansons, Determination of the constitutive equations for a magnetic fluid, J. Fluid Mech. 137 (1983) 187–216. [29] F.R. Cunha, E.J. Hinch, Shear-induced dispersion in a dilute suspension of rough spheres, J. Fluid Mech. 309 (1996) 211–223. [30] F.R. Cunha, H.L.G. Couto, Transverse gradient diffusion in a polydisperse dilute suspension of magnetic spheres during sedimentation, J. Phys. Condens. Matter 20 (2008) 204129. [31] F.R. Cunha, H.L.G. Couto, On the influence of the hydrodynamic interactions on the aggregation rate of magnetic spheres in a dilute suspension, J. Magn. Magn. Mater. 323/1 (2008) 77–82. [32] R.G. Gontijo, F.R. Cunha, Magnetic-induced migration in a sedimenting suspension of magnetic spherical particles, J. Nanosci. Nanotechnol. 12 (2012) 9286–9294. [33] Ph. Depondt, J.C.S. Lévy, F.G. Mertens, Vortex polarity in 2-D magnetic dots by Langevin dynamics simulations, Phys. Lett. A 375 (2010) 628–632 (2011). [34] H. Zhenghua, L. Xiang, L. Huilin, L. Guodong, H. Yurong, W. Shuai, X. Pengfei, Numerical simulation of particle motion in a gradient magnetically assisted fluidized bed, Powder Technol. 203 (2010) 555–564. [35] X. Deng, J.V. Scicolone, R.N. Davé, Discrete element method simulation of cohesive particles mixing under magnetically assisted impaction, Powder Technol. 243 (2013) 96–109. [36] S. Wang, Z. Sun, X. Li, G. Jinsen, L. Xingying, Q. Dong, Simulation of flow behavior of particles in liquid solid fluidized bed with uniform magnetic field, Powder Technol. 237 (2013) 314–325. [37] F.R. Cunha, G.C. Abade, A.J. Souza, E.J. Hinch, Modeling and direct simulation of velocity fluctuations and particle-velocity correlations in sedimentation, J. Fluids Eng. ASME 124 (2002) 957–968. [38] G.C. Abade, F.R. Cunha, Computer simulation of particle aggregates during sedimentation, Comput. Meth. Appl. Mech. Eng. 196 (2007) 4597–4612. [39] A.J.C. Ladd, Dynamical simulations of sedimenting spheres, Phys. Fluids A 5 (2) (1993) 299-10. [40] A.J.C. Ladd, Numerical simulation of particulate suspensions via a discretized Boltzmann equation. Part 2. Numerical results, J. Fluid Mech. 271 (1994) 311–339. [41] P.P. Ewald, Die Berechnung Optisher und Elektrostatischer Gitterpotentiale, Ann. Phys. Lpz 64 (1921) 253–287. [42] A.O. Ivanov, O.B. Kuznetsova, Magnetic properties of dense ferrofluids: an influence of interparticle correlations, Phys. Rev. E 64 (2001) 041405.1–041405.12. [43] C.W. Oseen, Uber die Stokessche Formel und uber die verwandte Aufgabe in der Hydrodynamik, Ark. Matematik, Astronomi och Fysik 6 (No.29) (1910) 1910. [44] A.B. Basset, A Treatise on Hydrodynamics, Deighton Bell, Cambridge, United Kingdom, 1888. [45] Y.D. Sobral, T.F. Oliveira, F.R. Cunha, On the unsteady forces during the motion of a sedimenting particle, Powder Technol. 178 (2007) 129 141. [46] F.R. Cunha, Hydrodynamic dispersion in suspensions, Department of Applied Mathematics and Theoretical Physics, Cambridge University, Cambridge, UK, 1995 (Ph.D Thesis). [47] D.V. Berkov, L.Yu. Iskakova, A.Yu. Zubarev, Theoretical study of the magnetization dynamics of nondilute ferrofluids, Phys. Rev. E 79 (2009) 021407. [48] J.F. Brady, G. Bossis, Stokesian dynamics, Annu. Rev. Fluid Mech. 6 (1988) 111–157. [49] D.A. McQuarrie, Statistical Mechanics, Harper & Row, 1976. [50] A. Einstein, Investigations on the Theory of the Brownian Movement, Dover, New York, 1956. [51] C.W.J. Beenakker, Ewald sum of the Rotne–Prager tensor, J. Chem. Phys. 85 (1986) 1581–1582. [52] J.L. McWhirter, G.N. Patey, Nonequilibrium molecular dynamics simulations of a simple dipolar fluid under shear flow, J. Chem. Phys. 117 (6) (2002) 2747–2761. [53] Z. Wang, C. Holm, H.W. Mller, Boundary condition effects in the simulation study of equilibrium properties of magnetic dipolar fluids, J. Chem. Phys. 119 (2003) 379–387. [54] D.V. Berkov, N.L. Gorn, Stochastic dynamic simulations of fast remagnetization processes: recent advances and applications, J. Magn. Magn. Mater. 290–291 (2005) 442–448. [55] N.A. Usov, Yu.B. Grebenshchikov, Superparamagnetic relaxation time of a singledomain particle with a nonaxially symmetric double-well potential, J. Appl. Phys. 105 (2009) 043904. [56] R.G. Gontijo, Micromechanics and Microhydrodynamics of Magnetic SuspensionsPh.D Thesis University of Brasília, Brazil 2013, p. 262 (in Portuguese). [57] J.M. Nitsche, G.K. Batchelor, Break-up of a falling drop containing dispersed particles, J. Fluid Mech. 340 (1997) 161–175.