Simulations of behavior of magnetic particles in magnetic functional fluids using a hybrid method of lattice Boltzmann method, immersed boundary method and discrete particle method

Accepted Manuscript

Simulations of behavior of magnetic particles in magnetic functional fluids by using hybrid method of lattice Boltzmann method, immersed boundary method and discrete particle method Yasushi Ido , Hirotaka Sumiyoshi , Hiroaki Tsutsumi PII: DOI: Reference:

S0045-7930(16)30123-2 10.1016/j.compfluid.2016.04.019 CAF 3158

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

2 November 2015 6 March 2016 15 April 2016

Please cite this article as: Yasushi Ido , Hirotaka Sumiyoshi , Hiroaki Tsutsumi , Simulations of behavior of magnetic particles in magnetic functional fluids by using hybrid method of lattice Boltzmann method, immersed boundary method and discrete particle method, Computers and Fluids (2016), doi: 10.1016/j.compfluid.2016.04.019

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Highlights  New hybrid code of LBM-IBM-DPM with magnetic interactions for MR fluids is developed.  Three-dimensional microstructures of magnetic particles in MR fluids are visualized  Microstructures in MR fluid flows are shown considering fluid-particle interactions.  Requirement of forming sheet-like structures of chain-like clusters is reported.

ACCEPTED MANUSCRIPT

Simulations of behavior of magnetic particles in magnetic functional fluids by using hybrid method of lattice Boltzmann method, immersed boundary method and discrete particle method Yasushi Idoa*, Hirotaka Sumiyoshib, Hiroaki Tsutsumia a Department of Engineering Physics, Electronics and Mechanics, Graduate School of Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan b Makita Corporation, 3-11-8 Sumiyoshi-cho, Anjyo 446-8502, Japan *corresponding author: [email protected]

CR IP T

ABSTRACT The microstructure formation of magnetic particles in magnetorheological (MR) fluids has an important role on the macroscopic properties such as viscosity. To simulate the behavior and microstructure formation of magnetic particles, hybrid simulations were conducted by combining the lattice Boltzmann, immersed boundary and discrete particle methods, and by taking account of magnetic interactions between magnetic particles under an external magnetic field. The microstructure formation of magnetic particles in an MR

AN US

fluid between two parallel plates, where one plate moves with constant velocity and the plate remains stationary, is calculated. A sheet-like structure appears when relatively strong shear and strong magnetic field are applied. In contrast, a fibrous-like structure moves in the flow when the shear is small. Keywords: MR fluid, Lattice Boltzmann method, Immersed boundary method, Discrete particle method,

M

Cluster formation, Microstructure

1. Introduction

ED

Magnetorheological (MR) fluid is a type of magnetic functional fluid that reacts to a magnetic field. Magnetic functional fluids are regarded as mixtures of a base fluid and magnetic particles. Ferrofluids are

PT

stable suspensions of fine magnetic particles with an average diameter of about 10 nm [1] and MR fluids are fluids mixed with micrometer-size magnetic particles [2-4]. Magnetic particles in an MR fluid form chain-like clusters in the direction of the magnetic field [5-7]. The microstructure formation of magnetic

CE

particle clusters is an interesting topic in the fields of physics, chemistry, and engineering. It is important to know the microstructure for both an understanding of the basic properties of MR fluids and for the development of industrial applications, because macroscopic properties such as viscosity are dependent on

AC

the microstructure in the fluid. In some applications, such as dampers [8-10] and polishing [11], the properties due to the microstructure formation are exploited. However, it is difficult to directly observe the three-dimensional microstructure in a fluid directly by using microscopy because the images obtained from a microscope are two-dimensional. Therefore, simulation is a useful tool for the investigation of such three-dimensional microstructures. Many researchers have studied the microstructure of magnetic particles in MR fluids using computer simulations such as the discrete particle method and Brownian dynamics. Zhu et al. suggested the phenomenological model of nucleation-controlled microstructures in a dilute MR fluid [12]. Mohebi and Jamasbi reported that chain-like clusters of submicron-size magnetic particles form wall-like structures and labyrinthine patterns [13]. Molecular dynamics simulations performed by Gullery and Tao [14] indicated three different structures, the bct (body-centered tetragonal) lattice looking along the direction of the applied magnetic field, chains and the liquid state. Two-dimensional microstructure development was computed by Ly et al. [15]. Ido et al. reported the distribution of both magnetic particles

ACCEPTED MANUSCRIPT

and nonmagnetic particles in MR fluids in the presence of a magnetic field, and it was shown that the nonmagnetic particles are also rearranged slightly in the field direction due to the motion of magnetic particles to form chain-like clusters, and especially the axial direction of capsule nonmagnetic particles could be easily changed to the direction of the magnetic field [16,17]. The distribution of particles in the fluid was calculated using simplified Stokesian dynamics, and the influence of the ratio of the nonmagnetic particle diameter to that of magnetic particles on the distribution and microstructure formation was examined. Sherman et al. conducted real scale simulations to calculate chain structure formation in an MR fluid with more than one million particles [18]. Li and Peng analyzed the microstructure of magnetic particles in MR fluids with different-sized particles and they showed that the ratio of particle size in bidisperse suspensions has an effect on the shear stress [19]. Most of these simulations were performed to

CR IP T

study stable microstructures in static MR fluids. Even when the fluid flow was considered, the flow was given as the steady flow and effect of suspended particles on the fluid was not taken into account. However, to use MR fluids in industrial applications such as dampers, it is important to determine the microstructure of magnetic particles in a flow of MR fluid with consideration of the interactions between particles and the fluid.

In a recent study reported by Zhang et al. [20], a new scheme that involved the combination of three

