Numerical analysis for the normal force of magnetorheological fluids

As one of the most important mechanical parameters of magnetorheological fluids (MRF), the normal force has been studied by a large number of researchers using various experimental equipment and methods. However, it is difficult to reveal the effects of particle microstructure characteristics on the normal force of MRF by experimental methods, especially for such factor as the particle size distribution. To advance this research area, a numerical method for calculating the normal force of MRF is proposed in this study. The accuracy of the simulation model is verified with the experimental results measured by a plate-plate magnetorheometer. The results show that the simulated values agree well with the experimental data, which indicates the feasibility of calculating the normal force of MRF using numerical methods and provides a new research idea for a more intuitive and comprehensive analysis of the normal force characteristics of MRF. Besides, the simulated data show that the increase in magnetic field intensity and particle volume fraction will sufficiently enhance the normal force of MRF and improve the response time of MRF. For the MRF with the same particle volume fraction, a decline in the average particle diameter will increase the normal force. Moreover, increasing the dispersion of particle size of MRF particles can also improve its normal force.


Introduction
Magnetorheological fluids (MRF) is a type of intelligent fluid, and its rheological properties can be controlled by an external magnetic field. It is mainly composed of nonmagnetic carrier liquid, magnetic particles and some additives used to improve the performance of MRF [1,2]. In the natural state, MRF exhibits the general characteristics of the Newtonian fluid. As it is subjected to an external magnetic field, the suspended magnetic particles in the carrier liquid will form chain, column and network structures within a few milliseconds [3]. The presence of these microstructures drastically changes the rheological properties of MRF [4]. First, the shear yield stress will be increased by several orders of magnitude, which has been well used in the field of vibration control. Numerous energy absorption and damping devices have been developed, such as magnetorheological (MR) vibration absorbers and MR dampers [5][6][7]. Second, the formation of particle microstructures also changes its mechanical properties toward the magnetic field. Due to the squeezing effect of free particles on the particle chains, a normal force along the direction of the magnetic field is generated [8]. Recently, developments and applications of normal force of MRF have also obtained satisfying results in the field of squeeze mode MR dampers and MR valves [9,10].
The current studies on the normal force of MRF are basically based on experimental methods, generally using a plate-plate rheometer to test the normal force of MRF under different conditions. These studies dealt with the effects of shear rate, strain amplitude, temperature, particle concentration, MRF layer thickness and other factors on the dynamic and static normal force of MRF [11][12][13][14][15][16][17]. Simultaneously, some theoretical models were also developed, and the results showed that the theoretical predictions of the MRF normal force agreed well with the experimental results [8,13]. Besides, some researchers have studied the normal stresses of special magnetorheological materials such as magnetorheological gels, magnetorheological foams and magnetorheological greases, and tried to give explanations from a microscopic point of view corresponding to the experimental results [18][19][20]. Recently, Li used a combination of theoretical calculations and experimental tests to preliminarily analyze the effects of the normal force on the magnetorheological finishing [21].
However, these aforementioned studies can only analyze the normal stress characteristics of MRF at the macroscopic level, but cannot correlate the essential causes of these mechanical characteristics to changes in the microstructure of particles. This limits the deeper understanding of the normal stresses characteritics of MRF. In contrast, the molecular dynamics simulation is very effective for solving such problems. In recent years, the analysis of the shear motion patterns of MRF using numerical computational methods has received a lot of attention from researchers [22][23][24]. However, up to now, there has been little research on the analysis of the normal stress of MRF using molecular dynamics simulation methods.
Therefore, a numerical method is proposed to analyze the normal force of the MRF with 10%, 15% and 20% particle volume fraction in this study. Firstly, according to Newton's second law, the kinetic equations of the magnetic particles suspended in MRF and the contact model between the wall particles and suspended particles are established. The time step of 0.1 μs is used to simulate the particle movements in the entire simulation area, and the total normal force of MRF is obtained by accumulating the normal forces of all wall particles. Secondly, the normal force of the MRF with the same particle volume fraction is tested using a plate-plate magnetorheometer. Then the experimental data is used to compare with the numerical result. Finally, the proposed numerical method is used to analyze the effect of particle size and particle size distribution on the normal force of the MRF.

Materials
The non-magnetic carrier liquid used in the simulation model and experiments is the silicone oil with a viscosity of 0.1 Pa · s and density of 975 kg · m 3 . The suspended particles in MRF are carbonyl iron powders (CIP) with a density of 7860 kg · m 3 and an average particle size of 5 μm. As shown in figure 1, the relative permeability μ r of CIP is measured by a vibrating sample magnetometer (VSM). figure 2 illustrates the relationship between the relative permeability of CIP and the magnetic flux density. Besides, we introduced the monitoring wall particles in the simulation process, these particles have no actual physical significance, and are presented only as sensors for measuring the normal force of MRF.

