Aggregation and flow behavior of magnetic particles in microchannel flow governed by a transverse magnetic field

Magnetic particles (MPs) can be used as carriers for intentionally transporting medium with the effect of external magnetic field. To reveal the aggregation and transportation regularity of magnetic particles in microchannel flow, the three-dimensional (3D) and two-phase lattice Boltzmann method combined by both the immersed boundary method and the discrete element method is established and the investigation on the evolution of chain-like structure composed of MPs in a microchannel flow governed by transverse uniform magnetic field is investigated. It is revealed MPs can maintain the chain-like structure in the channel flow due to the magnetic dipole force, and the final postures of MP chains are determined by the coupling effect of the base fluid and the chains. The average velocity of the MPs moving in the form of chains is smaller than that moving as discrete elements and the velocity of near-wall fluid is increased. Being a cohesive force, if magnetic dipole force is not large enough, the chain-like structure will be broken due to the push of base fluid with velocity gradient. It is also found that there exists a critical value of the intensity of magnetic dipole force for MPs to maintain the chain-like structure in a microchannel flow with the given inlet velocity. This work demonstrates that, under transverse uniform magnetic field, the discrete MPs can be transported as an entirety in microchannel flow on the promise that the magnetic dipole force is large enough to hold the MPs in the flow field with velocity gradient.


Introduction
Current challenges in healthcare, biology, chemistry, etc require integrated platforms, e.g. microfluidic chips, which can handle complex samples in a multiplexed fashion involving minimal cost and training [1]. The manipulation of motion, concentration and disaggregation of nano-or micro-particles plays an important role in microfluidic system, such as that in microelectromechanical systems (MEMS/NEMS) [2][3][4]. In magnetic field based MEMS/NEMS devices, magnetic particles (MPs) can be used as carriers for intentionally transporting medium with the effect of external magnetic field. Besides the common physical properties, such as dielectric constant, electrical and magnetic properties of MPs [5,6], the migration of MPs in base fluid has been paid increasing attentions in these years. Comparatively more contributions are about the flow of magnetorheological fluid affected by magnetic field, among which the mass flow and heat transfer problems have been investigated continuously [7][8][9][10]. These work show that the ferrite, Al 2 O 3 and CuO nanoparticles, spherical or non-spherical [11], can all be manipulated by magnetic field. On the other hand, these contributions are based on the lattice Boltzmann method (LBM) which provides a highly efficient solution to the nanofluid flow problems.
Satoh introduced the viscosity-modifying method for generating random forces in lattice Boltzmann simulations [12]. Later the hybrid-type method of lattice Boltzmann and Brownian dynamics as a simulation technique was developed for a magnetic particle suspension [13]. Recently the LBM has been developed and achieved considerable success in the simulation of unsteady flow with arbitrary boundary shapes [14]. The Koo-Kleinstreuer-Li (KKL) model provides a strategy to reveal the Brownian motion influence on nanofluid

Numerical simulation method
The lattice Boltzmann equations The lattice Boltzmann method, which is valid for estimating the drag force on spherical particle [27], is employed to investigate the microchannel flow of carrier fluid containing MPs. The equation of LBM was derived by He et al [28] with the BGK collision model, which holds where f is the density population function, f eq is the corresponding equilibrium function, τ f is the momentum relaxation time, c k is the lattice velocity, and the subscript k denotes the lattice direction (figure 1); Δx and Δt are the lattice spacing and time step, respectively; r denotes the coordinate vector; F is the external force on the lattice point. For a three-dimensional flow, the D3Q19 lattice with 19-speed, as shown in figure 1, is appropriate and can be used in this work.
The local equilibrium distribution function is obtained from the Maxwell-Boltzmann equilibrium distribution, which is a function of the density and the flow velocity, as shown by where c is the sound speed, ω k is the weighting coefficient and ω 0 =1/3, ω k =1/18 for k=1 to 6 and ω k =1/36 for k=7 to 19. The lattice velocity for the D3Q19 lattice is Model of two-phase flow through a 3D microchannel Figure 2 illustrates a model of two-phase flow through a rectangular microchannel with the dimensions of X, Y and Z in the Cartesian coordinate system of x-y-z. The base fluid, which is considered as a Newtonian fluid, flows into the channel from the inlet boundary at velocity, v in , and outflows from the opposite boundary which is an open boundary. The MPs, being the solid phase, is driven by base fluid and affected by a transverse uniform magnetic field.
The left side of the channel is set as a velocity inlet while the right side is an open boundary. The upper and lower walls are set as rebound format boundaries [29]. To reduce amount of computation, only a part of microchannel flow is considered and the two surfaces (y=1 and y=Y) are set as periodic boundaries. In this work, we mainly concern the transportation of MPs in the microchannel flow governed by transverse uniform magnetic field, so Y is valued as two or three times of the diameter of MP. At the beginning of simulation, MPs are seeded in the local region near the inlet boundary and furthermore, in order to shorten simulation time, the y-coordinates of all MPs are set as y=0.5Y.