AN US

different methods, the lattice Boltzmann method [21,22], the immersed boundary method [23, 24] and the discrete element method [25], was developed to simulate the behavior of particles in fluids. Interactions between particles and the base fluid are considered in these simulations, and the fluid motion is determined taking into account the force due to the motion of particles. When the immersed boundary method is adopted, the body in the fluid is treated as an elastic deformable object with high stiffness [23]. The immersed boundary method was developed by Peskin [23] to investigate blood flows of the heart, and this method has

M

been applied to suspension flow problems. The body in the fluid is discretized as small segments with boundary points and each boundary point can be slightly distorted from its original position by the motion of the surrounding fluid. The distortion of the boundary points means that restoring forces of Hook’s law acts

ED

on the body, and these forces are calculated as interaction forces between the body and the fluid. In this paper, a new hybrid code based on a combination of the lattice Boltzmann method, the immersed

PT

boundary method and the discrete particle method, is developed to consider the magnetic interactions between magnetic particles. To evaluate the code, drag force acting on a fixed spherical particle in a uniform flow is calculated and the parameter of the immersed boundary method is determined. The behavior of

CE

magnetic particles in shear flows of MR fluid between two parallel plates in the presence of a magnetic field is simulated. The effects of the volume fraction of magnetic particles, the strength of the shear force and the

AC

intensity of the applied magnetic field on the behavior of magnetic particles are investigated.

2. Simulations

The hybrid simulation scheme, which combines the lattice Boltzmann, the immersed boundary and discrete particle methods with consideration of the magnetic interactions, is developed to simulate the behavior of magnetic particles in MR fluids and take into account interactions between the base fluid and suspended magnetic particles. The motion of the base fluid is simulated using the lattice Boltzmann method, and the motion of magnetic particles in MR fluids is calculated using the discrete particle method. Interactions between the base fluid and the magnetic particles are described using the immersed boundary method. 2.1. Lattice Boltzmann method

ACCEPTED MANUSCRIPT

The motion of the base fluid is calculated using the lattice Boltzmann method [21,22]. Numerical region is discretized using a regular fixed lattice with virtual particles that move toward the lattice with definite velocity. The flow velocity and the density of the fluid are obtained by calculating the time development of the local distribution function of virtual fluid particles on the lattice points. The time development of the local distribution function 𝑓𝑖 (𝑥⃗𝑖 , 𝑡) is derived from the following equation: 1 3 𝑒𝑞 𝑓𝑖 (𝑥⃗ + 𝑐⃗𝑖 𝛥𝑡, 𝑡 + 𝛥𝑡 ) − 𝑓𝑖 (𝑥⃗, 𝑡) = [𝑓𝑖 (𝑥⃗, 𝑡) − 𝑓𝑖 (𝑥⃗, 𝑡)] + 𝑤𝑖 𝑓⃗𝑖𝑒 ⋅ 𝑐⃗𝑖 𝜏 2

(𝑖 = 1, 2, ⋯ , 15) , (1)

CR IP T

where 𝑐⃗𝑖 is the lattice velocity in the i-direction, 𝜏 is the relaxation time, 𝛥𝑡 is the time step, 𝑓⃗𝑖𝑒 is the external force due to the interaction between the fluid and the particles. The first term of right-hand side of Eq. (1) is the collision term and the second term expresses the external force. In the BGK (Bhatnagar-Gross-Krook) model [26] for calculation of the collision term, variation of the distribution 𝑒𝑞

function is determined using the relaxation time 𝜏 and the local equilibrium distribution function 𝑓𝑖 . The local equilibrium distribution function is the distribution function when the stream and collision processes

AN US

of fluid particles are in a state of equilibrium in finite area. The local equilibrium distribution function is obtained from the Maxwell-Boltzmann equilibrium distribution [21], which is a function of the density and the flow velocity, as shown by:

𝑒𝑞

= 𝜌𝑤𝑖 [1 +

3𝑐⃗𝑖 ⋅ 𝑢 ⃗⃗ 9(𝑐⃗𝑖 ⋅ 𝑢 ⃗⃗)2 3𝑢 ⃗⃗ ⋅ 𝑢 ⃗⃗ + − ], 2 4 2 𝑐 2𝑐 2𝑐

M

𝑓𝑖

(2)

where c = √3𝑅𝑇, R is the gas constant, T is the absolute temperature, and wi is the weight function of the

CE

PT

ED

3D15Q model given by:

2 (𝑖 = 0) 9 1 𝑤𝑖 = (𝑖 = 1 to 6) . 9 1 { 72 (𝑖 = 7 to 14) (3)

AC

In this paper, the 3D15Q model is adopted and the lattice velocity, 𝑐⃗𝑖 is defined as: 0 1 0 0 −1 0 0 1 −1 1 1 −1 1 −1 −1 [𝑐⃗0 , 𝑐⃗1 , ⋯ , 𝑐⃗14 ] = [0 0 1 0 0 −1 0 1 1 −1 1 −1 −1 1 −1] . 0 0 0 1 0 0 −1 1 1 1 −1 −1 −1 −1 1

(4)

The relaxation time is set at 𝜏= 0.7 in the present simulations. The density and flow velocity are calculated using the distribution functions from the following equations:

𝜌 = ∑ 𝑓𝑖 , (5)

ACCEPTED MANUSCRIPT 𝑢 ⃗⃗ =

1 ∑ 𝑐⃗𝑖 𝑓𝑖 . 𝜌

(6) 2.2. Immersed boundary method In the immersed boundary method, interactions between the base fluid and particle are calculated for the Lagrangian boundary points arranged on the surface of a spherical magnetic particle. The gravity point of the spherical magnetic particle is the reference point of the immersed boundary method. Figure 1 shows the concept of the immersed boundary method. The reference point of the particle and the Lagrangian boundary points on the surface of the particle are independent from the fixed lattice Boltzmann grids, as shown in Fig.

𝐹⃗𝑗𝑙 = −𝑘𝜉⃗𝑗𝑙 = −

𝜉𝑗𝑙 𝜅 ⃗⃗⃗⃗⃗ , 𝜖𝑑 ‖𝜉⃗⃗⃗⃗⃗ 𝑗𝑙 ‖

CR IP T

1(a). The force acting on the fluid from the Lagrangian boundary point l is given by:

(7)

AN US

where k is the positive spring constant, 𝜅 is the force scale factor, 𝜖𝑑 is the stiffness scale factor, 𝜉⃗𝑗𝑙 is the relative position vector of the position where the boundary point moves with the flow velocity at the boundary point against the position where the boundary point moves with the rigid motion of the particle, i.e., the position where the point moves with velocity and angular velocity of the particle [27]. Figure 1(b) illustrates the relative position vector 𝜉⃗𝑗𝑙 . The Lagrangian boundary points on the spherical particle are

AC

CE

PT

ED

M

arranged using the method proposed by Uhlmann to describe the spherical magnetic particles [28].

(a)

(b) Fig. 1. Concept of the immersed boundary method. The reference point and boundary points are indicated in

ACCEPTED MANUSCRIPT

(a). The left-hand side of (a) is 2D image, while the right-hand side of (a) is a 3D image. In right-hand side of (a), boundary points are arranged using the Uhlmann method. The basic concept and relative position vector are described in (b). 2.3. Discrete particle method The discrete particle method is based on the equations of motion. In the present simulations, a magnetic dipole model for spherical magnetic particles and a simplified Stokesian model [29] are adopted. The motion of j-th magnetic particle with mass 𝑀𝑗 at time t and position 𝑟⃗𝑗 is described by the following equations of motion: 𝑁

𝑙=1

𝑘≠𝑗

𝑁

CR IP T

𝑑 2 𝑟⃗𝑖 𝑟𝑒𝑝 𝑀 𝑀𝑗 2 = ∑ 𝑘𝜉⃗𝑗𝑙 + ∑ 𝐹⃗𝑗𝑘 + 𝐹⃗𝑗 , 𝑑𝑡

(8)

𝑑𝜔𝑗 𝑀 ⃗⃗𝑗𝐻 + ∑ 𝑇 ⃗⃗𝑗𝑘 ⃗⃗𝑗𝑟𝑒𝑝 , 𝐼𝑗 = ∑(𝑋⃗𝑗𝑙 − 𝑥⃗𝑗 ) × 𝑘𝜉⃗𝑗𝑙 + 𝑇 +𝑇 𝑑𝑡 𝑘≠𝑗

AN US

𝑙=1

(9)

where 𝑋⃗𝑗𝑙 is the position vector of the l-th boundary point of the j-th particle, 𝑥⃗𝑗 is the position vector of the 𝑟𝑒𝑝 reference point on the j-th particle, 𝐼𝑗 is the inertial moment of the particle, and 𝐹⃗𝑗 is the repulsive force

M

acting on the particle, which is derived by:

𝑟𝑒𝑝 𝑝 𝐹⃗𝑗 = ∑ 𝐹⃗𝑗𝑘 + 𝐹⃗𝑗𝑤 , 𝑘≠𝑗

ED

(10)

PT

𝑝 where 𝐹⃗𝑗𝑘 is the repulsive force due to contact between the i-th particle and the j-th particle, and 𝐹⃗𝑗𝑤 is the

repulsive force due to contact between the particle and the wall surface. In some previous works, the force from the Leonard-Johns potential [30,31], exponential expression [13,14], or DLVO (Derjaguin, Landau,

CE

Verwev, and Overbeek) theory [16,17] was used to describe the repulsive force, however, the Hertz contact

AC

theory is applied to derive the repulsive force in the present study because we assume that the developed 𝑝 code will be used in simulations for applications such as polishing. The repulsive force 𝐹⃗𝑗𝑘 is given by: 3

1

1

1

𝑝 𝑝 ⃗⃗𝑗𝑘 ⋅ 𝑛⃗⃗) 𝑛⃗⃗ − 𝑘𝑡𝑝 𝛿𝑛2𝛿⃗𝑡 − 𝜂𝑡 |𝛿⃗4 | 𝑉 ⃗⃗𝑓𝑗𝑘 , 𝐹⃗𝑗𝑘 = (−𝑘𝑛 𝛿𝑛2 − 𝜂𝑛 𝛿𝑛4 𝑉 𝑡

(11)

where 𝑛⃗⃗ is the normal unit vector against the contact surface, 𝛿𝑛 is the displacement in the normal direction due to contact between particles, 𝛿⃗𝑡 is the relative displacement vector in the tangential direction ⃗⃗𝑗𝑘 = 𝑉 ⃗⃗𝑗 − 𝑉 ⃗⃗𝑘 is the relative velocity vector, 𝑉 ⃗⃗𝑗 and 𝑉 ⃗⃗𝑘 are the velocity vectors of j-th at the contact point, 𝑉 ⃗⃗𝑓𝑗𝑘 is the relative velocity vector in the tangential direction at the particle and k-th particle, respectively, 𝑉 𝑝

𝑝

contact point, 𝑘𝑛 and 𝑘𝑡 are the coefficients of elastic deformation in the respective normal and tangential directions at the contact point between magnetic particles, which are obtained by:

ACCEPTED MANUSCRIPT √2𝐷𝐸𝑝

𝑝

𝑘𝑛 =

3(1 − 𝜎𝑝2 )

, (12)

𝑝

𝑘𝑡 =

2√2𝐷𝐺𝑝 , 2 − 𝜎𝑝 (13)

where 𝐸𝑝 is the longitudinal elastic modulus of the magnetic particle, 𝜎𝑝 is the Poisson ratio of the magnetic particle, and 𝐺𝑝 is the modules of rigidity of the magnetic particle. 𝜂𝑛 and 𝜂𝑡 are the coefficients of friction in the respective normal and tangential directions at the surface of contact point, and they are

CR IP T

defined by:

𝑝

𝜂𝑛 = 𝛼√𝑀𝑘𝑛 ,

𝑝

AN US

𝜂𝑡 = 𝛼√𝑀𝑘𝑡 ,

(14)

(14)

where M is the mass of the magnetic particle, and α is the non-dimensional constant parameter for the viscous damping force, which is set as α = 2.0 in these simulations. The repulsive force 𝐹⃗𝑗𝑤 , is obtained by: 2

1

𝑝

(15)

ED

M

𝐹⃗𝑗𝑤 = 𝑘𝑛𝑤 𝛿𝑛3 𝑛⃗⃗ + 𝑘𝑡𝑤 𝛿𝑛2 𝛿⃗𝑡 ,

𝑝

where 𝑘𝑛 and 𝑘𝑡 are the coefficients of elastic deformation in the respective normal and tangential

PT

directions at the contact point between the magnetic particle and the wall surface. These coefficients are

AC

CE

given by:

−1

𝑘𝑛𝑤

4√𝐷 1 − 𝜎𝑝2 1 − 𝜎𝑤2 = ( + ) 3 𝐸𝑝 𝐸𝑤

, (16)

𝑘𝑡𝑤 =

8√𝐷𝐺𝑝 , 2 − 𝜎𝑝 (17)

where D is the diameter of the particle, 𝐸𝑤 is the longitudinal elastic modulus of the wall, 𝜎𝑤 is the Poisson ratio of the wall. It is assumed that deformation of the wall in the tangential direction can be neglected. 𝑀 The force 𝐹⃗𝑗𝑘 is the magnetic interaction force between magnetic particles, which is given by: 𝑀 𝐹⃗𝑗𝑘 =

𝑟⃗𝑗𝑘 𝑟⃗𝑗𝑘 3 1 ⃗⃗⃗𝑗 ⋅ 𝑚 ⃗⃗⃗𝑘 ) − 5(𝑚 ⃗⃗⃗𝑗 ⋅ 𝑟⃗𝑗𝑘 )(𝑚 ⃗⃗⃗𝑘 ⋅ 𝑟⃗𝑗𝑘 ) 3 + [(𝑚 ⃗⃗⃗𝑗 ⋅ 𝑟⃗𝑗𝑘 )𝑚 ⃗⃗⃗𝑘 + (𝑚 ⃗⃗⃗𝑘 ⋅ 𝑟⃗𝑗𝑘 )𝑚 ⃗⃗⃗𝑗 ] } , 4 {(𝑚 𝑟𝑗𝑘 𝑟𝑗𝑘 4𝜋𝜇0 𝑟𝑗𝑘 𝑟𝑗𝑘 (18)

ACCEPTED MANUSCRIPT

where 𝜇0 is the magnetic permeability in a vacuum, 𝑚 ⃗⃗⃗𝑗 is the magnetic dipole moment of the j-th particle, ⃗⃗𝑗𝐻 is and 𝑟⃗𝑗𝑘 = 𝑟⃗𝑗 − 𝑟⃗𝑘 is the relative position vector. The magnetic torque due to the applied magnetic field 𝑇 given by: ⃗⃗𝑗𝐻 = 𝑚 ⃗⃗ , 𝑇 ⃗⃗⃗𝑗 × 𝐻 (19) 𝑀 ⃗⃗ is the applied magnetic field vector, and the magnetic interaction torque 𝑇 ⃗⃗𝑗𝑘 where 𝐻 is expressed by:

1 3 ⃗⃗⃗𝑗 × 𝑚 ⃗⃗⃗𝑘 − 2 (𝑚 ⃗⃗⃗𝑘 ⋅ 𝑟⃗𝑗𝑘 )𝑚 ⃗⃗⃗𝑗 × 𝑟⃗𝑗𝑘 ] . 3 [𝑚 𝑟𝑗𝑘 4𝜋𝜇0 𝑟𝑗𝑘

CR IP T

𝑀 ⃗⃗𝑗𝑘 𝑇 =−

(20)

The position of the particle is developed by solving the equations of motion, Eqs. (8) and (9).

AN US

3. Flow around a fixed spherical particle 3.1. Analytical model

To verify the simulation code and determine the appropriate setup for the immersed boundary method to describe a spherical particle in a flow, the drag force of a fixed spherical particle in a uniform flow is examined using the developed method. The size of the numerical region is 21𝐷 × 11𝐷 × 11𝐷, where D is the number of grids that correspond to the diameter of the fixed particle [32,33]. The particle, whose diameter is

M

described as five grids, is fixed at 5.5D from the inlet boundary on the center line of the numerical region. A uniform flow of constant velocity U is given as the velocity inlet boundary condition at the entrance of the numerical region. At the side boundaries, velocity inlet boundary conditions are also imposed and the given

ED

flow velocity has only x-component, of which the value is U. A free outflow condition is imposed at the exit of the numerical region. The local distribution functions are given at the boundary surface to satisfy the

PT

boundary condition. The force acting on the spherical particle is calculated and the drag coefficient is obtained. The Reynolds number based on the size of the particle 𝑅𝑒 = 𝐷𝑈/ν, where ν is the kinetic viscosity, is small, and thus the flow around a particle is a Stokes flow. Simulations are performed for several

CE

Reynolds numbers of 0.04, 0.1, 0.2, 0.5 and 1.0, and the results compared with theoretical results, such as the Stokes law and the Oseen approximation. The number of the Lagrangian boundary points on the surface of the spherical particle is 402 and these boundary points are located using the Uhlmann method based on

AC

the arrangement of equilateral triangle panels [28]. The relaxation time 𝜏 is set at 0.7 and the density of the fluid 𝜌, is 1 for lattice Boltzmann simulations. 3.2. Simulation result and drag coefficient The drag force acting on a fixed spherical particle is calculated using the developed simulation code. More than five grids corresponding to the diameter of the particle are sufficient to describe the behavior of a magnetic particle in the fluid, as Satoh et al. reported for their two dimensional lattice Boltzmann simulations [34]. Drag coefficients are calculated for diameters D, from 3 to 10, and the five grids for the diameter of a particle are determined as appropriate for simulation of the particle behavior in the fluid. Thus, five grids are assumed for the diameter of a magnetic particle and a suitable parameter for the immersed boundary method is determined from the simulation results. The parameters for immersed boundary method in Eq. (7) are determined as follows. The force acting on the fluid from a small spherical

ACCEPTED MANUSCRIPT

particle in a flow of low Reynolds number is expressed by the Stokes law, i.e., 𝐹⃗ 𝑣 = −3𝜋𝜌𝜈𝐷(𝑢 ⃗⃗ − 𝑑𝑟⃗/𝑑𝑡), where ν is the coefficient of kinetic viscosity, which is given as ν = (𝜏 − 1/2)/3 in this lattice Boltzmann method. This force corresponds to the total force from the boundary points on the particle. Therefore, we approximately derive 𝑁𝑏 𝜅 /𝜖𝑑 = 3𝜋𝜌𝜈𝐷 from Eq. (7). The stiffness factor can be determined from the simulation result of Re=0.04 and 𝜖𝑑 = 1.425 is obtained. From this result, the force scale factor can be determined as 𝜅 = 1.114 × 10−2 for 402 boundary points arranged on the surface of the particle. Figure 2 demonstrates the drag coefficients for validation of these parameter settings. Comparison of the simulation

AN US

CR IP T

results with the theoretical results indicates good agreement.

Fig. 2. Simulation results for the drag coefficient of a spherical body in a uniform steady flow. The open circle indicates the simulation results, and the solid line shows that using the Stokes law, 𝐶𝐷 = 24/𝑅𝑒, while the

M

dashed line indicates the result from the Oseen approximation, 𝐶𝐷 = 24(1 + 3𝑅𝑒/16)/𝑅𝑒.

ED

4. Behavior of magnetic particles in MR fluid shear flows 4.1. Analytical model for simple shear flows of MR fluids

PT

Figure 3 shows a schematic of the analytical model used to simulate the behavior of magnetic particles in a MR fluid between two parallel plates. The upper plate starts to move with constant velocity U in the x-direction and a uniform magnetic field is simultaneously applied in the z-direction, which is perpendicular

CE

to the wall surface. Periodic boundary condition are imposed in the x and y directions, while no-slip condition is imposed at the surfaces of both the lower and upper plates. All magnetic particles have the same diameter and the basic numerical region is 56 × 56 × 106 grids for lattice Boltzmann simulations, where the

AC

diameter of magnetic particles corresponds to five grids. The size of the basic numerical region is determined as the minimum appropriate size for the microstructure formation of magnetic particles after comparison the simulation result of this size with the results of larger and smaller sized basic numerical regions. The size was chosen as the minimum size that can be obtained qualitatively the same phenomenological structures. The number of Lagrangian boundary points on the surface of the spherical magnetic particle is 402, as in the previous simulations. The relaxation time 𝜏 is 0.7 and the density of the fluid is 1.0. Magnetic particles are arranged randomly in the initial state. The volume fractions of magnetic particles in the fluid are 0.08, 0.16 and 0.32. The Reynolds numbers, defined by 𝑅𝑒 = 𝑈𝐿𝑧 /ν, are 0, 0.004, 0.02, 0.04, 0.2, and 0.4, i.e., the constant velocity of the upper plate is changed. The magnetic flux densities of the applied uniform magnetic field normalized according to B0 are 0, 12.5 and 25.0, where B0 is the constant conversion factor. The conversion factor is the ratio of the magnetic flux density in the lattice Boltzmann method against that in the discrete particle method. In principle, it is difficult to determine the value of B0 because the magnetic

ACCEPTED MANUSCRIPT

interaction force is not directly substituted to the lattice Boltzmann equation, i.e., the external force, which is substituted to the lattice Boltzmann equation, is calculated using the immersed boundary method. To determine the conversion factor is the further subject. The magnetic particles are assumed to be iron particles with diameters of 1.2 μm and the base fluid is kerosene. The time is represented by the time ∆𝑥/𝑐, where ∆𝑥 is the grid size, and c is the particle velocity in the lattice Boltzmann method. The microstructure

AN US

CR IP T

formation and behavior of magnetic particles in the fluid is investigated.

Fig. 3. Analytical model for magnetic particles in an MR fluid between two parallel plates.

M

3.2. Behavior of magnetic particles in MR fluid shear flows

Figure 4 shows the time development of the microstructure formation of chain-like clusters in the MR fluid under a magnetic field when the upper plate is stationary. Immediately after application of the

ED

magnetic field, the magnetic particles begin to move and form chain-like clusters that are arranged in the direction of the magnetic field. After the time 2 × 104, the contact coefficient, as defined in Ref. [16], takes

PT

an almost constant value of 0.86, and the microstructure then becomes steady-state. After the time 105, almost all of the particles move to form the chain-like clusters, as shown in Fig. 4(d). First, small chain-like clusters are formed, which then combine with each other to form long chain-like clusters. These long

CE

chain-like clusters then combine with each other to form larger structures, i.e., fibrous structures, as shown in Fig. 4(d-1). This result is the same as that previously reported [16,17]. Figures 5 and 6 show the time development of the microstructure formation of magnetic particles in MR

AC

fluid shear flows. Figure 7 shows the time development of the number of clusters in the basic numerical region. When particles contact each other, they become members of the same cluster. From Fig. 7, the variation of the number of clusters becomes quite small after the time 6 × 104 , which indicates that the flow field becomes steady before this time. Magnetic particles interact magnetically with each other and form chain-like clusters in the direction of the magnetic field. However, the final microstructure is dependent on the Reynolds number. When the Reynolds number is small, chain-like clusters form fibrous structures and these structures move slowly with the fluid flow, as shown in Fig. 5. The shear of the flow is not sufficient to break chain-like clusters, i.e., the magnetic interactions between magnetic particles are stronger than the shear force. When the Reynolds number is small, the shear flow is weak, so that the chain-like clusters maintain their structure in the direction of the magnetic field because the magnetic interactions between magnetic particles are strong enough to maintain the connection and hold the chain-like clusters in the field direction. From side views, the axes of the chain-like clusters are almost in the field direction. In contrast,

ACCEPTED MANUSCRIPT

completely different structures appear when the Reynolds number is relatively high. From Fig. 6, a sheet-like structure in the direction of flow appears as the final microstructure. These snapshots confirm that the final microstructure changes with an increase of the Reynolds number, and a sheet-like structure appears more clearly with increased Reynolds number. In addition, the axes of the chain-like clusters are more inclined with an increase in the Reynolds number. The simulation results indicate that reconnection and rearrangement occur around the central area between two parallel plates. This phenomenon suggests that the chain-like clusters are repeatedly broken and rebuilt near the center due to the gap between the velocity of the particles near the upper plate and the velocity of the particles near the lower plate. Figure 8 illustrates the relationships between the velocity profiles of the flow and the velocities of the magnetic particles. In Fig. 8(a), magnetic particles move with the fluid flow in the absence of magnetic field,

CR IP T

and the flow calculated by simulation is identified as Couette flow, where the velocity is proportional to the z-position. In contrast, magnetic particles do not follow the fluid flow in the presence of magnetic field as shown in Figs. 8(b) and (c). Fibrous microstructures of magnetic particles move with the fluid flow and the microstructure is basically maintained in a weak shear flow under an applied magnetic field. However, the formed structure moves with rotation and deformation due to the applied shear and interactions between clusters, so that the velocity of magnetic particles is different. On the other hand, when the sheet-like

AN US

structure forms, most of the chain-like clusters that compose the sheet-like structure maintain their axes in the field direction; therefore, most of the particles have the same velocity. The velocity of the particles changes rapidly with respect to the velocity near the wall surface of the plate. When a magnetic field of relatively strong intensity is applied, magnetic particles located in the central area have almost the same velocity because magnetic interactions between magnetic particles are sufficiently strong to keep the chain-like clusters in the field direction. Chain-like clusters in the sheet-like structure are slightly inclined

M

due to the shear; however, the axes of the chain-like clusters still remain in the field direction. A fluctuation of the flow velocity is also noticeable due to the magnetic particle motions. Figure 9 shows the effect of the Reynolds number on the microstructure formation of magnetic particles

ED

when a weak magnetic field is applied. A sheet-like structure does not appear clearly in the presence of a weak magnetic field, even when the Reynolds number is 0.4. In addition, magnetic interactions between

PT

magnetic particles become weak, so that the inclination angle of axes of chain-like clusters becomes larger in relation to the field direction with an increase of the Reynolds number. These results indicate that the appearance of the sheet-like structure requires two conditions: a strong applied magnetic field and a

CE

relatively high Reynolds number. Figure 10 shows the histograms of magnetic particle distributions. From Fig. 10, the particle are almost uniformly distributed in the fluid without a flow, while the anisotropic distributions appear in shear flows. Aggregation of magnetic particles in the limited area appears clearly as

AC

shown in Fig. 10(d-2), and these aggregation correspond to the sheet-like structures. Figure 11 shows the average velocity profile of the flow and the effect of the applied magnetic field intensity on the flow velocity. The average velocities are average of the velocities between the time 2 × 104 from the time 8 × 104. The flow is accelerated near the lower plate, while the flow is deaccelerated near the upper plate. This is because the fluid flow is affected by the motion of both magnetic particles and the plates due to the fluid viscosity. The axes of the chain-like clusters remain in the field direction: therefore, the fluid near the upper plate is pushed back by the clusters, while the fluid near the lower plate is pushed forward by the clusters. Figures 12 and 13 show the effect of the volume fraction of magnetic particles on the microstructure formation in the presence of a magnetic field. Formation of the sheet-like structure is observed when the Reynolds number is 0.4, while the fibrous structures are formed when the Reynolds number is 0.04. When the volume fraction of magnetic particles is 0.32, the sheet-like structure does not appear clearly; however,

ACCEPTED MANUSCRIPT

most of the chain-like clusters are rearranged in the flow direction. This result indicates that the volume fraction of magnetic particles is not a requirement for the formation of a sheet-like structure. Figure 14 shows the average velocity profile of the flow along the central line of the basic numerical region, which indicates the influence of the volume fraction of magnetic particles on the flow velocity; the gradient of flow velocity around the central area between the two parallel plates is small, while the gradient of flow velocity near the wall is large. We have discussed the microstructure formation of magnetic particles in the presence of a magnetic field qualitatively. In further studies, quantitative analyses are expected to be realized through further

PT

ED

M

AN US

CR IP T

development of the hybrid simulation code.

CE

Fig. 4. Snapshots of magnetic particle distributions in MR fluid. The Reynolds number is 0, the magnetic flux density is 25.0, and the volume fraction of magnetic particles is 0.16. The times are (a) 0 (initial state), (b) 1.0 × 104, and (c) 1.0 × 105 (steady state). The snapshots show views from the (1) z-direction (top), (2)

AC

y-direction, and (3) x-direction.

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 5. Snapshots of magnetic particle distributions in MR fluid shear flow. The Reynolds number is 0.04, the applied magnetic flux density is 25.0, and the volume fraction of magnetic particles is 0.16. The times are (a)

ED

1.0 × 104, (b) 4.0 × 104, and (c) 1.0 × 105 (steady state). The snapshots show views from the (1) z-direction

AC

CE

PT

(top), (2) y-direction, and (3) x-direction.

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 6. Snapshots of magnetic particle distributions in MR fluid shear flow. The Reynolds number is 0.4, the applied magnetic flux density is 25.0, and the volume fraction of magnetic particles is 0.16. The times are (a)

ED

1.0 × 104, (b) 4.0 × 104, and (c) 1.0 × 105 (steady state). The snapshots show views from the (1) z-direction

AC

CE

PT

(top), (2) y-direction, and (3) x-direction.

Fig. 7. Time development of the number of clusters in the basic numerical region. The volume fraction of magnetic particles in the fluid is 0.16, and the applied magnetic flux density is 25.0.

(b)

AC

CE

PT

ED

M

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(c) Fig. 8. Velocity profile of the flow and x-component of the velocity of magnetic particles in the flow direction (a) in the absence of a magnetic field, and (b,c) in the presence of a magnetic field. The solid line indicates the flow velocity of Couette flow (theory), the dashed line shows the flow velocity profile along the central line of the numerical region, and the open circles represent the velocity of magnetic particles. The Reynolds numbers are (a) 0, (b) 0.04, and (c) 0.4. The time is 1.0 × 105 (steady state), and the magnetic flux density is 25.0.

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 9. Snapshots of magnetic particle distributions in MR fluid shear flow. Steady state and with a Reynolds

M

number of (a) 0.04, (b) 0.2, and (c) 0.4. The magnetic flux density is 12.5, the volume fraction of magnetic particles is 0.16, and the steady state. The snapshots show views from the (1) z-direction (top), (2) y-direction,

AC

CE

PT

ED

and (3) x-direction.

ACCEPTED MANUSCRIPT

Fig. 10. Histograms of magnetic particle distributions. (a) shows the initial state and the others show the steady states. The Reynolds numbers are (b) 0, (c,e) 0.04, and (d,f) 0.4. The volume fraction of magnetic

ED

M

AN US

(a)

CR IP T

particles is 0.16 and the applied magnetic flux densities are (b,c,d) 25.0 and (e,f) 12.5.

(b)

Fig. 11. Average velocity profile of flow along the central line between two parallel plates. The volume

AC

CE

PT

fraction of magnetic particles is 0.16, and the Reynolds numbers are (a) 0.04 and (b) 0.4.

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 12. Snapshots and histograms of magnetic particle distributions in MR fluid shear flow. The Reynolds number is 0.04, the magnetic flux density is 25.0, and the steady state. The snapshots show views from the

M

(1) z-direction (top), (2) y-direction, and (3) x-direction, and the histograms show along the (1) z-direction (top), (2) y-direction, and (3) x-direction. The volume fractions of magnetic particles are (a,d) 0.08, (b,e) 0.16,

AC

CE

PT

ED

and (c,f) 0.32.

ACCEPTED MANUSCRIPT

Fig. 13. Snapshots of magnetic particle distributions in MR fluid shear flow. The Reynolds number is 0.4, the magnetic flux density is 25.0, and the steady state. The snapshots show views from and the histograms shows along the (1) z-direction (top), (2) y-direction, and (3) x-direction. The volume fractions of magnetic

ED

M

AN US

(a)

CR IP T

particles are (a,d) 0.08, (b,e) 0.16, and (c,f) 0.32.

(b)

PT

Fig. 14. Average velocity profile for flow along the central line between two parallel plates. The applied magnetic flux density is 25.0, and the Reynolds numbers are (a) 0.04, and (b) 0.4. The key legend indicates

CE

the volume fractions of the magnetic particles.

4. Concluding remarks

AC

A new hybrid simulation code that combines three simulation methods, the lattice Boltzmann method, the immersed boundary method and the discrete element method, and which also considers magnetic interactions between magnetic particles, has been developed to calculate behavior of magnetic particles in MR fluid flows. The behavior and microstructure formation of magnetic particles in simple shear flows of MR fluid between two parallel plates have been simulated. A sheet-like structure appears in flows with a relatively high Reynolds number under strong applied magnetic field. In the presence of weak magnetic field, fibrous structures of chain-like clusters of magnetic particles move in the fluid flows.

References [1] Rosensweig RE. Ferrohydrodynamics, Cambridge University Press 1985; 7-8. [2] Kciuk M, Turczyn R. Properties and application of magnetorheological fluids. J Achievements in

ACCEPTED MANUSCRIPT

Materials and Manufacturing Engineering 2006;18:127-130.

[3] Baranwal D, Deshmukh TS. MR-Fluid Technology and Its Application –A Review. Int J Emerging Tech Advanced Engineering 2012;2:563-569. [4] Spaggiari A. Properties and applications of Magnetorheological fluids. The Italian research on smart materials and MEMS 2013;23:57-61. [5] Oyama T, Kitakami O, Fujita T, Shimada Y, Nishiyama H. Cluster Structure and Rheology of Magnetorheological Fluids. JSME annual meeting 2000;4:273-274 (in Japanese) [6] Dominguez-Garcia P, Melle S, Pastor JM, Rubio MA. Scaling in the aggregation dynamics of a magnetorheological fluid. Phys. Rev. E 2007;76:051403. [7] Yongzhi L, Xinhua L, Hao L. Dynamic Simulation for Three-Dimensional Micro-Structure of

CR IP T

Magnetorheological Fluids. Advances in information Sciences and Service Sciences 2012;4:55-61. [8] Yao GZ, Yap FF, Chen G, Li WH, Yeo SH. MR damper and its application for semi-active control of vehicle suspension system. Mechatronics 2002;12:963-973.

[9] Sapinski B. Magnetorheological dampers in vibration control of mechanical structures. Mehcanics 2009;28:18-25.

[10] Chen C, Liao WH. A self-sensing magnetorheological damper with power generation. Smart Material

AN US

Structure 2012;32:0251014.

[11] Hong KP, Cho YK, Shin BC, Cho MW, Choi SB, Cho WS, Jae JJ. Magnetorheological (MR) Polishing of Alumina-Reinforced Zirconia Ceramics Using Diamond Abrasives for Dental Application. Materials and Manufacturing Processes 27;2012:1135-1138.

[12] Zhu Y, Haddadian E, Mou T, Gross M, Liu J. Role of nucleation in the structure evolution of a magnetorheological fluid. Phys Rev E 53;1996;1753-1759.

M

[13] Mohbi M, Jamasbi N. Simulation of the formation of nonequilibrium structures in magnetorheological fluids subject to an external magnetic field. Phys Rev E 54;1996:5407-5413. [14] Gullery GL, Tao R. Structures of a Magnetorheological Fluid. Int J Modern Phys B 2001;15:851-858.

ED

[15] Ly HV, ITO K, Banks HT, Jolly MR, Reitich F. Dynamic simulation of the temporal response of microstructure formation in magnetorheological fluids. Int J Modern Phys B 2001;15:894-903.

PT

[16] Ido Y, Inagaki T, Umehara N. Numerical simulations of distributions of magnetic and nonmagnetic particles in MAGIC fluids. Magnetohydrodynamics 2008;44:83-91. [17] Ido Y, Inagaki T, Yamaguchi T. Numerical simulation of microstructure formation of suspended particles

CE

in MR fluids. J Phys Condensed Matter 2010;22:324103. [18] Sherman SG, Paley DA, Wereley NM. Massive parallel simulations of chain formation and restructuring dynamic in a magnetorheological fluid. Proc ASME 2011 Conf Smart Materials, Adaptive

AC

Structures and Intelligent Systems 2011:SMASIS2011-5188. [19] Li H, Peng X. Simulation for the Microstructure and Rheology in Bidisperse Magnetorheological Fluids. J Comput Phys 2012;7:1405-1412. [20] Zhang H, Tan Y, Shu S, Niu X, Trias FX, Yang D, Li H, Sheng Y. Numerical investigation on the role of discrete element method in combined LBM-IBM-DEM modeling. Computers & Fluids 2014;94:37-48. [21] McNamara GR, Zanetti G. Use of the Boltzmann Equation to Simulate Lattice-Gas Automata. Phys Rev Lett 1988;61:2332-2335. [22] Chen S, Doolen G. Lattice Boltzmann method for fluid flows. Annu Rev Fluid Mech 1998;30:329-364. [23] Peskin CS. Immersed boundary method. Acta Numerica 2002;11:1-39. [24] Mittal R, Iaccarino G. Immersed boundary methods. Annu Rev Fluid Mech 2005;37:239-261. [25] Luding S. Introduction to Discrete Element Methods. European J Environmental and Civil Engineering 2008;12:785-826.

ACCEPTED MANUSCRIPT

[26] Chen S, Chen D, Mattaeus WH. Lattice Boltzmann model for simulation of magnetohydrodynamics. Phys Rev Lett 1991;67:3776-3779. [27] Feng ZG, Michaelides EE. The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems. J Comp Phys 2004;195:602-628. [28] Uhlmann M. An immersed boundary method with direct forcing for the simulation of particulate flows. J Comp Phys 2005;209:448-476. [29] Klingenberg DJ, von Swol F, Zukoski CF. Dynamics simulation of electrorheological suspensions. J Chem Phys 1989;91:7888-7895. [30] Enomoto Y, Oba K, Okada M. Simulation study on microstructure formations in magnetic fluids. Physica A 2003;330:496-506.

CR IP T

[31] Niu XD, Shu C, Chew YT, Peng Y. A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows. Phys Lett A 2006;354:173-182.

[32] Shi Y, Sader JE. Lattice Boltzmann method for oscillatory Stokes flow with applications to micro- and nanodevices. Phys Rev E 2010;81:036706.

[33] Wu J, Shu C. An improved immersed boundary-lattice Boltzmann method for simulating three-dimensional incompressible flows. J Comput Phys 2010;229:5022-5042.

AN US

[34] Satoh A, Chantrell RW. Basic study on the lattice Boltzmann method for application to many particle dispersions: a uniform flow past a two-dimensional circular particle. Proceedings of the JSME Fluids

AC

CE

PT

ED

M

Engineering Conference 2009:555-556 (in Japanese).