Simulation technique
As shown in figure 3, the initial distribution area of MRF is limited to a cube with a side length of 150 μm, and the size of the simulation region is determined after considering the authenticity and efficiency of the simulation algorithm. We should ensure that the simulation results can sufficiently reflect the kinetic properties of each magnetic particle and entire flowing deformation of MRF. Simultaneously, it should be avoided that the computation time is too long due to too many particles. Generally speaking, when the side length of the simulation cube reaches fifteen times of the average particle size, the dimensions can meet the above requirements [22]. The four boundaries of the x and y directions are set to be the periodic boundary conditions to ensure that the particle volume fraction in the simulation region remains unchanged. The upper and lower surfaces in the z-axis direction are composed of wall particles, which are represented by tightly arranged soft spherical particles with a diameter of 2 μm. The positions of all wall particles remain constant during the simulation. These wall particles serve to limit the range of motion of the magnetic particles while sensing the normal force generated by the MRF.
The influences of gravity and Brownian force are usually neglected in the simulation procedure because these two forces are several orders of magnitude smaller than the magnetic forces between particles and the drag forces of the carrier liquid on particles. Besides, the simulation temperature is assumed to be constant in the simulation. Therefore, we should consider the interaction forces between the particles, and the forces between the particles and carrier fluid. The kinetic equations of magnetic particles are proposed based on the Newton's second law.ˆ( where m i is the mass of magnetic particle i, r i is the particle displacement, F i m is the magnetic force on particle i. n F i is the drag force of the carrier liquid on particle i, η = 0.1 Pa · s is the viscosity of the continuous phase, F i r is the repulsive force generated when particle i contacts another particle. F i m is expressed as  here, μ 0 = 4π × 10 −7 T/m · A is the vacuum permeability. χ = μ r − 1 is the magnetization rate of the magnetic particles, which can be calculated from the data shown in figure 2. H 0 is the external magnetic field applied on the MRF, with corresponding magnetic flux density ranging from 0.1 T to 0.5 T, and the magnetic field is parallel to the z-axis direction in the simulation region. R i and R j are the radii of particles i and j, r ij is the distance between the particles i and j, and θ is the angle between the magnetic field direction and line that connects the centercenter points of the particles i and j.r ij is the unit vector of r ij and k is the unit vector of magnetic field. For equation (1), F 0 is derived from the equilibrium state of the magnetic force and repulsion force when two particles are just in contact. The formula is 0 are the magnetic moments of particles i and j. The contact force F i w between the wall particles and magnetic particle i is described by the soft sphere model as where, k w is the adjustable contact parameter, which is equal to 50 in our model. s is the contact area of two particles, and e n is the normal vector of the two particle contact surfaces.
Since equation (1) is difficult to obtain an analytical solution, the Verlet-Velocity algorithm, commonly used in molecular dynamics, is introduced to decompose the particle motion. The specific iterative formulas are, see [25] According to the above equations, the spatial position of a magnetic particle in MRF at the moment t + δt can be calculated by giving the initial particle position ( ) t r , velocity ( ) t v , acceleration ( ) t a and the time step δt. In this study, the simulation time step is set to 0.1 μs. Furthermore, the force F i m (see equation (2)) and acceleration ( ) d + t t a on the particle at the moment t + δt can be obtained. See equation (6), the particle velocity ( ) at the moment t + δt is obtained before entering the next round of the iteration process.

