Microstructure Simulation and Constitutive Modelling of Magnetorheological Fluids Based on the Hexagonal Close-packed Structure.

This paper presents a new constitutive model of high particles concentrated magnetorheological fluids (MRFs) that is based on the hexagonal close-packed structure, which can reflect the micro-structures of the particles under the magnetic field. Firstly, the particle dynamic simulations for the forces sustained by carbonyl iron powder (CIP) particles of MRFs are performed in order to investigate the particles chain-forming process at different time nodes. Subsequently, according to the force analyses, a hexagonal close-packed structure, which differs from the existing single-chain structure and body-cantered cubic structure, is adopted to formulate a constitutive model of MRFs with high concentration of the magnetic-responsive particles. Several experiments are performed while considering crucial factors that influence on the chain-forming mechanism and, hence, change the field-dependent shear yield stress in order to validate the proposed model. These factors include the magnetic induction intensity, volume fraction and radius of CIP particles, and surfactant coating thickness. It is shown that the proposed modeling approach can predict the field-dependent shear yield stress much better than the single-chain model. In addition, it is identified that the shear yield stress is increased as the particle volume fraction increases and surfactant coating thickness decreases. It is believed that the proposed constitutive model can be effectively used to estimate the field-dependent shear yield stress of MRFs with a high concentration of iron particles.


Introduction
Magnetorheological fluids (MRFs), a new type of intelligent material, is usually composed of base carrier fluid and magnetic particles that are dispersed in the base carrier fluid. The flow state of MRFs can be rapidly governed by the imposed magnetic field. Some advantages of MRFs, such as easy controllability, quick response, and reversible change, enable it to display its huge potential applications in diverse fields including automotive suspension system [1], self-powered magneticfield sensor [2], and vibration control structures [3,4]. It is well known that the macroscopic properties of materials are determined by their microstructures. Thus, the changes in the microstructure of MRFs will directly affect its macroscopic performance. A large amount of research works has been reported regarding the performance and constitutive model of MRFs since it has been discovered in 1948 by Rabinow. The flow state of MRFs between two parallel fixed plates under the influence of a constant magnetic field perpendicular to the flow direction, and a constitutive model of incompressible MRFs were investigated in [5]. Jolly et al. have established a standard isolated-chain model that is feasible only when the direction of magnetic dipolar interaction of adjacent particles is along the applied magnetic field in view of the system consisting of infinite chains of particles in MRFs [6]. However, the volume fraction of magnetic particles in MRFs must be low enough to ensure the calculation accuracy, leading to the limitation of this model. Varela-Jiménez et al. proposed a semi-empirical constitutive model in view of the flow state transition of MRFs under the magnetic field in order to investigate MRFs with high volume fraction of magnetic particles [7]. It has been shown that this model properly reflects the relationship between the shear yield stress and magnetic induction strength, volume fraction of carbonyl iron powder (CIP) and other factors. Pan Sheng et al. deeply investigated the relationship between the shear yield stress and the strength of magnetic field, temperature, and composition of MRFs and then constructed a new shear stress model [8,9]. What is more, an optimum percentage of particles in the carrier fluid is obtained to improve efficiency [10]. As early as 1997, R. Tao et al. have mentioned that the arrangement of particles in MRFs might form a dense hexagon when the electric field and magnetic field are compatible [11].
On the other hand, several experimental works have been done to find the relationship between the shear yield stress and temperature at different shear rates under inhomogeneous magnetic field. Gao Chunfu addressed that most of the constitutive models of MRFs proposed so far have been formulated based on the theory of magnetic dipole moment, according to which the magnetic particles under the external magnetic field are regarded as magnetic dipoles. The classical theoretical models of magnetic dipole only usually take the influence of the adjacent magnetic chains on the target particle into account and can only analyze a single particle separately to obtain its magnetic dipole moment in the particle chain [12]. Taking the magnetic dipole theory as a foundation, Chen et al. calculated the interaction force between CIP particles, and then established the dipole boundary theoretical model that has been verified by experiments [13]. Ma et al. put forward a distance weighting factor model of the chain-structure utilizing the ampere molecular current hypothesis [14]. This model describes the interaction force between the magnetic particles, and the Monte Carlo method is adopted to simulate the microstructure of magnetic particles with the aim of exploring the influences of microstructure and macro-mechanics behaviors of materials. The simulating coordinates and magnetic moments of CIP particles are extracted in order to solve the shear stress value in theory. Enhanced magnetorheological performance of highly uniform magnetic carbon nano particles was explored by Seungae Lee et al. to investigate the influence of particle size on MR properties [15]. There is still a blank in the study of constitutive model for the chain structures of MRFs in high concentration considering the influence of magnetic particle chains on each other, as evident from the above literature survey from both theoretical and empirical aspects.
Consequently, the main technical contribution of this work is to propose a new constitutive model that is based on the hexagonal close-packed structure that can reflect the micro-structures of the particles in high concentrated MRFs. In order to achieve this goal, firstly the particle dynamics analyses of the forces sustained by CIP particles in MRFs are investigated, showing the motion model of CIP particles. Subsequently, the simulations of the chain-forming process of MRFs at different time nodes under imposed magnetic field are provided. According to the above force analyses, a hexagonal close-packed structure, which differs from the existing single-chain structure and bodycentered cubic structure, is proposed, which takes the effects of adjacent chains under the aggregation of carbonyl iron particles into consideration [16]. Subsequently, the constitutive model of MRFs is established by the particle dynamics method that is based on the hexagonal close-packed structure. A series of experiments on the influencing factors, such as the magnetic induction intensity, volume fraction and radius of CIP particles, and surfactant coating thickness, are carried out in order to validate the proposed constitutive model. From the comparative work between the predicted and measured results, it is shown that the proposed modeling approach can be effectively applied to understand the chain structures of the magnetic particles and, hence, provides a guideline to develop high performance of MRFs.