The immersed boundary method (IBM)
Compared with the complex models for the hydrodynamic force and torque experienced by a spherical particle [27,30], the IBM, which is relatively simpler, is used in this work to deal with the boundary interaction between MP and base fluid. As shown in figure 3, some lattice points are inside the solid boundary while some are outside the solid boundary.
In IBM the collision between the lattice points inside and outside the solid boundary can be performed in the Ladd's half-step bouncing format, i.e., collision takes place at the midpoint between a pair of lattice points. Thereby in the simulation the circular solid boundary is substituted by the simplified boundary denoted by the symbol '☐'. The half-step bouncing format can be illustrated by figure 4,where x in and x out are the two lattice points inside and outside the solid boundary respectively. x b is the midpoint between x in and x out , which presents the boundary point. At time t the populations at points x in and x out are ¢ f k (x in , t) and ¢ -f k (x out , t), respectively, then after Δt/2, the collision at x b gives After another half time step of Δt, the new populations at x in and x out become The resultant force, F S , and torque, T S , arising from surrounding fluid can be respectively calculated by where the subscript, δ, is the number of points on the simplified boundary denoted by '☐', i presents the number of related lattice directions, and x c presents the position of particle center.

External forces on MPs
The external force exerting on MP, F MP , is a resultant force of the Stokesian force arising from base fluid, F S , (equation (6)), magnetic dipole forces among MPs, F m , the Brownian force, F B , the Van der Waals, F VW , and the gravity and buoyancy of MPs, etc. In addition to the torque arising from base fluid, T S , the torques generated by the magnetic dipole potential and magnetic field, T m and T mf , also influence the motion of MPs.
The magnetic dipole force between two MPs, F m , can be calculated using equation (8) [31].
n n t n t n t t n t n n t n   Brownian motion is significant for submicron size particles. This random motion is modeled as a white noise process, which is dependent on the physical properties of the medium, such as the temperature and density. The model for Brownian motion is given mathematically by equation (9) [32].
where m p represents the mass of particle, n i (t) are Cartesian components of vector n B (t), μ is the dynamic viscosity of the fluid, T is the absolute temperature, k B is the Boltzmann constant, S is the ratio of the density of particle to the density of fluid, ρ is the density of the fluid, and G i is a Gaussian random variable with zero mean and unit variance. C c is the Cunningham slip correction factor, which is given as where R is the radius of the particle and λ is the mean free path of the medium.
The van der Waals force, F VW , between two spheres is a function of separation r ab , which yields where A is the Hamaker constant for nanoparticle in water, and R a and R b are radii of two spheres. F VW is along the line connecting two MPs. Other external forces, such as the gravity and buoyancy are neglected in the simulation.
T m and T h can be calculated according to equations (11) and (12), respectively.
The DEM is applied to analyze the motion of MPs, i.e., the velocity and displacement of each MP are calculated using the explicit central difference scheme. The normal and tangential contact forces between two MPs are processed based on the Hertz soft sphere model [34] and the Mindlin-Deresiewicz contact theory [35], respectively. The two-phase flow of MRF in microchannel is an unsteady process and the simulation includes two main procedures, one is the LBM simulation of base fluid flow by taking the MPs as static solid obstacles and the other is the DEM simulation of MPs which will pass through the lattices. The convergence of LBM simulation of base fluid flow is judged according to equation (13) for velocity.