Results and discussion
4.1. Uniform particle size To obtain the response time of MRF to generate stable normal force, the magnetic energy of the particle system is calculated firstly. The magnetic energy of the magnetic particle i can be expressed as where B ji is the magnetic field applied on the particle i. According to equation (7), the magnetic energy of the whole particle system is proposed as When the particle diameter is 5 μm, magnetic flux density is 0.1 T and particle volume fraction is 10%, as shown in figure 4, the energy-simulation time curve can be obtained from the simulation. The results show that the magnetic energy of the particle system reduces quickly with the increase of simulation time. The change of the system magnetic energy is inconspicuous when the simulation time reaches 2 ms, which indicates that the particle system has reached a stable state. Then, the simulation ends when the time step is superimposed to this value. Additionally, the simulation cutoff time for all other simulation cases in this study was chosen to be 2 ms, which was proposed by comparing a large number of simulation results. It is found that the particle size, volume fraction and the increase of the excitation magnetic field will reduce the response time for the particle microstructures of MRF to reach stability. Furthermore, the MRF with an average particle size of 3 μm and a volume fraction of 10% has the longest response time in all the simulation models, but it is also less than 2 ms. Therefore, using the cutoff time of 2 ms can ensure that the normal force of the MRF obtained in the study is reliable. Figures 5, 6 and 7 show the normal force response of MRF in magnetic fields of different intensities for volume fractions of 10%, 15% and 20%, respectively. The reason for selecting these three particle volume fraction is mainly because the random distribution conditions of the initial positions of particles limit the position initialization process of high concentration particles. The maximum concentration of the MRF that is easy to fill with randomly distributed particles is about 20%. Although the concentration of particles can be increased by some regular distribution methods, the distribution of particles also loses randomness. Comparing the curves in this three Figures, it can be seen that the normal force of MRF will increase significantly with the increase of magnetic flux density and volume fraction, and the relationship between them is nonlinear. Additionally, the response time of MRF will also decrease due to the enhancement of the magnetic field.
As shown in figure 8, to verify the accuracy of the simulation results shown in figures 5 to 7, the static normal forces of MRF with volume fractions of 10%, 15% and 20% and an average particle size of 5 μm are tested using an Anton Paar MCR-302 plate-plate magnetorheometer. The height of the cylindrical MRF layer between the two test plates is 0.15 mm, and the effective radius is 10 mm. For comparison, the numerical results are  . It can be seen from the figure 9 that the numerical results basically coincide with the experimental data, which indicates that the numerical method proposed in this study is effective for calculating the normal force of MRF in a static magnetic field. However, there are still some errors. It is possible that these errors arise from the experiment itself because the MRF layer is very thin, and even though we used a syringe to add the MRF sample, some small volume changes and differences in the sample addition position can affect the experimental results. We had to perform many replicate experiments to reduce the experimental error. Besides, we used a simplified dipole model in the numerical algorithm and assumed that the particle sizes are the same, which could cause some errors in the simulation results.

Polydisperse particle size
It is well known that the particle size of the CIP is not uniform, but shows a certain distribution pattern. It can be known from some previous studies that the distribution of particle size has a certain influence on the properties of MRF [22,24]. However, the effects of particle size distribution on the normal force of MRF are difficult to investigate experimentally, hence this problem will be explored for this section by using the numerical method.  Normally, the Gaussian distribution function can be used to describe the particle size distribution as where a 0 = 5 μm is the average particle diameter. σ 0 is the standard deviation, whose values are are set to be 0.01, 0.1, or 0.2 in equation (9). When the particle volume fraction of MRF is set to be 10%, according to equation (9), the particle numbers in each diameter range can be obtained as shown in figure 10. figure 11 shows the microstructures of MRF with different particle size distributions when the magnetic flux density is 0.1 T and simulation time is 2 ms. The values of σ 0 for figures 11(a), (b) and (c) are 0.01, 0.1, and 0.2, respectively. According to the previous study, it can be known that the polydisperse particle size will make more thick and twisted particle chains formed in MRF [26], which leads to an increase in the external magnetic moment output of the particle chains and thus causes an enhancement of the normal force of MRF. As shown in figure 12, it can be seen that as the standard deviation σ 0 rises from 0.01 to 0.2, a small increase in the normal force of MRF occurs. However, since the particle volume fraction is relatively small (f = 0.1) for these simulation models, there is enough space in the carrier fluids for the magnetic particles to form straight singlechain structures, hence the effect of the particle size distribution pattern on the normal force is still relatively limited. On the other side, the influence of the average particle size a 0 on the normal force of the MRF is much     Figure 12. The normal forces of MRF with different particle size distributions when the particle volume fraction is 10%. more significant compared to the particle size distribution pattern. It can be seen from figure 12 that the normal force of the MRF will increases significantly as the average particle size a 0 decreases when the distribution parameter σ 0 and particle volume fraction of MRF are constant. This is due to the fact that there are more particle chains formed in the MRF.
Increasing the particle volume fraction of the MRF to 15% and 20% under the same simulation conditions, as shown in figures 13 and 14, it can be found that the effects of the variation in the particle size distribution (parameter σ 0 ) on the normal force of MRF becomes more pronounced. There are two main reasons for this phenomenon, the first is that the increase of particle volume fraction enhances the normal force generated by MRF, which makes the effects of the particle size distribution on the normal force to be amplified synchronously. Secondly, the increase of particle concentration makes it easier to form complex particle microstructures in MRF, and the polydisperse particle size facilitates the generation of the bending particle chains and squeezing forces between particles. Figure 13. The normal forces of MRF with different particle size distributions when the particle volume fraction is 15%. Figure 14. The normal forces of MRF with different particle size distributions when the particle volume fraction is 20%.

Conclusions
The study proposes a numerical method that can be used to compute the normal force of MRF, and the accuracy of the numerical model is verified by conducting the corresponding experiments using a plate-plate rheometer. The results show that the rise of external magnetic field strength and particle volume fraction will considerably enhance the normal force of MRF. Meanwhile, a stronger magnetic field will improve the response time of MRF. For the MRF with the same volume fraction, a decline in the average particle diameter will increases the normal force. Additionally, the normal force generated by MRF with polydisperse distributed particle size is greater than that of MRF with uniformly distributed particle size.