Methods
It is necessary to analyze the forces on CIP particles and its influence factors, and some equations of motion need to be established, to accurately study the microstructure of MRFs. After performing the analyses of the forces and the motion state, the microstructure of MRFs under the magnetic field can be simulated and an appropriate microstructure model is provided.

Particle Dynamics Analysis
The complex microstructure of MRFs and influence factors on the evolution process are difficult to be observed by experiments. Numerical simulation technology can overcome these shortcomings. As a step to completely develop the numerical simulation, this work deals with particle dynamics investigation on the microstructure of MRFs after comprehensive considerations [17][18][19]. The particle dynamics method (PDM) is a new type of simulation that was developed from molecular dynamics combined with discrete element method, of which the essence is to analyze the forces and motion state of particles based on the Newton's law. Meanwhile, the positions of all particles in the system are calculated and updated continuously with the changing time, and, at last, the stable microstructure of the particles system is obtained. The magnetized CIP particles are affected not only by magnetic force, but also by Brownian motion, repulsion force between CIP particles, gravity, buoyancy force of the base carrier fluid, and viscous resistance. The Brownian motion is overcome by CIP particles to form chains due to the large applied magnetic field, so Brownian motion can be neglected when the particle dynamics analysis is conducted. The gravity of CIP particles and the buoyancy force of base carrier fluid are also ignored because they are considerably smaller to the magnetic interaction force between magnetic dipoles under the action of magnetic field. Therefore, only magnetic force, viscous resistance, and repulsive force are considered in this work. There are two assumptions when the PDM is used; the magnetic particles are ideal sphere-shaped with the same size and the external imposed magnetic field is uniformly distributed.

Magnetic Field Force
To a certain extent, the magnetic field formed by magnetic dipole i is affected by that formed by other magnetic dipoles, which is always ignored by some of the proposed approaches. The magnetic interaction force Fi m of CIP particles is given by [20]: Where, a is the radius of CIP particles, H the magnetic field intensity, rij the center distance of the CIP particle i and j, which is not less than 2a,  the deviation angle between the central line of particle i and j and the direction of the magnetic field,  is the vacuum susceptibility with a value of 4 × 10 −7 ,  the relative magnetic susceptibility of particles, ̅ the unit position vector of the particle i and j, and ̅ three-dimensional the unit vector in the magnetic field direction, respectively.