Grid independency and program validation
The main parameters used in the LBM simulation are nondimensionalized as where Δt is the relaxation time, ρ b and ρ p are the densities of base fluid and MP respectively. Both the dimensionless lattice side length, Δx * , and relaxation time, Δt * , equal to 1.
Since the IBM and half-step bouncing format are used in the simulation, the key grid parameter that influences the precision is the lattice number with respect to MP diameter. A LBM code for a microchannel flow with a MP being fixed in the channel center is established to study the grid independency, and the main parameters include d p =1×10 −4 m, X=10d p , Y=3d p , Z=4.5d p and v in * =0.1. The dimensions, such as d p , X, Y and Z will not change, but the lattice length, Δx, is reduced continuously to see the change of Stokesian force, F S , with the increasing d p /Δx. It is found that when d p /Δx10, the variation magnitude of F S is slight, therefore in the following simulations, the lattice length is conformed to the requirement of d p /Δx=10. Also based on the fixed single MP model, the Stokesian force, F S , exerting on the fixed MP were sampled at different inlet velocities and compared with that obtained by using the commercial code FLUENT, as show in figure 5. It can be seen that the dimensionless Stokesian force, F S * , obtained by using the program developed in this work is almost linearly increasing with the increase of inlet velocity and it is a little smaller than that obtained by using FLUENT. But the drag force obtained by using the Stokesian equation (F S =3πμd p Δv) is much larger than the two numerical results. For the complex two-phase flow field, the Stokesian equation is not precise enough to estimate the drag force on solid particles.

Results and discussion
Simulation conditions where N is number of MPs and ΔX=x 2 −x 1 . Lowering Re in 3D simulations by a factor of R (the dimensionless particle radius in LB simulations) resulted in a close match between numerically computed and experimentally obtained particle velocities [36], therefore the largest inlet velocity is set as 0.1 in the following simulations, corresponding to Re=30.
Under the simulation conditions illustrated above, the motion of MPs mainly occurs in the plane of y=Y * /2, so in the following results, the 2D flow fields corresponding to the x-z surface when y=Y * /2 will be extracted and shown.