Repulsive Force
When the CIP particles move under the action of magnetic field, a collision generating repulsive force occurs between the particles. The repulsive force, which can refrain the particles from position overlapping, increases exponentially with increasing displacement of CIP particles during the collision. The concept of material parameter is introduced for quantifying the rapidity of this increase of repulsion force. The simulation focuses on the iteration of the particle position, which is not subject to the motion state of particle itself. Thus, the moment that is generated by the particle rotation is neglected in the analysis in order to simplify the calculation. A relatively concise expression for calculating the repulsion force Fi r is given by [20]: Where,  denotes the material parameter and m is the mass of the CIP particle. The direction of repulsion force is along the line of two colliding particles centers. Suppose that the CIP particles are rigid spheres without elastic-plastic deformation. The parameter A is given for expressing zero interaction force, when two colliding particles are just in contact. At the same time, it is necessary to consider the collision between the vessel wall and the particles. The wall repulsion force F r i(wall) along the magnetic field direction is given by [21]: where, (4) and zi represents the ordinate of particle i in the simulant space, Lz is the coordinate length of Z axis in the simulant space. The equation of wall repulsion force perpendicular to the magnetic field is derived from Equation (3) by transforming the longitudinal coordinates of particles into horizontal coordinates and the longitudinal length into horizontal length, respectively.

Viscous Resistance
The base carrier fluid of MRFs is deionized water and the viscous resistance Fi d of CIP particles from deionized water is generally described by the Stokes equation [22] when moving in the base carrier fluid without shear: Where  is the motion velocity of CIP particles and  is the zero magnetic field viscosity of base carrier fluid. Under the magnetic field, MRFs undergoes the shear with a rate of . The shearing action is perpendicular to the direction Z of the applied magnetic field along the X direction. Subsequently, the viscous resistance Fi V of CIP particles from the base carrier is achieved from: In the above equation, ̅ represents the unit vector in the X direction of the simulant space.

Establishment of Motion Model
On the base of aforesaid analyses, it can be included that the movement of CIP particles in MRFs is comprehensively affected by many forces, among which the magnetic field force, the repulsion force and the viscous resistance of the base carrier fluid are the main important elements that affect the movement, the chain-forming process, and the microstructure of MRFs. The CIP particle is acted upon by the magnetic field force, the repulsion force, and the viscous resistance. Therefore, Newton's second law is introduced to establish the kinetic equation of CIP particles without shear action [23], which is written as: Where t is time and ai is the radius of CIP particle i. The kinetic equation of CIP particles under shear action is given by: Specially, in a rotating machinery, the CIP particles of MRFs are subjected to centrifugal force Fcf, as follows: Where, R is the distance between the rotating machinery and the particle center. Accordingly, the kinetic equation of particle i undergoing the shear action in the rotating machinery is expressed as: Thereby, the following equation of CIP particles without shear action can be derived from the kinetic Equation (7): Where, D = 6 a. The time step used in the simulation is Δt, the term on the right side of Equation (11) is a constant and the left side is integrated at an infinitesimal time step. Subsequently, the derivation of particles position change at the time step Δt is given, as follows: Δsi is the variation of particle position at time step Δt, and v0 is the particle velocity at the beginning of time step Δt. When Δt is much smaller than m/D, it is deduced from Equation (12):

Microstructure Simulation
The chain-forming simulation result of CIP particles is conducted via MATLAB to obtain the microstructure and configuration of MRFs. This simulation is performed over a range of conditions that are given in Table 1. The CIP particles are randomly distributed in the base carrier fluid in the absence of magnetic field. The simulation starts with imposing the magnetic field on MRFs, and the CIP particles are rapidly magnetized into magnetic dipoles under the effects of magnetic force, repulsion force, and viscous resistance, collectively. The stress state of CIP particles is deduced based on the kinetic equations that were proposed earlier and then the displacement increment of CIP particles under constant iterations is figured out, according to which the steady state of the particle motion system is obtained, eventually. In this work, the volume fraction of CIP particles is set by 25%. The microstructure evolution of MRFs under magnetic field is simulated and the results are illustrated in Figure 1, the three-dimensional distribution of CIP particles position varying with time under magnetic field. Figure 1a shows the initial state of CIP particles with the random initial position given by the simulation system and zero initial velocity. Figure 1b-d is, respectively, the three-dimensional schematic diagrams of CIP particles distribution in MRFs under the external magnetic field at different time nodes, in which the unit of coordinate axis is meter. The CIP particles are randomly distributed, as it appears from Figure 1a. At first, apply the magnetic field vertically downward along the Z-axis direction. Subsequently, the magnetized carbonyl iron particles attract each other and they begin to move under the synthetic influence of various forces of the magnetic force, viscous resistance, and buoyant force. The particles gradually arrange into short and medium chains. Randomly distributed particles get put into order by degrees. Since the microstructure of MRFs is essentially transformed, the corresponding macroscopic characteristics of MRFs are simultaneously changed from the liquid to the solid-like state. Short and medium branched chains and isolated chains in Figure 1b gradually aggregate with time into long chains, as shown in Figure 1c,d, which is consistent with the chain-forming process observed by previous experiments [24]. Besides, Figure 1d shows the arrangement and aggregation results of CIP particles at 1.68 ms after applying the magnetic field. It is observed that the particles finally form the chains configuration. As time goes on, the motion system of MRFs achieves the steady state. The microstructure dynamic analyses of MRFs and the results of chain-forming in the stable state are obtained, as previously discussed in this paper. In the following subsection, the microstructure model of MRFs under shear action in the steady state is discussed and the constitutive model of MRFs is established.

The Proposal for the Hexagonal Close-Packed Structure
The single-chain structure is mostly adopted in the incipient constitutive models of MRFs, with the merits of simple structure and strong versatility. However, in practical application, the volume fraction of CIP is generally high, and this results in the agglomeration of particle chains, which gets more serious with increasing volume fraction [18,25]. Therefore, it is inaccurate to describe the microstructure of a high concentration MRFs with the single-chain model, because it leads to the failure of MRFs constitutive model derived from this model to simulate the shear yield stress of MRFs appropriately. Some studies have been proposed in order to resolve the shortcomings of single-chain model. One example is the body-centered cubic structure model which can simulate high concentration MRFs. According to the body-centered cubic structure, the CIP particles are infinitely distributed in space in matrix arrangement. However, it is essential to build a more effective model due to the particularity that the universality of body-centered cubic structure is not high when the volume fraction of CIP varies.
It can be seen from Figure 2a that each point can be regarded as a chain of CIP particles. By triangulation, randomly distributed CIP particle chains in the direction perpendicular to the magnetic field forms pentagons, hexagons, or heptagons. Among which, the hexagons are the majority. This situation is more evident in high concentration MRFs. Accordingly, the irregular polygon is simplified to regular a hexagon to facilitate the calculation. Thus, in this work, the hexagonal closepacked structure is proposed. Figure 2b is the schematic drawing of the hexagonal close-packed structure. In this structure, each CIP chain has six adjacent chains and the distance between two chains is equal. The distance between adjacent chains is greatly affected by the volume fraction of CIP particles: the distance decreases with increasing volume fraction and it increases with decreasing volume fraction. In each chain, the magnetic dipole is impacted by other magnetic dipoles in this and adjacent chains. On this basis, the constitutive model of MRFs is to be constructed in this work. The salient feature of this method is to take a single chain as a unit [26] and then consider the influence of adjacent chains to adapt the hexagonal close-packed structure to the concentration change of MRFs.

Analyses of the Constitutive Model
In the previous studies, a single chain is separated from others to only discuss its own influence and ignore that from the adjacent chains. In this work, the effects of magnetic dipoles in the single chain and the adjacent chains on the single chain itself are considered through the hexagonal closepacked structure, and then a new constitutive model of MRFs considering the adjacent chains is formulated in order to resolve this problem.