Transportation of single ideal chain of MPs
Under uniform magnetic field, the MPs can form chains parallel to the line of magnetic force [37], which has been observed by many researches. Thereby in the following simulations, to save simulation time, we directly set a straight chain of MPs in the channel near the inlet.
To deeply understand the motion behaviors of the chain-like structure of MPs driven by a microchannel flow of base fluid, we firstly immersed a straight chain of MPs near the inlet of the channel with the length of X * =400 and depth of Z * =80, as shown in figure 6.
Since the chain is located ideally in the base fluid at the beginning, it can keep the posture being perpendicular to the flow direction. The integrity of chain is ensured due to the magnetic dipole forces which can adjust the relative positions of MPs in order to match the Stokesian forces arising from base fluid. The final velocity of the chain, v * ch , is found to be 0.049 and it is a little smaller than the inlet velocity of 0.05. Since the two ends of the chain are close to the two walls, the velocity of the near-wall fluid is increased. In general on this condition, the chain of MPs is almost straight and evens the surrounding velocity distribution obviously. Figure 7 shows the motion of the single chain with the doubled inlet velocity of v * in =0.1. It can be seen that the increase of the inlet velocity not only increases the chain velocity but also makes the chain be curved. The final velocity of the chain, v * ch , is 0.09 and it is also smaller than the inlet velocity. On the larger inlet velocity condition, the chain is not straight anymore and is not perpendicular to the flow direction. Due to the calculation accuracy error, the flow field is not symmetrical about the central line of the channel and the chain   The effect of the chain on the velocity distribution in base fluid can be seen in the figure 8 which shows the velocity profiles at different cross-sections when t * =900. Curve-1 presents the velocity at the cross-section of x * =180 which locates at an upstream position relative to the chain. The overall profile shows that it is not obviously disturbed by the chain except for that it is not symmetrical about the central line. Curve-3 shows the velocity profile of base fluid near the channel inlet, and in this region, the velocity at the center of channel is little smaller than that in the adjacent fluid. The velocity profile of the downstream base fluid is denoted by Curve-2 and there exists a short platform in the curve. The most especial velocity profile locates in the upstream flow just near the chain, as shown by Curve-4. It has a wide region with uniform velocity which means the chain creates an average velocity distribution of the bulk flow around itself, and meanwhile there exist the larger velocity gradients in the near-wall fluid.
On low Reynolds number condition, it is the Stokesian force, F S , that pushes or drags MP to move together base fluid. Meanwhile the inherent magnetic dipole forces among MPs hold the chain-like structure. The mechanical model for any couple of MPs in a formed chain in base fluid with velocity gradient is illustrated in figure 9.
If there are only two MPs which are relatively at rest, there must exist such a relationship between F S and F m as shown by Considering F mab and F′ mab are a couple of action and reaction and F S can be estimated by F S =3πμΔvd p , we have where v MPs is the velocity of the two MPs, v a and v b are the velocities of local fluids interact with the two MPs respectively. Thereby It can be seen that the velocity of the MP group falls in between v a and v b . Therefore F Sa is in the opposite direction relatively to F Sb . So it can be deduced that for a curved chain of MPs with a fixed posture, its velocity must be smaller than the mainstream speed of the base fluid while larger than the velocity of near-wall fluid. As a result, the MPs in the middle of the chain are pushed by the base fluid while the near-wall fluid is pushed by the MPs at the ends of the chain.
From equation (14) it can also be seen that if the velocity difference between two adjacent MPs increases, their relative displacement along the x-axis will increase and F m increase correspondingly. Since there exists the largest x-component of F m [38], once the relative displacement of the two adjacent MPs is larger than the critical distance, the two MPs will depart from each other.
According to the force analysis, it can be revealed that the two factors that influence the chain-like posture are the velocity gradient of base fluid and the intensity of magnetic dipole force (or the magnetic moment). To investigate the relationship between the two factors, the intensity of magnetic dipole force is reduced to see whether the chain can still maintain its intensity under the same inlet velocity. As a result of the intensity of magnetic dipole force being decreased to be 0.5M V , the curvature of the chain is increased due to the increased relative displacement of MP couples in the first 100 time steps, as shown in figure 10. Along with the development of the flow, the velocity gradient increases, and two MPs at the ends of the chain break from the chain and fall behind because there the velocity gradient is the largest (t * =300). In the flowing steps, without the drag effect of the near-wall fluid, the curvature of the short chain composed of five MPs is reduced and almost move steadily along the flow direction.
If the intensity of magnetic dipole force is reduced further to be 0.95×10 5 A·m −1 , the chain lost its integrity and disperse system behaves like a common nanofluid, as can be seen in figure 11. At t * =900, only two MPs form a group which locate at the channel center where the velocity gradient is the smallest.
The results shown by figures 6, 10 and 11 indicate that when the inlet velocity of base fluid is constant, the intensity of magnetic dipole force plays the dominant role in the integrity of chain-like structure. Figures 12, 13 and 14 show the evolution of chain-like structure formed by randomly dispersed MPs, corresponding to different intensities of magnetic dipole force such as 0.19×10 5 , 0.95×10 5 and 1.9×10 5 A·m −1 respectively. In these cases the inlet velocity is 0.1 and the length of the channel is prolonged to be 640 and the simulation will be terminated once any MP approaches x * =580. Firstly it can be seen that the overall velocity of the MPs is reduced with the increase of M V . It takes 1500 time steps for the MPs in the first case (M V =0.19×10 5 A·m −1 ) to get to x * =580 while it costs 2700 time steps for the third case (M V =1.9×10 5 A·m −1 ). In figure 12 no obvious chains can be found in the channel, and as the motion goes on, the MPs diffuse continuously due the velocity gradient of the base fluid along the y-axis.

Transportation of randomly dispersed MPs
When the intensity of the magnetic dipole force is increased to be 0.95×10 5 A·m −1 ( figure 13), some chains form in the base fluid but it can be seen the long chain which touches the upper and lower wall boundaries formed before t * =900 (denoted by the dashed rectangle) breaks into two parts before t * =1500. Due to the drag effect of the near-wall fluid, those short chains that attach the wall boundaries are not perpendicular to the flow direction. After t * is larger than 1500, no obvious change of the entire structure is observed.
As the intensity of magnetic dipole force is large enough, almost all the MPs contribute to the formation of chains, as can be seen in figure 14. After the chains have fully developed, their postures and relative positions will  not change obviously as can be seen from the postures from t * =1900 to t * =2700. Compared with the single chain shown in figure 7, the chains shown below is more curving and that is mainly caused by coupling effect of the based fluid and MP chains.  where N′ is the number of MPs locating in the volume of ΔX×Y×Z which moves with the leading MP. If other MPs can not move together with the leading MP, so ε′ will decrease. Like shown in figure 14, there may be a few individual MPs that can not catch up with the main part which is composed of chains with different lengths, so those individual MPs can be ignored while calculating ε′. Figure 16 illustrates the variation of ε′ with time steps corresponding to different values of M V . From the variation of ε′ it can be seen that if the intensity of magnetic dipole force is not strong enough, MPs will disperse continuously due to their velocity difference. When M V =1.52×10 5 A·m −1 , there exists a small decrease for ε′ in the former 500 time steps and then ε′ remains constant and that means most of the MPs can move in a entirety. Thereby the critical intensity of magnetic moment can be determined, i.e., M V =1.52×10 5 A·m −1 . Similarly, the critical diameter and inlet velocity can also be determined according to the simulation result.