Force Analyses of Particles in a Single Chain
The dipole theory and the Maxwell stress tensor theory represent many studies on the constitutive model of MRFs [27,28], and it is known that the dipole theory is relatively concise and practical. Thus, this work builds the shear yield stress model of MRFs on account of it. The assumptions are made as followings: (1) the size and material of all CIP particles are the same; (2) the particles shape is approximately a sphere. With these premises, CIP particles are regarded as magnetic dipoles under the action of magnetic field. The magnetic energy model of a magnetic dipole in a single chain is constructed with the interaction between two magnetic dipoles into account. According to the magnetic dipole theory [29], the magnetic energy that is produced by arbitrary magnetic dipole j at i can be expressed, as follows: Where mi is the magnetic dipole moment of the i-th magnetic dipole. The mi can be expressed, as follows: Where, V is the volume of a single CIP particle and M is the magnetization. Equation (14) indicates that the magnetic energy of the magnetic dipole j at a random position is inversely proportional to the third power of the distance from the magnetic dipole i. The interaction between magnetic dipoles is relatively strong along the magnetic field direction. Additionally, the force of the two interacting magnetic dipoles both in the same chain is greater than that in two different chains because the distance between the magnetic dipoles in the same chain is smaller than that in the adjacent chain. Suppose that there is an infinite amount of CIP particles in a single chain. Subsequently, by inserting Equation (15) (16) Where, n is the number of magnetic dipoles in a single chain, t the surfactant coating thickness of a CIP particle,  the net distance between the adjacent CIP particles in the same chain, and  is the deviation angle of the magnetic dipole under shear, which are specifically shown in Figure 3a.
The influence coefficient of the magnetic dipole in the same chain is presented by the following equation: The space occupied by one single CIP particle.

Analyses of the Influence from Adjacent Chains
For the constitutive model of high concentration MRFs, there is a need to consider the magnetic energy that is produced by the magnetic dipole not only at i in the same chain, but also in the adjacent chains. It is generally considered that the distance between the adjacent chains is greatly affected by the volume fraction of CIP particles. Taking the hexagonal close-packed structure as the basis, the distance  between CIP particles chains can be deduced. The space occupied by a single CIP particle is a hexagonal prism, as shown in Figure 3b. Then, the volume ratio of the CIP particle to the hexagonal prism is the volume fraction  of the CIP particles of MRFs, which is given by: Where, VL represents the volume of the hexagonal prism unit. At the same time, it can be seen that the height of the hexagonal prism is exactly h, and it follows that the volume of the hexagonal prism unit can be expressed, as follows: Subsequently, the bottom area SL of the hexagonal prism is obtained by solving the above two Equations (18) and (19). According to the top view of the hexagonal prism space occupied by CIP particles shown in Figure 4a, the following equation is achieved that contributes to deducing the side length p of the bottom of regular hexagonal prism: Subsequently, the distance  between CIP particle chains can be written as: The central distance a between CIP particles in the same chain is taken as 2.2a in order to simplify the calculation. The ratio between dipolar and repulsion force of two particle chains along the magnetic field should not be too large or too small in order to obtain the optimum shear yield stress.
 should be bigger than 2.2a which is proved in previous studies and practical experience, because the ratio can reach the value 10 (0.1) when the distance between particle chains increases to  / 2a = 1.1 (however, decreases to /2a = 0.9). [14,20]. The upper limit of volume fraction of the presented model comes out 45.4% by taking a as 2.2a into Equation (21). Moreover, the volume fraction of CIP particles in MRFs is usually less than this value due to the restriction of zero magnetic field viscosity and other conditions. Therefore, the proposed model meets the practical application conditions. It is verified from Equation (16) that the influence of magnetic dipoles in adjacent chains is much smaller than that in the same chain because the distance between magnetic dipoles in adjacent chains is larger than that in the same chain. Only the influence of the chains that directly adjoin to the chain of the magnetic dipole i is considered, while the influence of farther chains is ignored to predigest the constitutive model. The expression of magnetic energy of the magnetic dipoles in adjacent chains to the magnetic dipole i is derived by Equation (15) where, Ei2 and Ei3 represent, respectively, the magnetic energy along the shear direction produced by the magnetic dipole, which is far from and closer to the magnetic dipole i after deviation under the shear action, Ei4 and Ei5 represent, respectively, the magnetic energy along the shear direction produced by the magnetic dipoles on both sides, which are far from and closer to the magnetic dipole i after deviation under the shear action. It is noted that the magnetic dipole only moves along the shear direction. As expressed in Figure 4b,  and  represent the angle between the offset magnetic dipole i and the far and close magnetic dipole in the vertical direction, while  and  separately denote the angle between the offset magnetic dipole i and the far and close magnetic dipoles on both sides in the vertical direction. Figure 4c illustrates the angles between the front and rear magnetic dipoles and the offset magnetic dipole along the shear direction. To obtain those angles, the edge lengths of two rightangled sides of the right-angled triangle are first solved according to Figure 4c. Subsequently, applying the Pythagorean theorem yields the following result:

Constitutive Model of MRFs
As explained above, Ei the total magnetic energy of magnetic dipole i is expressed as: Subsequently, the magnetic energy that a single chain in unit length in the hexagonal closepacked structure has is solved as: Where, L represents the length of a single chain in a hexagonal closed-packed structure. The shear yield stress of the single chain is obtained by taking the derivative of Equation (31) [30]: Note that the shear yield stress of MRFs is comprised of the sum of shear yield stress that is formed by every single chain within the unit area [31]: (34) Where, N is the number of chains of CIP within unit area and Vm is the MRFs volume. However, the CIP particles in this paper are hypothesized to be ideal, considering them to be equally sized and sphere-shaped and uniformly distributed in MRFs. However, the non-uniformity of CIP particles distribution and difference of the particle shape and size do exist, thus, the error that is caused by the non-uniformity of CIP particles distribution should be taken into consideration when the MRFs constitutive model is established. The parameter k is introduced to express the effect of in homogeneity of CIP particles on the shear yield stress of MRFs; the updated MRFs shear yield stress is as follows: The complete form of Equation (35) is complicated and inconvenient, so, to simplify the expression, the simplified form is written as Equation (35). The function f() in Equation (35) is derived in the Appendix A.
The unit of data measured of MRFs external magnetic field is magnetic induction intensity, which is substituted for the magnetic field intensity H in Equation (35), which is not measurable by instrument [32]. When considering the magnetization of paramagnetic particles, the magnetic induction intensity B is expressed by the following formula: The Froelicb-Kennelly material model [33] is adopted to write the relative magnetic susceptibility because of the nonlinear and saturated magnetization of CIP particles under magnetic field:

Experiments
The shear yield stress of MRFs is affected by many factors, as mentioned above. It can be seen that magnetic induction strength, size, and volume fraction of CIP particles and the surfactant coating thickness are the dominant influences on the shear yield stress of MRFs, according to the constitutive model of MRFs proposed. The following sections investigate the unknown parameters by the fitting method. In addition, the variable-controlling approach is used to explore these factors. Now some rules of experiments are set as follows： • The magnetic induction intensity is 0.4 T.
• The volume fraction of CIP particles is 34%.
• Read the shear yield stress when the shear strain is 0.31.
• The relationship between the surfactant coating thickness and distance of particles is expressed as 2t +  =  m.
• The radius of all CIP particles is 4 m.
• The temperature to conduct the experiments is 25 °C.
It is noted that, when the effect of a certain factor to the properties of MRFs is investigated, the other influenced factors are assumed to remain as unchanged according to the above settings.

The Preparation of MRFs for Experiment
MRFs is usually composed of additives and micron-sized magnetic particles that are uniformly distributed in water-based or oil-based carriers. The base carrier liquid is generally required to have a higher boiling point and lower freezing point to meet different application temperatures. The zero field viscosity of MRFs is usually greatly affected by the viscosity of base carrier which is generally proportional to its density, so the base carrier needs to have lower viscosity. When considering the sedimentation rate of the MRFs, the density of the base carrier should not be too low to avoid the density difference between the magnetic particles and the base carrier to be too large. In addition, in practical application, whether the oil-based or water-based carrier is chosen is on the basis of whether the material that the MRFs works on is water-soluble or oil-soluble. In practical application, many materials are water-insoluble and, according to Table 2, as below, water-based carrier fluid was selected for preparing the MRFs. In the magnetic field, the rheological properties of MRFs mainly depend on the magnetic particles. The shear yield stress of MRFs is the macroscopic reflection of the force between the magnetic particles under the magnetic field. Therefore, the magnetic particles should possess high permeability and magnetic saturation strength to ensure the high shear yield stress of MRFs. At the same time, it is necessary to have a low hysteresis so that the MRFs can quickly change into Newtonian fluids after the external magnetic field disappears, and the viscosity can quickly return to zero-field-state viscosity for circulating the MRFs. Under comprehensive consideration, CIP, meeting the above requirements and its simple manufacture, is selected for preparing the MRFs as the magnetic particles.
The solid magnetic particles usually disperse in the base carrier, but they are difficult to dissolve in the base carrier due to the large density difference between them, which is prone to cause the settlement of magnetic particles. Therefore, surfactant needs adding into MRFs to better the settlement stability. The surfactant in MRFs should have special molecular structure: one end is the functional group with high affinity for absorbing solid magnetic particles; the other end is the elastic group that is easily dispersed in the base carrier, which prevents the magnetic particles from approaching each other and agglomerating by space repulsion it produces, and at the same time reduces the sedimentation. The magnetic CIP particles in the MRFs are easy to be corroded and oxidized by the base carrier, which weakens the magnetic saturation strength of the magnetic particles. The surfactant coated on the surface of magnetic particles can effectively prevent it. Besides, the deionized water is chosen as base carrier mentioned above, thus sodium dodecyl sulfonate is selected as surfactant. Table 3 summarizes the compositions of MRFs. So far, MRFs can be prepared according to the selected compositions above and Table 4.