Limitation of the two-phase LBM simulation method
The simulation results indicate that due to the velocity gradient along the transverse direction with respect to the bulk flow direction, the MPs will disperse continuously if there are no strong enough magnetic dipole forces among MPs. So for some applications MPs are used as carriers, the precise control of the number of MPs can not be realized easily. But if the integrity of the chain-like structure can be ensued during the transportation process, MPs can be moved in the form of cluster. As described in the former sections, the final structure composed by the MPs is codetermined by several factors, such as the inlet velocity and the intensity of magnetic dipole force.
In this work, we just reduce the intensity of magnetic force to enlarge the difference between F S and F m but not increase the inlet velocity in consideration of the limitation of LBM simulation. Since the velocity field is a result of fluid-structure interaction, to determine the final chain-like structure of MPs still depends on numerical simulation. Therefore, it is very difficult to give a general analytical criterion for the posture of MPs under different flow conditions. For a pipe or channel flow, the Reynolds number, Re, is always taken as the similitude law that ensures simulation results can be in accordance with that in prototype. For the two-phase flow concerned in this work, the main parameters involved in the situation include . By taking Z, ρ and v as the fundamental parameters, four dimensionless Π groups can be obtained as  According to the similitude theorem, the LBM simulation results of a model can be used to predict the flow situation of a real microchannel flow of MRF on the condition that all the four Π groups for the model are equal to that of the prototype, respectively. Π 1 presents the ratio of the particle diameter to the characteristic dimension of channel, and this law can be confirmed by scaling the parameters of d p and w using the same ratio while modeling. Π 2 is the ratio of the density of MP to that of the base fluid, that can be easily satisfied by using the same materials with that of prototype. Π 3 requires that the Reynolds number of real flow should equal to that of the model and in this work, * n = v Z Re in which is termed as the lattice Reynolds number. In LBM simulation v in * is valued from 0 to 0.1, and it is not directly related to actual inlet velocity. If v in * is determined, the dimensionless kinematic viscosity ν can be determined. The expression of Π 4 indicates that the magnetic dipole force F m exerting on any MP should be determined when Re, ν and ρ b are given. But F m is not a constant and will vary with the posture of MPs, it is very difficult to model a microchannel flow of MRF in terms of Π 4 .
Therefore, it is impossible to use the results of one flow model to predict the characteristic of another flow. That is to say, we have to establish a new model for a new flow. If the difference between MP diameter and characteristic length of channel is increased, the grid number will be increased that results in the increased computational effort. For a real MRF flow in microchannel, the diameter of MPs is much smaller than that used in this work, and the number of MPs in a chain will increase correspondingly. As explained in section 3.3, it is very difficult to establish a complete flow model to investigate the chains of MPs. Thereby, on the premise that the velocity distribution in a microchannel flow of MRF is known, some Couette flow models with different shear rates can be established to simulate the motion of MPs in different regions of the flow. Then MPs can be modeled according to the real dimension and the structure evolution of long chains can be predicted by grouping all the results of Couette flow models.

Conclusions
In this work, the 3D and two-phase LBM integrated by the immersed boundary method and the discrete element method is developed to simulate the aggregation and motion regularity of MPs in microchannel flow governed by a transverse uniform magnetic field.
The ideal single chain of MPs simulation results show that the chain-like structure is affected by both inlet velocity and intensity of magnetic dipole force. Once the Stokesian force is larger than the resultant magnetic dipole force on a MP, the MP will depart from the chain and then the chain will be shortened. This regularity is also valid for the multi-chain system, and the final postures of the chains are determined by the coupling effect of the base fluid and MPs. The velocity of the entire chain is a little smaller than the bulk velocity but larger than the velocity of the near-wall fluid, that means the chain-like structure hinders the bulk flow of base fluid. The average velocity of MPs will be reduced if they move in the form of transverse chains.
The investigations illustrate that, under transverse uniform magnetic field, the discrete MPs can be transported as an entirety in microchannel flow on the promise that the magnetic dipole force is large enough to hold the MPs in the flow field with velocity gradient. It is found that there exists a critical value of the intensity of magnetic dipole force for the MPs to maintain the chain-like structure in a microchannel flow with given inlet velocity.