Shear Stress-Strain Relationship
The shear yield stress of material is the maximum shear stress to resist deformation, which is reflected as the maximum stress point on the shear stress-strain curve. For MRFs, it is considered that the chain structure reaches its maximum deformation under the maximum shear stress. This group of experiment simulates the relationship between the shear stress and strain taking the strain as variable while the other factors unchanged. Figure 6 depicts the shear stress-strain relationship of MRFs in view of its chain structure. This result tells the point of maximum shear stress and the corresponding strain with a value of 0.31. The trend of the shear stress-strain relationship of other volume fraction of CIP particles, such as 25% and 40%, is also the same as the experimental result of 34%; the strain is almost often 0.31 when the shear stress reaches its maximum. Therefore, the shear stress is taken as the shear yield stress of MRFs in the subsequent experiments when the shear strain is 0.31.

Effect of Magnetic Induction Intensity
This group of experiment is designed not only to research the relationship between the shear yield stress and magnetic induction intensity, but also verify the strong adaptability of the constitutive model to the change of volume fraction of CIP particles. Thus, the groups of MRFs with different concentration are prepared to carry out this experiment and obtain the experimental data of MRFs shear stress at different concentration varying with the magnetic induction intensity. This group experiment is divided into three parts differing in the volume fraction of CIP particles. The magnetic induction intensity is taken as a variable in each group, while the other factors remain unchanged. Compare the experimental result with the Simplified Single-chain Model [34]. Figure 7a-c illustrate the curves of relationship between different volume fraction of CIP particles and magnetic induction intensity, respectively. In Figure 7, the 'Experimental Model' is obtained by fitting the experimental data that are measured in this work by conducting groups of experiments with MRFs prepared according to Section 4.1. The abscissa values are substituted into the constitutive model deduced above to obtain the corresponding ordinates, and the 'Constitutive Model' is formed by these points. The same is true for the 'Simplified Single-chain Model'. Note that the 'Simplified Single-chain Model' is a simplified form, which definitely differs from the classic "Single-chain Model", of force analysis of one single chain, according to the Section 3.3.1. It is simplified as compared with the "Constitutive Model", which considers the effects of adjacent chains. Figure 7a is the fitting diagram of known model and Figure 7b is the simulation result of unknown model.
From Figure 7a-c, it is validated that the constitutive model constructed in this work well fits the trend of the shear yield stress as a function of the magnetic induction intensity. The saturation point of the shear yield stress moves backward with the increasing volume fraction of CIP particles, which is consistent with the previous research. Additionally, it is obviously seen from Figures 7b,c that the simulation of single-chain model to high concentration MRFs is not good showing relative larger error than the proposed constitutive model. It is remarked that the constitutive model that was built in this work can predict the shear yield stress well under magnetic saturation of MRFs, which is, in general, hard to achieve from experimental tests.

Effect of CIP Particles Volume Fraction
The volume fraction of CIP particles is taken as the variable while other parameters are kept to be unchanged to perform this group of experiment. Based on the constitutive model, the relationship between the shear yield stress and the volume fraction of CIP particles of MRFs is simulated and shown in Figure 7d. It is clearly observed that there is a linear relationship between shear yield stress and CIP volume fraction: the MRFs shear yield stress increases with the increasing CIP volume fraction, but the trend slowly becomes smooth, which means that the slope reduces. When the CIP particles deviate under the action of shear, the offset particles will approach to the CIP chain along the direction of deviation. According to the previous analyses of this work, at this moment, the force that is exerted on the offset particles by the chain along the deviation direction is greater than that exerted by the chain in reverse direction. Moreover, with the increasing CIP volume fraction, the distance between the particle chains decreases, which leads to the stronger influence of adjacent chains. Thus, the increasing trend of shear yield stress of MRFs decreases with increasing the volume fraction of the CIP particles.

Effect of CIP Particle Radius
The effect of particle radius on the shear yield stress is studied in this group of experiments. With the fitting parameters determined and other conditions remaining unchanged, the changes of MRFs shear yield stress with CIP particle radius is simulated based on the model proposed and shown in Figure 8a. It is seen that the shear yield stress obviously increases with the increasing radius of CIP particles. Specifically, the magnetic interaction enhances with the increasing CIP particle radius, which leads to the increase of MRFs shear yield stress. At the same time, this trend is subjected to the increased center distance of particles, causing the diminished magnetic interaction. When the particle radius is relatively small, this increasing trend is more obvious; when the particle radius relatively is large, this trend slows down. The MRFs sedimentation stability is also affected by the changing radius of CIP particles, which is illustrated by the increase of sedimentation ratio with an increasing particle radius. Therefore, practically, the particle size should be carefully selected to meet the imposed requirements of MRFs, because a bigger particle size does not mean a better result of the properties of MRFs.

Effect of Surfactant Coating Thickness
Generally, the layer of surfactant coated on the CIPs is taken as even and the thickness of surfactant coating is regarded as an invariant for simplifying the calculation and facilitating the simulation [35]. In fact, with the increase of the surfactant volume fraction, the thickness of surfactantcoated particles increases. It is noted that when the surfactant volume fraction increases to a certain level, the coating thickness will stop increasing and the rest surfactant will have a negative influence on the properties on MR fluids. Therefore, an optimum volume fraction of surfactant is of great significance. The coating thickness is used to express the amount of surfactant that is adsorbed on or around the CIP particles. In this group of experiments, in the circumstances of other variables keeping invariant, the surfactant volume fraction is enhanced in order to thicken the surfactant coating. Subsequently, the thickness change of surfactant coating with MRFs shear yield stress is simulated and shown in Figure 8b. It is seen that the shear yield stress slowly decreases with an increasing surfactant coating thickness. This is attributable to the increscent distance between CIP particle centers that results from the increase of surfactant coating thickness. Although the surfactant takes up a small proportion in MRFs, it has great influence on the sedimentation stability. The sedimentation ratio of particles in MRFs decreases with the increase of surfactant coating thickness, which further leads to the reduction of MRFs shear yield stress. Through practical experience, the surfactant coating of CIP particles after reaching a certain thickness usually no longer increases. Therefore, the surfactant volume fraction should be selected based on that of the CIP particles, not being too large.

Conclusions
In this research, a constitutive model for predicting the shear yield stress of high concentrated MRFs has been proposed and its effectiveness is validated through comparative works with a singlechain model and experimental result. The microstructure of the particle chains of MRFs with high concentration was modeled by choosing the hexagonal close-packed structure. In the validation of the proposed model, magnetic induction intensity, volume fraction, size of CIP particles, and surfactant coating thickness affecting shear yield stress were adopted in order to investigate the influencing weight to the properties of MRFs. When one experimental group is under testing, the other factors have been kept to be unchanged to accurately observe the effect trend of each factor. It has been demonstrated through both simulation and experiment that the proposed constitutive model can predict the shear yield stress of high concentrated particles based MRFs much better than the single-chain model, which has been conventionally and frequently used so far. In other words, the presented constitutive model highly coincides to the experimental results showing a very similar trend of the shear yield stress of MRFs.
It should be asserted that the proposed model based on the hexagonal close-packed structure of CIP chains is simple and novel. In addition, the constitutive model that is presented in this work provides a qualitative tool for readers to understand the microstructure and mechanical behavior of CIP particle chains in MRFs with high concentration. However, the model itself still exists complexity in part with several ideal assumptions that are made in the process of mathematical modeling. Therefore, some modifications of the proposed model with practical conditions are required to enhance the prediction accuracy and provide a theoretical basis for the application products of MRFs.