 1.1 Microscopic Description of the Dynamics
 1.1.1 Stokes Approximation
 1.2 Magnetic Interactions
 1.2.1 Mean Magnetization Approximation (MMA)
 1.3 Other Nonhydrodynamic Interactions
 1.4 Rheology
 1.4.1 Structure
 1.4.2 Preyield Regime. Static Yield Stress
 1.4.3 Postyield Regime. Shear Thinning Behavior
 1.4.4 Normal Stresses
 1.5 Conclusions
 Abbreviations
 References
Chapter 1: Introduction to Magnetorheological Fluids

Published:07 Jun 2023

Special Collection: 2023 ebook collectionSeries: Soft Matter Series
J. R. Morillas and J. de Vicente, in Magnetic Soft Matter
Download citation file:
Magnetorheological fluids are multiphase magnetizable suspensions with magnetic fieldcontrollable mechanical properties. In this introductory chapter we revisit the physics behind the rheological response of these particular materials making special emphasis on the influence of magnetostatic and hydrodynamic forces.
As happens to any kind of suspension, the macroscopic behavior of magnetorheological fluids (MRFs) is intimately connected to the phenomena that take place at the microscale. In this chapter, the theoretical basis necessary to describe and understand the particle microconfiguration, its dynamics and how these two are related to the MRF collective response are summarized. In the following, MRFs will be treated as suspensions of N magnetic monodisperse spherical particles (diameter d) in a continuous phase occupying a total volume V. This chapter mainly focusses on conventional MRFs (CMRFs), constituted by a dispersion of micronsized magnetizable particles in a Newtonian (nonmagnetic) liquid carrier. In these particular MRFs Brownian motion can be neglected.
1.1 Microscopic Description of the Dynamics
For any colloid or suspension (i.e. even those that are not magnetic), the size of the dispersed particles is large enough to consider the carrying continuous phase as a fluid with a welldefined viscosity η_{c}. When a dispersed particle moves within the continuous viscous phase, for example under the action of an external force, it will distort the fluid velocity and pressure fields. Thus, the new velocity and pressure fields will drag neighboring particles changing their positions. With this in mind, any particle acts over its neighbors through a fluidmediated force, the hydrodynamic one F⃗_{h}, that would not exist if the particles were in vacuum.
Generally speaking, hydrodynamic forces are not pairwise analytical interactions and hence they cannot be written solely in terms of the relative positions or dynamic state (i.e. velocity and/or acceleration) of the particles. Instead, since they are fluidmediated forces, they depend on the dynamics of the surrounding fluid. Consequently, the fluid flow problem must be solved first.
In most cases, the carrier fluid is found to be incompressible (i.e. it is a liquid). Thus, the massbalance (continuity) equation reduces to:
where v⃗ is the fluid velocity field. In addition, if the fluid satisfies Newton's viscosity law (i.e. the fluid is Newtonian), the hydrodynamic pressure p and velocity fields are coupled by the momentumbalance (Navier–Stokes) equation:
where ρ_{c} is the fluid density and f⃗_{ext} is the body force (per unit volume) exerted by any external agent. Furthermore, if the flow is not isothermal, an additional energybalance differential equation must be solved. For incompressible fluids this equation is uncoupled (i.e. it does not depend on velocity or pressure fields) and can be solved a posteriori. Nevertheless, through this chapter temperature dependencies and heat transfer phenomena are not considered. As a result, the energy equation does not need to be solved.
In order for the problem to be wellposed, eqn (1.1) and (1.2) must be completed with boundary conditions. These depend on the specific geometry under study. However, typically it is necessary to impose uniform velocity fields as farfield conditions or periodicity (for unbounded media), a prescribed fluid velocity (the socalled noslip condition) or a pressure on the confining walls (for bounded media).
Finally, the coupling between the carrier fluid and the particle movement is accounted for by imposing noslip conditions on all particle surfaces, that is, v⃗ must meet the particle velocity on its surface:
here u⃗ is the translational velocity of the particle, ω⃗ is its angular velocity and r⃗_{S} is a vector joining each point of the particle surface (where v⃗ is evaluated) with the particle center.
The particle velocity is given by Newton's second law:
where ρ_{p} is the particle density, V_{p} = πd^{3}/6 is the particle volume, Ĩ is its inertia tensor and F⃗ (T⃗) is the total force (torque) acting over the particle. This force (torque) can be split into the hydrodynamic force (torque) F⃗_{h} (T⃗_{h}) plus contributions coming from nonhydrodynamic effects F⃗_{nh} (T⃗_{nh}), that is, all those forces (torques) that the particles would experience if they were not immersed in a fluid.
Regarding the hydrodynamic interaction, the total force exerted by the fluid over the particle is given by the integral, over the particle surface S_{p}, of the fluid total stress tensor t̃:
where δ̃ is the secondorder identity tensor and Ẽ is the strain rate tensor which can be written in terms of the velocity field as Ẽ = [∇v⃗ + (∇v⃗)^{T}]/2. Therefore:
here n̂ is the surface normal vector pointing to the fluid. Due to the continuous nature of the fluid, this does not only exert a punctual force over the particle but a continuous distribution over its surface. As a consequence, the hydrodynamic interaction is not reduced only to a total force; the fluid also exerts a moment. In its more general form, the first moment of the hydrodynamic force can be written as follows:
where the symbol is understood as the dyadic product of the vectors t⃗_{n} and r⃗_{S} (note that given two vectors a⃗ and b⃗, the secondorder tensor would read in index notation a_{i}b_{j}). However, in fluid dynamics it is more useful to split M̃ in its symmetric and antisymmetric parts. The former part is called the stresslet:
and it does not have any effect on the dynamics of rigid particles (see below). On the contrary, it can be shown that the antisymmetric part Ã = (M̃ − M̃^{T})/2 is just the torque due to F⃗_{h}:^{1 }
and thus, responsible for the particle rotation.
As can be seen, the coupling between particle and fluid motion is in both directions. The velocity of the fluid v⃗ depends on the particle translational u⃗ and angular ω⃗ velocities through the noslip condition (eqn (1.3)) and the particle velocities depend on the fluid velocity, through the hydrodynamic force and torque (eqn (1.7)–(1.10)), that accelerates/decelerates the particles.
Eqn (1.1)–(1.10) govern the dynamics of the (Newtonian) continuous and (solid) dispersed phases using first principles. However, the resultant differential equation system is extremely complicated to solve. Firstly, just regarding the fluid dynamics, the Navier–Stokes equation is a nonlinear partial differential equation and approximations are needed to solve it. Secondly, a manybody problem with nonlinear and irreversible features results when a dispersed phase is incorporated to the fluid. In addition, due to the absence of Brownian motion, the problem becomes deterministic and requires tracking the explicit movement of a vast number of particles precluding any analytical treatment.
Tackling the problem from a numerical point of view is not easy. Most common Computational Fluid Dynamics (CFD) tools are meshbased techniques where the fluid domain has to be discretized (in mesh elements). In this way, the governing partial differential equations (eqn (1.1) and (1.2)) can be transformed in algebraic ones to be solved by a computer. However, in MRF problems, the fluid domain continuously evolves when the particles change their position. Consequently, at every time step, both the computational domain and the mesh should adapt to the new particle configuration. Furthermore, in order to evaluate the time derivative in eqn (1.2) at the present time step, the velocity field in the previous time step has to be projected onto the new mesh (by interpolating the velocity field).^{2 }
Taking this into consideration together with the fact that: (i) the mesh has to be dense enough around every particle (to properly compute the force and torque acting on them), (ii) the great difference between the involved lengths scales (in a CMRF, the typical particle diameter is around 1 µm but the sample size is usually around 0.3–1 mm) and (iii) the large number of unknowns (pressure plus three velocity components) one can see that solving eqn (1.1)–(1.10) is not straight forward and some approximations are required.
1.1.1 Stokes Approximation
In many cases, the study of MRFs is done under conditions where both fluid and particle inertia are negligible in comparison to the viscous dissipation at both sample and particle length scales. Therefore, the Navier–Stokes equation is reduced to the the Stokes equation:
In this expression, external body forces are assumed to be conservative and therefore can be included in the pressure term (LHS). Consequently, p will not stand for the absolute pressure and derived forces exerted by the fluid (see eqn (1.6) and (1.7)) will not account for the effects of that body forces. As can be seen, the Stokes equation is linear, and it does not contain any explicit dependence on the time. Consequently, although a new mesh is required at each time step, the flow field in the previous steps is not needed. It only depends on the current particle configuration and particle velocities (that are introduced as boundary conditions, that is, known values from previous step).
Although the use of eqn (1.11) instead of eqn (1.2) allows cessation of projecting and time deriving the flow field, solving the MRF dynamics by meshbased CFD techniques still implies a large computational effort due to the large number of nodes present in the meshes. During the evolution of the system the highest velocity gradients, and with them the highest forces (see eqn (1.6) and (1.7)), are localized in the gaps between solid surfaces and confining walls (if applied). In order to accurately resolve the velocity gradients, these gaps should be densely meshed, increasing the computational cost. The situation is aggravated with the number of particles and geometry complexity restricting meshbased CFD simulations to very simple systems with a small number (tens) of particles and/or 2D domains.^{3–5 }
This problem can be partially overcome using techniques that do not require adaptive meshes such as the Lattice Boltzmann Method (LBM)^{6,7 } where an initial mesh can be defined without further rearrangement because the flow field is only solved on those nodes that are not occupied by the particles. Notwithstanding, a dense mesh is still needed to solve plausible field gradients. Finally, meshfree methods such as Smooth Particle Hydrodynamics (SPH)^{8,9 } or Dissipative Particle Dynamics (DPD)^{10,11 } have been used to study MRFs as well. The basic difference with respect to meshbased CFD methods lies in regarding the continuous phase not as a continuum but as an ensemble of discrete ‘fluid elements’ at the mesoscale. These interact with each other in a way that, once time and spaceaveraged, are able to reproduce the fluid macroscopic behavior. Nevertheless, these meshfree methods base their accuracy on the ‘fluid element’ size, that is, the smaller the ‘fluid element’ the greater the resolution of the flow field around the particles. Thus, although they are more versatile and some computational effort can be saved, the necessity of a very fine discretization (of the continuous phase instead of the computational domain) to resolve the flow problem does not make them the best choice when dealing with MRFs either.
From the previous discussion, it seems that solving explicitly the flow field issue around the particulate phase is an unaffordable task from a computational point of view. To circumvent this problem, other methods have been proposed to take advantage of the linear properties of the Stokes equation. As a first approximation, any flow field, regardless of its complexity, can be expressed around a given point r⃗_{0} by its Taylor series:
where and it should be understood that the velocity V⃗_{0}, the shear rate tensor Ẽ_{0} and the vorticity are evaluated on r⃗_{0}. As can be seen, according to eqn (1.12) any flow can be approximated, until first order, as the sum of simpler flow fields (Figure 1.1): a homogeneous flow field V⃗_{0}, a pure straining flow field Ẽ_{0}·(r⃗ − r⃗_{0}) and a pure rotational flow field .
Linearity implies that the effect of the total flow field over the immersed particles can also be split into the individual effects that each of the previous simpler flow fields has on the particles separately. In this way, for any flow field, one can evaluate the associated total hydrodynamic force/torque over a particle by addition of the corresponding force/torque sustained by the particle under a homogeneous, pure rotational and pure straining flow.
With this, the problem reduces to solve the total flow field (eqn (1.1) and (1.11) together with the boundary condition, eqn (1.3)) around a spherical particle under these simpler flows. It can be shown that, applying eqn (1.7) and (1.10) at the particle surface, hydrodynamic force and torque over the particle can be written until first order as follows:^{1 }
Eqn (1.13) is commonly also known as the Stokes’ drag law. Several points can be highlighted from eqn (1.13) and (1.14):

The hydrodynamic force (torque) comes only from the homogeneous (pure rotational) flow field, that is, a spherical particle immersed in a homogeneous (pure rotational) flow does not experience any hydrodynamic torque (force).

The pure straining flow does not exert any force or torque over the particle.

The flow field and particle movement are supposed to be noninertial. Thus both force and torque can be expressed as a function of the relative velocity, that is, it is always possible to find a reference frame where the particle or the fluid are quiescent.

A linear dependence exists between the force (torque) and the translational (rotational) velocity because the Stokes equation is linear.
From the previous discussion, it could be inferred that the pure straining flow is not affected by the presence of the particles as no net force and/or torque are developed. However, this result arises from the spatial symmetry of the straining flow field around a spherical particle. Actually, particles are really under a stresslet only due to the pure straining flow component:
but, as solid bodies, they can bear any stresslet without straining (Ẽ = 0 inside the particles), and consequently S̃ does not play any role in their movement. Note also that, because of this, the stresslet in eqn (1.15) is not expressed in terms of the relative strain field as it happens for the force and the torque (eqn (1.13) and (1.14)).
Until now, it has only been considered that the continuous phase has a negligible inertia. Taking into consideration the small size of the particles, it can also be supposed that the particles are inertialess, and hence, the total force, F⃗, and torque, T⃗, acting on them must be zero. According to eqn (1.4) and (1.5), this implies that the hydrodynamic force/torque is always balanced by the nonhydrodynamic forces/torques: F⃗_{h} = −F⃗_{nh} and T⃗_{h} = −T⃗_{nh}. Introducing these conditions in eqn (1.13)–(1.15) yields:
Note that in this new set of eqn (1.16)–(1.18), force and torque values do not depend on the kinematic state of the continuous phase as happens in eqn (1.13)–(1.15). Now the forces and torques are the nonhydrodynamic ones that usually depend on the particle positions.
Eqn (1.16)–(1.18) govern the particle dynamics in the suspension (whenever inertia is neglected), but it should be born in mind that, in these expressions, the quantities V⃗, and Ẽ are not simply constants. They contain the homogeneous, pure rotational and pure straining parts of the total flow at the particle position, that is, the flow externally imposed plus the perturbation flow field created by each dispersed particle at that position. Focusing only on the translational motion, it is clear that the homogeneous part of the flow field created by particle α, V⃗_{α}, will depend on its translational velocity and nonhydrodynamic force, V⃗_{α} = f⃗_{α}(u⃗_{α},F⃗_{nh,α},…).
However, if particle α is also rotating at ω⃗_{α} and experiencing a nonhydrodynamic torque T⃗_{nh,α}, it will create another flow field that will have homogeneous, pure rotation and pure straining components, all of them dependent on ω⃗_{α} and T⃗_{nh,α}. In the same sense, if the particle α is under a straining flow, it will create another flow field with its corresponding components that depends, in this case, on the stresslet S̃_{α}.
Therefore, when all homogeneous components are joined together, it yields V⃗_{α} = f⃗_{α}(u⃗_{α},F⃗_{nh,α},ω⃗_{α},T⃗_{nh,α},S̃_{α}) with the main advantage that the previous dependence is linear in all variables due to Stokes flow properties. Now extending the previous reasoning to rotational and straining components induced by the rest of the particles we arrive to:
Since each function (f⃗_{α}, g⃗_{α}, h̃_{α}) is linear in its variables, when they are introduced in eqn (1.16)–(1.18) they can be rearranged so that the particle dynamics governing equations can be written as a matrix equation:
where the bold symbols consist of the corresponding magnitude of each particle appended one after another, that is, F⃗_{nh} = (F⃗_{nh1}, F⃗_{nh2},…F⃗_{nhN}). Note that with this notation, each element of is not a scalar but a tensor of the proper order to satisfy the linear relationship between the dynamic and kinematic variables. The terms V⃗_{∞}, and Ẽ_{∞} are usually known because they are imposed macroscopically, hence the solution is expressed in terms of the relative velocities.
is known as the grand resistance matrix. It consists of the proper combination of the f⃗_{α}, g⃗_{α} and h̃_{α} functions to relate the particle kinematics (translational and rotational velocities) and fluid strain to the dynamics acting over a given particle. Importantly, both ^{12 } and the LHS of eqn (1.20) only depend on the relative position of the particles, therefore the kinematics of each particle can be obtained by evaluating both terms in the current position of the particles and computing . Once the velocities are solved, the next particle configuration can be obtained by integration over time. This scheme constitutes the basis of the Stokesian Dynamics (SD) technique.^{13 }
As it can be seen, Stokes flow linearity allowed the construction of a matrix system to solve the particle dynamics but nothing has been said about how to compute the elements within . Frequently, these are obtained by computing the flow field around a pair of particles that is subjected to the three aforementioned components of a general flow: homogeneous, pure rotational and pure straining. Depending on the interparticle distance this total flow field is solved using a reflection/perturbation method (large distances above 2d), the lubrication approximation (short distances, from 1d to 1.01d) or tabulated values coming from simple CFD simulations for intermediate interparticle distances.^{14,15 } Note that according to these methods, the elements within really account for twoparticle interactions only. Nevertheless, it is stated that during computation, these elements are combined, and thus manybody interactions are finally considered.^{13 }
It is not the scope of this chapter to find the particular expressions for the elements within . Nevertheless, it is worth mentioning two of their basic features. On the one hand, at long interparticle distances R_{αβ}, the flow field and, with it, the hydrodynamic force induced by particle β at the position of particle α scales as ∝R_{αβ}^{−1}. Clearly, this is a longrange force that obliges one to consider force contributions even from particles that are widely separated.^{16 }
Conversely, at short R_{αβ}, the main hydrodynamic force is due to the lubrication. For the case of two particles approaching at a constant velocity U along their common axis the force is repulsive and scales as ∝U/(R_{αβ} − d). As can be seen, the lubrication force diverges with the surfacetosurface distance, meaning that particle contacts take an infinite time to happen. Actually, contacts may occur in real systems because of surface roughness and the existence of nonhydrodynamic interactions.
The computation of may be a laborious process. However, it can be implemented numerically without a high computational cost since it starts from already known/tabulated functions. Moreover, the resolution of particle dynamics using SD does not imply any approximations (whenever inertia can be neglected) regarding hydrodynamic interactions but, at the same time, does not require explicitly solving the continuous phase dynamics. As a result, it becomes a molecularlike treatment of the problem where only particle–particle interactions (dependent on the relative positions) are considered.
Although molecularlike SD simulations require a smaller computational effort in comparison to meshbased CFD ones, they are still computationally expensive as a dense matrix (remember that hydrodynamic interactions are longranged and thus couple the motion of all dispersed particles regardless their separation) inversion is needed for each time step. Looking at eqn (1.20), it can be seen that in order to follow the system evolution with time, it is necessary to solve only the first two rows concerning translational and rotational velocities of the particles. Therefore, from a computational point of view, the more expensive step is to invert a matrix of dimension 6N × 6N, where N is the number of simulated particles and 6 refers to the total number of components of the force and torque together. Note that the last row in eqn (1.20) concerning stresslets is not needed if in doing so, only the particle configuration is pursued.
Dense matrix inversion constitutes a strong limitation to the implementation of SD simulations to study MRF systems with only hundreds of particles, usually in 2D.^{17,18 } If SD are applied outside these borders, some kind of approximations in the hydrodynamics interactions, neglecting their short or longrange effects, for example, are needed.^{19–22 }
As commented above, available methodologies described so far to solve hydrodynamics (based on the Stokes or the Navier–Stokes equations) are too expensive from a computational point of view for modelling MRFs. This is essentially because the number of simulated particles to properly resemble real systems, and do averages, is typically above 1000. Thus, most frequently used approaches aim to solve particle and fluid dynamics separately, designing (or not) a coupling scheme to obtain consistent solutions for the whole system. In these approximations, particles move in the same spatial domain occupied by the fluid, that is, there is not a volume exclusion and particles do not displace the surrounding fluid during their motion.
Particle dynamics are governed again by Newton's second law (eqn (1.4)), usually neglecting rotational degrees of freedom. Although different expressions can be taken in order to consider neighboring particles, inertial flow and other effects on the hydrodynamic force,^{23,24 } the latter is typically given by the Stokes’ drag law:
where v⃗ is the fluid velocity at the particle position. This velocity together with the pressure fields are computed by solving continuity equation (eqn (1.1)) and a modified Navier–Stokes equation as follows:
Eqn (1.22) differs from eqn (1.2) in the term f⃗_{c}. It stands for the force exerted by all the immersed particles over the fluid and is responsible for the coupling between particle dynamics and fluid flow problems.^{25 } In order to satisfy Newton's third law we have:
where V_{c} is the fluid unit volume and the summation on α runs over those particles inside V_{c}. This scheme (usually known as twoway coupling) roughly considers the effect of a particle on the fluid that surrounds it reducing/increasing the fluid velocity to adjust it to the particle one. The manybody problem of particle dynamics is usually solved through molecularlike simulations (such as molecular dynamics or discrete element method) while the flow problem is addressed using meshbased CFD (in this case, V_{c} in eqn (1.23) is the mesh element volume) or SPH (V_{c} would be the cutoff distance related to each ‘fluid particle’), for example.
The twoway coupling scheme is capable of dealing with 3D systems containing a large number of particles (N ∼ 1000), but it can take a considerable computation time. In addition, through this section it has not been addressed other factors that could influence the particle dynamics and that may increase the required computational effort. Examples of these factors can be confining wall effects,^{26,27 } complex flow geometries,^{28 } nonspherical particles,^{29–31 } rotational degrees of freedom,^{32,33 } polydispersity^{34–36 } or complex nonhydrodynamic interactions that do not have to be pairwise.^{37–39 } Thus, in order to include and study these new features while keeping a reasonable computational time and cost, a final approximation is commonly used. It is the oneway coupling, where f⃗_{c} = 0. This scheme completely disconnects particle and fluid dynamics. Both problems are solved independently: the flow field is solved first without any immersed particle and then it is introduced in the particle dynamics problem through Stokes drag law.
1.2 Magnetic Interactions
In Section 1.1 it has been shown that theoretically it is possible to compute particle trajectories in a CMRF as long as the nonhydrodynamic forces acting over all dispersed particles are known at each time step. In this section only magnetostatic forces are addressed because they are the most relevant ones in magnetorheology. It will be seen that, in order to exactly compute the magnetostatic forces, the particle spatial configuration and how the total magnetic field, flux and magnetization distribute in the whole system are needed.
Ultimately, these distributions will depend on the magnetic behavior of the dispersed particles, and therefore on the particle material and size. Magnetic particles used in the formulation of CMRFs are usually made of (soft) ferro or ferrimagnetic materials such as iron (and its oxidized derivatives: magnetite and maghemite), nickel, cobalt or their alloys. The reason for this is that these materials exhibit a reasonably large magnetization under moderate magnetic fields, i.e. for a field strength of ∼100 kA m^{−1}, magnetization reaches values of the order of 1000 kA m^{−1} (iron), 300 kA m^{−1} (nickel) and 100 kA m^{−1} (cobalt).^{40 }
Regarding their size, particles within a CMRF have a typical diameter of the order of microns (above 1 μm). From an atomistic perspective, these can be seen as multidomain magnetic bodies: each particle can be divided into domains, separated by Bloch walls, where the atomic spins are strongly coupled, all pointing in the same direction. As a consequence, in every domain there is net magnetization. In the absence of an external magnetic field, the direction related to each domain within the particle is randomly distributed yielding a null particle magnetization.
When an external magnetic field H⃗_{ext} is applied (acting as a homogenous background field), domains aligned, or almost aligned, with the field start to grow at the expense of misaligned domains. As a result, the net particle magnetization becomes nonzero. Then, as H⃗_{ext} is further increased, favorable domains grow larger and magnetization increases till all atomic spins are aligned with H⃗_{ext}. At this moment the particle is fully saturated, that is, increasing H⃗_{ext} does not result in a higher magnetization level.^{41 }
The atomistic description of the particle magnetization mechanism is an extraordinarily complex task that is out of the scope of this chapter. Fortunately, this is not necessary as all these microscopic events can be accounted for when working with ‘macroscopic’ (understood at the particle scale) magnitudes. In this sense, the previous discussion can be simplified defining the particle magnetization as a function of the total magnetic field evaluated inside the particle (the socalled internal field H⃗_{int}) with some phenomenological features:^{42 } (i) it should saturate to a constant value M_{S} at high magnetic fields, (ii) it should increase linearly (slope of μ_{i} − 1) at low magnetic fields, (iii) it should be zero in the absence of the field and (iv) it should be collinear with the magnetic field. Note that the previous description captures a nonlinear response and assumes the absence of coercivity or anisotropy.
There are several models for the magnetizationfield constitutive equation. In the literature, the Fröhlich–Kennelly model is commonly chosen:
where μ_{r}(H_{int}) is the scalar and fielddependent relative permeability while M_{S} and μ_{i} (saturation magnetization and initial magnetic permeability, respectively) are the two parameters of the model. It must be noted again that H⃗_{int} stands for the total magnetic field inside the particle and this one does not meet with H⃗_{ext}.
Generally speaking, if a finite magnetic body is magnetized, its induced magnetization will act as a source of a secondary magnetic field H⃗_{sec} defined over the whole space, inside and outside the body:^{41 }
where the first (second) integral is done over the body volume (surface), R⃗ is position vector of the point where H⃗_{sec} is computed, relative to the considered differential element of the magnetic body during integration; ρ_{M} = −∇·M⃗ is the magnetization charge density, ϱ_{M} = n̂·M⃗ is the magnetization charge surface density and n̂ is the normal vector of the body surface (pointing outwards).
Clearly, H⃗_{sec} depends on its source magnetization (how it is distributed in the magnetic body), on the body shape as well as on the evaluation point (it is not homogeneous). Traditionally, when H⃗_{sec} is evaluated inside its own source body, it becomes to be known as the demagnetization field H⃗_{de} since it has opposite directionality to its source magnetization. In the following we will keep this convention, saving the name H⃗_{sec} only for the secondary field created by the body outside it.
Thus, for the case of an isolated body, it is easy to identify the only two contributions to the internal field as H⃗_{int} = H⃗_{ext} + H⃗_{de}. In Figure 1.2, the relative permeability for an isolated magnetizable spherical particle (H⃗_{de} = −M⃗/3) is shown as a function of H_{int} and H_{ext}. As it can be seen, it undergoes three differentiated regimes (the boundaries indicated by illustrative vertical lines): a linear regime at low fields where μ_{r}(H_{int}) → μ_{i}, a saturation regime at high magnetic fields where μ_{r}(H_{int}) → 1 and a transition regime connecting previous values, which is responsible for the nonlinear behavior.
Results shown in Figure 1.2 can be computed easily because only an isolated particle has been regarded. On its part, CMRFs are constituted of a large number of magnetized particles, therefore, a particle β in a given CMRF will be exposed to a local field that is the sum of H⃗_{ext} plus all the nonhomogenous secondary fields created by the rest of the magnetized particles H⃗_{sec,α}: . Bearing this in mind, the internal field of particle β is given now, as a function of the local field instead of the external one, by H⃗_{int,β} = H⃗_{loc} + H⃗_{de,β} with H⃗_{loc} evaluated at the particle β position. At this point two main problems are present. Firstly, the calculation of H⃗_{sec,α} due to the rest of dispersed particles is a hard task as it depends on the corresponding H⃗_{int,α} that, at the same time, depends on H⃗_{loc} but now evaluated on the rest of α particle positions and thus, containing the initially pursued H⃗_{int,β}. Secondly and more severe, as none of the CMRF particles are under a homogeneous field, their demagnetization fields will not fulfill H⃗_{de} = −M⃗/3 as in the case of an isolated particle. Indeed, the only thing that can be said is that the relationship between both magnitudes has to fulfill eqn (1.25) when it is evaluated inside each particle.
As can be seen, the computation of the magnetic field inside a CMRF is not trivial and requires a selfconsistent approach. Nevertheless, from first principles, the total field H⃗ at any point can be obtained from Maxwell equations:^{44 }
where D⃗ and B⃗ are the electric and magnetic field flux densities respectively, E⃗ is the electric field, ρ is the free charge density and j⃗ is the free current density. As it is known, fields and flux densities are not independent from each other, but they are related through the constitutive equations:
being ε_{0} and μ_{0} the vacuum permittivity and permeability respectively, P⃗ the material polarization and M⃗ the previously introduced material magnetization (for example, see eqn (1.24)). Typically, in CMRFs, the field changes in time are very slow allowing a quasistatic study. In addition, there are not free charges nor currents thus, eqn (1.26)–(1.29) are reduced to:
where V_{e} and V_{m} are scalar electric and magnetic potentials. As can be seen E⃗ is fully decoupled from H⃗, and, since the former does not play any relevant role, it will not be regarded in the following discussion. However, before discarding it, it is worth noting that the total analogy between E⃗ and H⃗, pointing that phenomena observed in a magnetic suspension must have an electric counterpart based on polarizable particles. Effectively those systems exist and are known as electrorheological fluids (ERFs). Although, from an experimental point of view, ERFs entail some crucial differences with respect to CMRFs (which made them less useful from an application point of view), ERFs can be theoretically seen as electric analogues when CMRFs are interrogated in their linear magnetic regime.^{45 }
Eqn (1.32) and (1.35) are applicable in the whole CMRF, that is, inside each dispersed particle and in the continuous phase. Just at the interphase between the two media, continuity equations apply:
where n̂ is the normal vector of the interphase and the subscripts p and c denote field/flux density in the particulate and continuous phase, respectively. Eqn (1.30)–(1.37) together with an appropriate boundary condition (for example, H⃗ = H⃗_{ext} far away from the CMRF if the CMRF is of finite size) are a closed system that allows us to solve the magnetostatic problem yielding magnetic field, magnetic field flux density and magnetization.
Once the field variables are known, the magnetic interaction, i.e. the force acting between magnetized particles, can be calculated. To do so, one can integrate the force density given by the wellknown expressions from Lorentz or Kelvin^{46 } over the particle volumes. However, for the typical methods used to solve the magnetostatic problem in CMRFs (see below), this kind of integration can imply some numerical issues related to the field gradients in the particle volumes and field discontinuities at particle surfaces.
At this point, it is more advantageous to compute the force following Maxwell's stress tensor formulation. This approach directly flows from the Lorentz force approach and expresses the force in terms of the total magnetic field and flux density (instead of their gradients) that are the most easily accessible quantities. In addition, it reduces the complexity of the computation since it is based on the integration of the proper tensor components over an arbitrary surface that must encompass the body over which force is to be computed.
Surface arbitrariness can be leveraged to choose a simple one, for example, lying completely in the continuous phase (which is a linear material without magnetic response) and not in the body surface (where field discontinuities could appear). Following Maxwell's stress tensor formulation, the magnetic force acting on the body F⃗_{m} is given by:
where ê denotes the direction of the force component that is to be computed, n̂ is the normal vector, pointing outwards, of the arbitrary surface S encompassing the body, and T̃ is Maxwell's stress tensor (do not confuse this secondorder tensor with the torque vector T⃗):
where B and H are the modulus of magnetic field flux density and magnetic field respectively. The use of this definition for T̃ implies that surface integration must be done in a linear medium.^{47 } In addition, that surface must not cross an interphase between media with different constitutive equations.
Finally, in some cases, it could also be useful to compute the energy of the CMRF as a whole in a volume V. In this case, the magnetostatic energy is defined as:^{48 }
The previous expression depends on the constitutive equation and therefore the integrand changes depending on whether dv is placed in the continuous phase or in the particulate one. If both phases are isotropic (that is B⃗ and H⃗ are collinear) the constitutive equations can be expressed as B⃗ = μ_{0}μ_{r}(H)H⃗, and therefore eqn (1.40) is rewritten as:
In the case of a linear media, μ_{r} ≠ μ_{r}(H) and the magnetostatic energy density is reduced to μ_{0}μ_{r}H^{2}/2 = BH/2. For nonlinear materials a more complex dependency W_{V}(H) is expected.
The energy computation is of interest by itself because it allows a comparison of the different particle configurations within the CMRF, discerning which is the most favorable (i.e. minimum energy state) and, hence, the most probable state at equilibrium. In addition, for conservative systems (note that according to eqn (1.35) this is the case), it also offers another way to compute forces by properly changing the configuration of the system:
being e the generalized coordinate related to the desired force component.
As can be seen, once the magnetostatic problem is solved, magnetic field and flux density are known at each point of the space so forces and/or energy can be calculated. Unfortunately, the necessity of applying eqn (1.36) and (1.37) at each interphase restricts the analytical solution to very simple problems consisting of two or three particles in high symmetry configurations. In cases where a large number of particles are taken into consideration, methods based on capacitance matrices,^{49 } multipole expansion^{50 } or boundary elements^{38 } have been applied in the literature. Nevertheless, the aforementioned solutions are restricted to linear materials (that is, not applicable to ferro or ferrimagnetic materials of interest in magnetorheology), and they do not yield a complete solution as this is given in terms of a finite number of multipoles or approximate particle positions at large interparticle distances.
When the materials involved are magnetically nonlinear, numerical methods are practically mandatory. Clearly, the Finite Element Method (FEM) is the preferred choice in magnetorheology.^{43,51–54 } This kind of numerical simulation aims to solve eqn (1.34)–(1.37) directly, including all multipoles, so they can be seen as virtually exact. FEM is a meshbased method, therefore it can become computationally expensive for systems involving very different length scales (as it was noted when the flow field problem was regarded). Due to this fact, FEM is not the definitive method to solve magnetostatic problems related to CMRFs.
1.2.1 Mean Magnetization Approximation (MMA)
Alternatively, one can decide to solve eqn (1.34)–(1.37) up to a first order approximation using the socalled ‘Mean Magnetization Approximation’ (MMA). This is based on a wellknown magnetostatic problem: an isolated spherical particle in a homogeneous external field H⃗_{ext} is magnetized homogeneously (i.e. its magnetization is constant through the particle volume) and creates, outside the particle, the field corresponding to a point dipole placed at the particle center:^{55 }
where R⃗ is the vector joining the particle center with the field point (i.e. the place where the magnetic field is evaluated) and m⃗ is the dipole moment. Since the particle is homogeneously magnetized, this is simply given by m⃗ = V_{p}M⃗ being M⃗ the magnetization:
where β(H_{int}) is the contrast factor, μ_{r,c} is the relative permeability of the continuous phase (required to be linear) and μ_{r,p} is the relative permeability of the particulate phase (linear or not). Taking this into consideration, the MMA states that, although any particle is surrounded by many others in a CMRF, they magnetize following eqn (1.44) but are evaluated at the local magnetic field H⃗_{ext} → H⃗_{loc}. Note that this approximation tries to incorporate multibody effects, but, at the end of the day, it is equivalent to supposing that the particles are still isolated, not under H⃗_{ext} but under another homogeneous field with the value H⃗_{loc}. In addition, eqn (1.44) shows that the magnetization mechanism comes from the permeability difference between the particulate and continuous phases. In particular, it demonstrates that nonmagnetic particles dispersed in a linear magnetic media also interact, at least, through dipoles. This allows one to consider inverse ferrofluids, dispersions of nonmagnetic particles in a ferrofluid, (where β < 0) as analogues of CMRFs (β > 0) whenever the ferrofluid used is linear.
The main advantage of the MMA is that it affords easy computation of the local field at any point in the CMRF because the secondary field due to any dispersed particle is simply H⃗_{sec,α} = H⃗_{d,α} (that now is analytically known through eqn (1.43)):
The summation in eqn (1.45) should run over all particles in the CMRF, however, since the magnetic field decreases with distance, a cutoff distance is commonly used to save computational cost. By doing this, only neighboring particles are considered.
Once the magnetic field is known at every point, the interparticle force F⃗_{d} can be computed using Maxwell's stress tensor, Lorentz or Kelvin's law. However, the force on a (field) point dipole m⃗_{α} due to a magnetic field is also known analytically. What is more, as the particles are replaced by point dipoles, any of them will sustain magnetic torque T⃗_{m} under the presence of the field that tries to align with it. Both magnitudes read as follows:
Nevertheless, as noted at the beginning of this section, the particle magnetization is supposed to be collinear with the field, thus CMRF particles are not expected to undergo this torque. Regarding the force, eqn (1.46) shows that it is zero unless the magnetic field changes in the space. Ideally, CMRFs are subjected to external homogeneous fields, thus the only magnetic force their particles experience should be due to the secondary field created by another (source) dipole m⃗_{β}. In that case, the force can be written as:
where R⃗_{αβ} is the vector joining the dipole m⃗_{β} with the dipole m⃗_{α}. Note that, in general, m⃗_{α} and m⃗_{β} do not have to be the same in magnitude nor direction. These will depend on the local field at the respective positions of both dipoles.
To get a deeper insight into the interparticle dipolar interaction, it is worthwhile particularizing eqn (1.48) for the case of two identical dipoles m⃗_{α} = m⃗_{β} = m⃗, writing eqn (1.48) in spherical coordinates (centered in the source dipole and supposed to be directed in the z direction):
It can be seen that the dipolar force is not central. It does not depend on the azimuthal coordinate but it is anisotropic; repulsive for 55° ≤ θ ≤ 125° and attractive in other cases. The polar component is negative for θ ≤ 90° (thus, it tries to move the field dipole over the source) and positive for 90° ≤ θ (in this case it moves the field dipole under the source). Only at θ = 90° is the polar component null but unstable, that is, any small perturbation that pushes apart the field particle from θ = 90° will provoke the field dipole to align with the source. At θ = 0° and θ = 180° the polar coordinate is also zero, however, in these cases, the position is stable: any perturbation at these positions will be balanced by a restoringlike force that will place the field dipole aligned with the source one. In addition, for a chosen interparticle distance R, it can be seen that the force in the dipole orientation direction, F⃗_{d}(θ = 0°) and F⃗_{d}(θ = 180°), is typically two times larger than in the normal direction, F⃗_{d}(θ = 90°). Thus, when the particles are freely suspended in the continuous phase, this behavior is translated in the formation of chainlike structures directed along the dipole (field) direction.
Finally, it is worthwhile remarking that the MMA predicts a magnetic force scale as:
Looking at Figure 1.2 and eqn (1.44), it can be inferred that, at low magnetic fields when the material behaves linearly with β constant, this force will be proportional to H_{ext}^{2}. On the contrary, at high magnetic fields when the material saturates and M^{2} tends to M_{S}^{2}, the force becomes independent of the external field.
Although the MMA offers a simple solution to the magnetostatic problem, it should be remembered that, strictly speaking, it is only valid in two limiting cases: (i) when the particles are far enough from each other to feel only a homogeneous field or (ii) when particles are fully saturated so that their magnetization is uniform through the particle volume. Nevertheless, case (i) is impossible to fulfill in CMRFs; as is pointed out in eqn (1.49), magnetic forces can be attractive and therefore, even at very low concentrations, they will tend to form particles aggregates where interparticle distances are small. As a consequence, particles will not feel a homogeneous field and the MMA would not be valid. For its part, the saturation regime in case (ii) is hard to reach experimentally in a controlled way, i.e. guaranteeing that all particles experience the same high external field, thus all of them are magnetized to the same (saturation) level.
Undoubtedly, the main drawback of using the MMA as a first order solution is that it does not take into consideration magnetization nonuniformity in the particle volume. This is, ultimately, responsible for the multipolar magnetic interactions. These multipoles enhance the magnetic field level in small regions between particles which leads to greater induced magnetization levels. Consequently, when the interparticle distance is small, multipoles play a critical role in the magnetic force computation. For example, in Figure 1.3 the magnetostatic force in an infinitely long isolated chain of particles is plotted as a function of the chain shear strain in affine deformation for a range of μ_{r,p}. Lines stand for the MMA solution including only dipoles and points stand for the complete solution, all multipoles included, computed by FEM. It can be seen that only when the magnetic particles are separated enough (high strains) or μ_{r,p} is small^{56 } (as in inverse ferrofluids and magnetic iron oxidedoped latex suspensions), the MMA can reproduce the complete solution. In other cases, the true magnetic interaction will be highly underestimated.
1.3 Other Nonhydrodynamic Interactions
As in any suspension, particles in a CMRF experience other nonhydrodynamic forces besides magnetostatic ones. Most of them are classic interactions in colloidal science, so we refer to wellknown treatises in this field for further details.^{57,58 } In this section we only summarize them to see whether they do or do not affect the particle dynamics when compared to magnetostatic forces.
Regardless of the kind of particulate or continuous phase, any pair of dispersed particles will interact through the London–van der Waals or dispersion forces. These are attractive interactions coming from the polarization of a particle atom induced by the atoms of the rest of the particles. For the case of two perfect spherical particles, this attractive interaction diverges as the surfacetosurface distance moves to zero pointing to an irreversible aggregation. Again, particle roughness poses a limit to the attractive force. For the case of CMRFs where d > 1 µm, supposing a typical roughness of the order of nanometers, the ratio of van der Waals to magnetostatic forces takes values of 10^{−5} or smaller.^{59 } As a result, van der Waals forces can be safely neglected when accounting for other nonhydrodynamic forces.
The fact that van der Waals forces favor particle aggregation should not be a problem in the presence of a magnetic field since this is precisely the purpose of applying it, i.e. to promote the directed selfassembly with the difference that van der Waals forces do not yield anisotropic structures. However, in the absence of a magnetic field it is desired that the CMRFs keep kinetically stable and with this aim in mind several routes are followed in the literature: core–shell structures, surfactant addition and surface charges among others. These elements introduce other interactions between the particles (electric doublelayer, osmotic and elastic potentials, depletion, etc.) that need to be considered when studying the mechanical behavior of the system. Their effects depend on a large number of factors (continuous phase permeability, electrolyte concentration, surfactant type and concentration, hydrophobicity of the continuous phase, particle size, etc.); hence, it is not easy to summarize their final role on the rheology of the MRFs.
Here, we will limit this to saying that because the main objective of these approaches is to counteract van der Waals forces, it is expected that an optimal route gives rise to a stabilization force (regardless of its origin) whose order of magnitude is similar to the van der Waals one. Thus, as a first approximation, its effects can be also neglected in the presence of a field.
Classical colloidal forces are due to particle–particle interactions but in a CMRF particles are subjected to body forces as well. The main one is the gravitational force that, including buoyancy effects, can be written as F⃗_{g} = g⃗(ρ_{p} − ρ_{c})V_{p}. The magnetostatic force is typically three orders of magnitude larger than the gravitational one so, a priori, the latter could be neglected. Nevertheless, in the case that the height of the macroscopic sample is above a critical size ∼2dF_{mag}/F_{g}, experimental and simulation studies have shown that gravity can collapse the field induced structures.^{60 }
In addition to previously discussed forces, any particle dispersed in a continuous phase will experience Brownian motion. This is a diffusive movement due to the continuous and random collisions between the dispersed particles and the continuous phase molecules. In order to see its importance in comparison to other forces, a typical force scale related to the Brownian motion can be defined as F_{B} ∼ k_{B}T/d where k_{B} is Boltzmann's constant (1.381 × 10^{−23} J K^{−1}) and T is the absolute temperature (do not confuse this with the torque vector T⃗ nor Maxwell's stress tensor T̃). Taking this scale into consideration, the specialized literature defines the socalled λ ratio (or coupling factor) to measure the relative importance of magnetostatics over thermal forces:
In CMRFs this ratio is very large, above 100, even for low external fields which leads one to discard Brownian motion in these systems.
As it can be inferred from the previous discussion, with the exception of hydrodynamic forces, any classical colloidal forces can be safely neglected in the study of CMRFs because the large particle size (above 1 µm) favors magnetostatic interactions. Indeed, strictly speaking, CMRFs are not true colloidal systems due to the absence of Brownian motion; they seem to share more features with granular suspensions, especially when high concentrations are considered (above volume fractions of 0.4). In the concentrated regime, particles suffer repeated collisions with each other and this makes contact forces (repulsive and frictional ones) become a critical agent, even more important than longranged hydrodynamics, in the macroscopic sample behavior.^{61 }
Different expressions have been proposed in the literature to model contact interactions, although current works in granular rheology usually choose Hertzian repulsion and Coulomblike friction forces for repulsion and frictional forces, respectively. These forces do not introduce a new force scale into the problem to compete against magnetic or hydrodynamic forces. On the contrary, they are conceived to model a hard sphere potential, thus they should be able to bear any applied load.^{61,62 } Ascertaining the actual role of friction in rheology is an open field of research nowadays, even in nonmagnetic suspensions, due to the entailed complexity from all experimental, simulation and theoretical points of view. Nevertheless, the parallelism between these granular suspensions and CMRFs would suggest that contact forces are also an ingredient to be taken into account in the complete description of CMRFs.
1.4 Rheology
Sections 1.1 and 1.2 deal with CMRFs from a microscopic point of view. By solving the governing equations, under given macroscopic constraints (particle concentration, applied flow field, external magnetic field) and initial conditions, it is possible to precisely know the dynamic state of the system, i.e. particle positions, velocities and fluid flow field at any time.
Although complete and detailed, this microscopic description does not inform, by itself, one about the collective behavior of the system, that is, about its macroscopic state where the CMRF is not seen as a suspension of magnetizable solid particles in a Newtonian fluid but as a continuous medium with macroscopic physical properties. To get these macroscopic properties and to relate them to the aforementioned constraints it is necessary to take a further step linking both micro and macrodescriptions.
Similar to other suspensions and colloids, CMRFs are constituted by an extraordinarily large number of interacting particles. As a result, the link between both kinds of descriptions must be done following Statistical Physics formalism. According to this, for each set of macroscopic constraints there is a group (called an ensemble) of compatible microscopic states, each of which having a different probability to really happen.
Eventually, any physical observable x depends on its microconfiguration, however, the macroscopically accessible value is given by its ensemble average 〈x〉: the sum of its value at each microscopic state of the ensemble weighted by the probability of the corresponding microscopic state. This average would imply knowledge of all compatible microscopic states together with their probability, and this is extremely complex. Hence, what is actually done is to average the pursued physical observable spatially over the system:^{1 }
This procedure is valid whenever all particles in the system are statistically equivalent, that is, each particle is surrounded by a configuration with the same macroscopic properties. In this way, each particle can be regarded as placed in a different microscopic state, but compatible with the imposed macroscopic one, which allows identification of a volume average with the ensemble average.
Since the mechanical state of a material is determined by its total stress tensor 〈t̃〉, this tensor together with 〈Ẽ〉 will be the physical observables that will be ascertained macroscopically:^{1 }
In eqn (1.53) it has been taken into account that the total volume occupied by the CMRF can be split between the volume occupied by the fluid V_{f} and the volume of each immersed particle V_{p}. In addition, since the carrier fluid is a Newtonian liquid, the total stress tensor within the first integral has been substituted by Newton's viscosity law (eqn (1.6)). Extending the first integral to the whole volume V, remembering that particles are rigid (i.e. Ẽ = 0 in V_{p}) and bearing in mind that 〈Ẽ〉 = Ẽ_{∞} (the macroscopically imposed flow field) we arrive at:
The stress tensor t̃ + pδ̃ inside the rigid particles is undetermined. However, supposing that there is not any nonhydrodynamic torque acting over the particles (which is fulfilled in CMRFs if magnetization is assumed to be collinear to the field), it can be written using the stresslet S̃ (eqn (1.8) and (1.9)) over the particles as:^{17,63 }
where S_{p} is the particle surface and ϕ = NV_{p}/V is the particle volume fraction. Eqn (1.55) is general, providing that Stokes flow without nonhydrodynamic torques are accepted, because no assumption has been done regarding the volume fraction. It indicates that the particle contribution to the macroscopic stress is just the sum, per unit volume, of the stresslets of all particles plus an isotropic term ∝〈δ̃*〉. The latter can be merged with the isotropic mean pressure −〈p〉δ̃ of eqn (1.54) since they do not play any role in the rheology of incompressible magnetic colloids.
The value of S̃_{α} can be computed from its definition (eqn (1.8) and (1.9)) once the flow field around each particle is known. However, remembering the theoretical development of Section 1.1.1, stresslets can be directly solved from the linear system presented in eqn (1.20). In particular, we get:
Remember that bold symbols (in the following omitted) stand for the stresslet associated to all suspended particles. It can be seen that the total stresslet has two contributions, the socalled hydrodynamic stresslet S̃_{h} (all those terms accompanying Ẽ_{∞}) and particle interaction stresslet S̃_{p} (proportional to F_{nh}). It is important to note here that in order to avoid a long expression, rotation contributions in eqn (1.56) have been omitted. However, since in CMRFs nonhydrodynamic torques are supposed to be absent, full expression (i.e. also containing the rotation contribution) will be still separable only in terms that exclusively accompany Ẽ_{∞} or F⃗_{nh}. At this point, some issues can be highlighted.
Firstly, the hydrodynamic stresslet S̃_{h} is not null even when F⃗_{nh} = 0 (i.e. when particles do not interact through nonhydrodynamic forces). The origin of this contribution lies on the infinite resistance of rigid particles to be strained, thus the only presence of rigid particles in the straining flow field contributes to S̃_{h}.
Secondly, the particle interaction stresslet S̃_{p} (in the following particle stresslet only, although both S̃_{h} and S̃_{p} come from the dispersed particles) clearly comes from nonhydrodynamic forces F⃗_{nh} experienced by the particles. However, it should be kept in mind that the origin of any stresslet is the distribution of the traction force t̃·dS⃗ over the particle surface. This indicates that S̃_{p} exists, not only because there exist nonhydrodynamic forces, but also because there is a continuous phase able to mediate the nonhydrodynamic interactions and transmit their effects. In other words, S̃_{p} would be zero even though F⃗_{nh} ≠ 0 in a static ensemble of particles (fixed in the space) able to interact through longrange forces (for example, magnetic ones) but not dispersed in a fluid. Furthermore, if we suppose the same ensemble but now dispersed in a quiescent fluid (for example, an MRF under an external magnetic field without a flow field), since u⃗ − V⃗_{∞} and Ẽ_{∞} are both zero, there would not be any stresslet either.
The two previous comments demonstrate that S̃_{p} cannot be the only contribution coming from nonhydrodynamic forces to the total stress tensor 〈t̃〉. By definition, 〈t̃〉·dS⃗ must account for the total force acting across dS⃗. If we suppose an imaginary surface halving the ensemble in the previous example consisting of longrange interacting particles not immersed in a fluid, all particles in one half would act over the other half with a net force, but this would not be taken into account since S̃_{p} = 0. To account properly for the total nonhydrodynamic force experienced by the particle in the total stress tensor, the socalled term^{64 } or thermodynamic stress must be included in the computation of 〈t̃〉, where:
wherein r⃗_{α} is the position vector of particle α and F⃗_{nh,α} is the total nonhydrodynamic force acting over it. Although the example proposed here only considers longranged nonhydrodynamic forces, eqn (1.57) is valid for hardsphere interactions (and, for example, contact forces) as well.
Putting together all previous contributions, the rheological constitutive equation of the suspension can be written as follows:
Here 〈p〉δ̃ is the mean pressure corrected by the isotropic term 〈δ̃*〉 derived during the stresslet computation (see eqn (1.55)). At first sight, eqn (1.58) already shows some characteristics that differentiate 〈t̃〉 from the Newtonian constitutive relationship:
Under the Stokes flow approximation, all contributions to 〈t̃〉 depend on the particle positions. Thus, in addition to interparticle forces, their structure also plays a role in the rheological behavior.
In the absence of flow and depending on the particle positions, due to the term, the system can develop shear stresses (nondiagonal elements). Thus, when an external stress is imposed trying to move the system from its equilibrium state, an internal restoring stress can appear to keep particles in their original configuration, i.e. there is a yield stress below that the system does not flow.
Not all terms are directly proportional to Ẽ_{∞}. On the contrary, some of them are proportional to Ẽ_{∞} while others are scalable by F_{nh}. These two sources of stress imply different stress regimes depending on which source is prevailing. As a main result, CMRFs show a marked shearthinning behavior.
Terms coming from particles and their interactions, , 〈S̃_{h}〉 and 〈S̃_{p}〉, are not isotropic and consequently normal stress differences can appear at rest or in simple shear flow.
In Sections 1.4.1–1.4.4, these topics are reviewed from theoretical, experimental and numerical points of view. When the system is strained, the discussion will be mainly focus on simple shear as it is the most studied flow mode in MRF research. In this case, the flow direction is supposed to be ŷ, the vorticity direction x̂ and the gradient direction, together with the external field direction ẑ (see Figure 1.4). With this choice, the flow field is v⃗_{∞} = zŷ, the strain rate tensor is:
and is the shear rate. Also, to simplify notation, the stress in the shear direction will be denoted as 〈t〉_{zy} = τ while deviatoric normal stresses (components of total stress tensor in the diagonal excluding the isotropic pressure contribution) will be denoted as σ_{xx}, σ_{yy} and σ_{zz}.
1.4.1 Structure
As indicated in Section 1.2.1, the magnetostatic interaction is anisotropic and contributes to the particle aggregation in the applied magnetic field direction. Traditionally, the fieldinduced structuration of an MRF is studied under no flow conditions in connection to the onset of a yield stress and with emphasis on determining the time scale for structure formation and associated order parameters.
In view of the theoretical analysis in Sections 1.1 and 1.2, imposing that the background flow is zero, the particle positions can be tracked by solving eqn (1.20). Since there is no flow and particle inertia is neglected, there is only one time scale in the structuration problem. As a first approximation, this time scale can be obtained as the time needed by a particle to displace a length equal to its diameter under the action of dipolar magnetostatic interaction (eqn (1.50)) and Stokes’ drag force (eqn (1.21)):
As observed in eqn (1.60), the time scale depends on the applied magnetic field strength (contained in the magnetization) and the carrier fluid viscosity. In view of this analysis, the final structure morphology of an MRF is independent of these magnitudes; t_{a} simply determines how quickly the structuration process occurs.^{65 }
Numerical studies (using magnetic dipolar forces, the Stokes’ drag approximation and a oneway coupling scheme without a background fluid velocity)^{66–68 } together with experiments (based on optical microscopy observations for 2D systems^{69–71 } or light scattering^{72–74 } and computerized tomography^{75 } for 3D systems) show that the final structure morphology and inner particle arrangement strongly depend on the volume fraction, sample confinement and way of applying the magnetic field.
When a uniaxial DC magnetic field is suddenly applied and then kept constant, particles reorganize and start to form doublets that quickly merge into short chains along the field direction. The onset of this doublet formation also depends on the particle concentration: the larger the concentration the smaller the initial interparticle distance and thus the smaller onset time.
In a second stage, these chains start to aggregate into more complex structures depending on the volume fraction and the confinement. Nevertheless, as the chain length increases, the aggregation rate is reduced due to the higher hydrodynamic resistance (the longer the chains the larger the resistance) and also the nature of the magnetic interaction between chains. The latter can be repulsive or attractive depending on chain length (which is determined by the confinement) and the mean interchain distance (which is controlled by the volume fraction). Generally speaking, low volume fractions favor isolated chains that repel each other giving rise to a final structure consisting of simple chains of one particle width. As volume fraction increases (above 0.05) the initial chains have a richer configuration spectrum (shifted or not with respect to each other in the confining direction, with their ends at different z positions, different chain lengths) that leads to the formation of thick interconnected columns/fibers.
The aggregation process described in the paragraph above is very rapid; its time scale is of the order of milliseconds for the chain formation and of the order of minutes for the subsequent fiber formation. However, these final aggregates are not equilibrium structures (in the case of interacting dipoles, the expected equilibrium structure is a tetragonal body centered lattice BCT^{76–79 }). Instead, the particles are kinetically arrested in open elongated flocs with a large number of defects. There are several possibilities to avoid the formation of kinetically arrested states, such as slowly increasing the magnetic field (instead of suddenly applying it) or using toggled fields.^{80 } The latter allows the structure to accommodate by diffusion during the short period when the external field is absent. This approach results in the longterm formation of disperse ‘drops’ consisting of particles that are arranged in a BCT lattice.^{81,82 }
1.4.2 Preyield Regime. Static Yield Stress
Undoubtedly, the yield stress τ_{0} is one of the most relevant figures of merit in an MRF.^{83 } Theoretically, it is defined as the critical stress below which a material does not flow but behaves as an elastic solid. Thus, if the applied stress is lower than τ_{0}, the sample will undergo a strain but once the applied stress is removed the sample will recover its initial state. On the contrary, if the applied stress is larger than τ_{0}, it is said that the sample leaves the elastic or preyield regime, breaks (also called sample failure) and starts to flow or to strain continuously.^{84 }
However, the previous definition together with its related phenomenology is rarely seen experimentally (e.g. Wang et al.^{85 }). In terms of the measured shear rate versus the applied stress, it is found that samples experience a drastic change in the shear rate (usually several orders of magnitude) when the applied stress is varied through a narrow interval (smaller than one order of magnitude) around τ_{0}. This narrow interval of stress is called apparent yield stress. Generally speaking, for stresses outside this interval (above or below), the sample shows an almost Newtonian behavior with the shear rate proportional to the applied stress.^{86 } As can be seen, the behavior at stresses above the aforementioned interval perfectly agree with the definition of τ_{0}. On the contrary, the appearance of the Newtonian behavior below τ_{0} is controversial as the MRF is not perfectly elastic.
This discussion and related controversy around τ_{0} measurement or existence are wellknown and not new and highlight the paradox between what flows or not: to properly demonstrate that a material is elastic for any stress below τ_{0}, it would be necessary that experiments of infinite duration were set up to see that the shear rate is effectively null in this regime.^{87,88 } Actually, what can be stated is that an apparent yield stress exists (in the sense of a dramatic change in the rheological properties) while the related τ_{0} is subjected to a time scale that has been previously chosen, which is typically the time scale of the experiment or problem under study. In this way, below the apparent yield stress, the sample can be seen as a solid albeit it creeps at very slow rate, i.e. very slow in comparison to the experimental time scale. Therefore, τ_{0} is a magnitude that depends on the measurement conditions. Nevertheless, its circumspect use, at least, as an engineering concept has been undoubtedly useful in material characterization.^{89 }
Yield stress ambiguity has given rise to several yield stress definitions. In this chapter, we will deal with two of them: the socalled static and dynamic yield stresses. Static yield stress will be discussed in this section and is identified with τ_{0}, that is, the lowest stress required for the onset of irreversible flow.^{89 } Leading on from this, the dynamic yield stress (τ_{y}) will be introduced in Section 1.4.3.
Motivated by torquetransfer applications in electromechanical systems, experimental studies on the yielding behavior of CMRFs have mainly focused on the dependence of the yield stress with the applied magnetic field strength. This dependence has been usually parameterized as a power law τ_{0}∝H$extm$ where the exponent m is not constant but changes with the field strength. At low fields (linear regime) some authors obtain m = 2.^{90–94 } Increasing the magnetic field strength, the exponent spans from 1.4 to 1.7^{94–103 } while further increasing the field strength leads to a (saturation) regime where the field dependence disappears, m = 0.^{98,104–106 } Note that these three intervals correspond to the same ones observed in the magnetic permeability (and thus magnetization) of a particle (Figure 1.2) in Section 1.2. This observation highlights the close relationship between yield stress, interparticle force and particle magnetization.
Apart from the magnetic field strength, the yield stress is also strongly dependent on the particle concentration. However, very few systematic works exist describing this dependence probably because of the complexity involved when dealing with a manybody problem both from an experimental and theoretical point of view.^{107,108 } Generally speaking, for CMRFs at lowmoderate fields, the yield stress increases linearly at low volume fractions (below ϕ = 0.1) and at a slightly faster rate for intermediate volume fractions (ϕ > 0.1).
Previous experimental trends have been investigated from both analytical and simulation points of view as well. Theoretically, hydrodynamics can be neglected in the preyield regime (i.e. absence of flow), therefore the strain suffered by the sample takes a very long time to relax, so that the only nonisotropic contribution to the total stress tensor is given by eqn (1.57). This depends on the nonhydrodynamic forces exerted over the particles and their positions, i.e. their microstructure. Hence, the goal is to find the rheological constitutive equation that relates the shear stress generated by the structure (restoring stress based only in the nonhydrodynamic forces) when it is shear strained an amount γ, τ = f(γ).
A typical example of this kind of curve is plotted in Figure 1.5a where different regimes can be delineated. Starting from the left, at sufficiently low strains, particle positions are a little distorted and consequently the interparticle force change is linear. This results in an elastic behavior where the dependence between the restoring stress and the applied strain is also linear (its slope is defined as the shear modulus G): as the structure is strained from the (meta)stable state, γ = 0, to γ > 0, a larger restoring stress appears and tries to turn the structure back to the original state. Then, the stress stops being linear, reaches a maximum τ_{0} = f(γ_{0}) and starts to decrease: in this second stage, greater strains find smaller and smaller restoring stresses, this implies that the structure cannot remain stable and breaks. In this context, it is easy to identify the aforementioned maximum τ_{0} with our definition of static yield stress understood as the onset of flow. The related strain γ_{0} will be called here the ‘yield strain’ despite the fact that some references save this name for the limit of the linear elastic behavior.^{109 } Analogously to the yield stress, the yield strain denotes the limit where the structure can recover its original shape once the applied stress has been removed.
As it was said in Section 1.3, in CMRFs nonhydrodynamics forces are completely dominated by magnetostatic ones. Thus, the restoring stress will be given by the magnetic field distribution inside the CMRF (Section 1.2). However, as described in that section, the complexity of the fieldinduced structure together with the length scale mismatch between the structure (sample length scale around 1 mm) and its constitutive elements (particles length scale around 1 μm) prevent us from fully solving the magnetic field distribution, and with it the magnetostatic forces, neither theoretical nor computationally. Because of this, two main approximations are assumed in the literature to get some insight on the yield stress. On the one hand, magnetic interactions are supposed to be known, while the structuration under shear is computed. On the other hand, particle positions are imposed and both magnetostatic forces and the total stress are solved.
The first approximation is just the oneway coupling scheme referred in Section 1.1: magnetic forces are supposed to be dipolar (eqn (1.48)) while hydrodynamic forces are reduced to Stokes’ law (eqn (1.21)). From the total force acting over each particle, their positions can be tracked by solving Newton's second law (eqn (1.4)). Different from structuration studies in the quiescent state, now it is necessary to move the particles in order to evaluate how the stress changes with the sample shear strain. To do this, exceedingly small shear rates are applied to drag the particles, and the strain is computed as γ = t (where t is the simulated time). Obviously, this contradicts our definition of preyield regime (where there is no flow or any dependence on time), however, if the drag force coming from the imposed shear rate is small enough in comparison to the magnetostatic force (ratio of 10^{−3} or smaller^{36,110 }), the computed shear stress versus γ still shows the expected behavior consisting of a linearly increase followed by a maximum, that is identified with the static yield stress.
Generally speaking, although this approximation is capable of mimicking the structures experimentally seen in CMRFs, it underestimates the experimental yield stress values, probably because magnetic multipoles are not included in the computation.
In the second approximation, the MRF structure is imposed and the force/stress computed from magnetostatic theory (Section 1.2). Traditionally, models within this group can be classified as being macro and microscopic depending on whether they account for particle positions inside the field induced aggregates or not. In Sections 1.4.2.1 and 1.4.2.2, the main features behind both models are discussed.
1.4.2.1 Macroscopic Models
In the macroscopic approaches the particles are supposed to arrange following the applied field direction in fusiform (cylinders or ellipsoids) or sheetlike aggregates whose internal structure (precise position of the constituent particles) is unknown but described by the aggregate internal volume fraction ϕ_{a} (commonly approximated by the randomly closed packing for spheres ϕ_{a} = 0.64). In this case, the stress curve τ = f(γ) is calculated starting from the magnetostatic energy (eqn (1.40)) that for linear materials can be rewritten as (see pages 95–98 in Rosensweig^{111 }):
where m_{z} = V(1 − μ_{zz}^{−1})H_{ext} is the dipole component in the field direction written in terms of the permeability tensor μ_{zz} of the suspension. This depends on the aggregate shape and strain. In the case of sheet like aggregates and cylinders, it can be written as follows:^{112,113 }
Here μ_{‖} and μ_{⊥} are the permeability components in the parallel and perpendicular directions to the aggregate axis respectively. Both are independent of the field (linear material) and the strain but they depend on the aggregate shape, contrast factor β and sample volume fraction. Rather than focusing in deriving these dependencies (the full expression computed according to the Maxwell Garnett theory can be found in Bossis et al.^{112 }) here they are merged in a general function g to directly study τ. Supposing isolated aggregates that do not interact with each other:
where the dependence of g on the aggregate shape is captured through β and the dependence on γ must be included again to cover the case of ellipsoidal aggregates. Nevertheless, it has to be said that dependences on γ and ϕ are weak. Consequently, the yield stress will happen at the maximum of the function 2γ/(1 +γ^{2})^{2} that is around 0.58.
Two main factors control the (yield) stress. The first is that it is proportional to the external magnetic field squared due to the linear magnetic behavior previously assumed. Second is that despite the fact that aggregates are supposed not to interact with each other, the stress is not a linear function of the volume fraction. The explanation for this is that as ϕ gets closer to ϕ_{a}, the aggregates tend to fill the whole sample volume. Consequently, at ϕ = ϕ_{a} the whole sample is a gapspanning homogeneous aggregate (since internal structure is not taken into account) that does not change its energy under the external field as it is strained (i.e. there is not any difference between an unbounded homogeneous body that is strained or not).
The comparison between the macroscopic models and experimental data on the yield stress of CMRFs is not satisfactory.^{112 } Similar to what happened with the first approximation using the oneway coupling scheme, the theoretical predictions are below the experimental data. Clearly, a good match was not expected at intermediatehigh fields when particles are partially or fully saturated as the model was thought for use in linear media. However, in addition, the agreement is not good at low fields when the particles fulfill this requisite. The main problem is the absence of interactions between neighbor aggregates and of multipolar contributions in the model. The latter will enhance the local field in the particle poles and gaps, hence increasing the magnetic field in these regions. On the contrary, the model assumes that the field inside of the sample is homogeneous and determined by μ_{zz} which cannot approximate the real field distribution of the field in the aggregates. Despite this finding, macroscopic models do agree very well when compared with experiments on inverse ferrofluids.^{107,112 } The fact that ferrofluids have a much smaller permeability than ferromagnetic particles used in CMRFs allows neglecting multipoles (see Figure 1.3) and therefore the macroscopic model properly works because it is not crucial to know the exact position of every particle. The agreement with inverse ferrofluids is good not only regarding the field dependence but also regarding the volume fraction dependence. It can also be seen an experimental maximum in the yield stress with the volume fraction that is also anticipated theoretically by eqn (1.63).^{107 } Besides, it must be said that these macroscopic models have been also derived for the case of nonlinear materials (thus involving a numerical resolution) and grounding the computation of the force on Maxwell's stress tensor.^{114 } In this case, the agreement with experiments on CMRFs is comparable with that provided by microscopic models exposed in Section 1.4.2.2.
Finally, it is also worth mentioning a last approximation for the preyield behavior based on the continuum theory at small strains. In this context, the CMRF is seen as a homogeneous body with a fieldindependent anisotropic permeability tensor but now independent of the aggregate shape. In fact, the model does not assume the existence of aggregates but anisotropy in the magnetic permeability.^{115 } Of course, the preferred direction is that along the fieldinduced aggregates that changes when the body is sheared. By computing how the permeability tensor changes with the strain (something that is accessible experimentally), it is possible to compute the rheological constitutive equation of the material. Although the model is limited to linear materials and small strains, the major advantage behind this formalism is that it does not require a microscopic description of the sample, that is, it does not depend on the particle shape/size, magnetization mechanism (dipolar or multipolar terms) or plausible forces between particles (non/hydrodynamics ones) that played a role in the eventual particle structure.
1.4.2.2 Microscopic Models
In microscopic approaches the position of each particle is imposed for every applied strain. This allows computing of the magnetic field distribution from the Maxwell equations, (eqn (1.34) and (1.35)), from this the magnetic force between particles (eqn (1.38)) for each strain and finally the constitutive relationship τ = f(γ) to identify the yield stress from the maximum of the curve. A priori, any particle configuration could be solved, however, this is a hard task, even for a few particles, due to the necessity of imposing the continuity conditions (eqn (1.36) and (1.37)). Only in the case of linear materials, methods based on multipolar expansions that always include a finite number of multipoles L (Lth order solution) have been implemented.^{49,50 }
In the case of nonlinear materials such as CMRFs, the most frequent choice is to use FEM to numerically solve the Maxwell equations without any further approximation.^{54,116–120 } However, the high computational cost associated to FEM has restricted the study to strong symmetric systems; in particular, periodic arrangements strained under affine shear mode where particles maintain constant their positions in the z and x directions while the y component is changed as y = γz (see Figure 1.5b). The main advantage of those arrangements is that they can approximate the CMRF only computing the magnetic field distribution around one or two particles in a unit cell while the effects of neighbors (rest of dispersed particles in the CMRF) can be accounted for using the appropriate boundary conditions on the magnetic field (typically, mirror symmetries and/or periodic boundary conditions).
In view of the structures experimentally reported (see Section 1.4.1), the natural choice for the periodic arrangement would be the 3D tetragonal lattice with body centered basis that evolves to a monoclinic one when it is sheared. However, only a few works really simulate a 3D lattice.^{54,121–123 } The majority implement axisymmetric structures that hardly resemble the real particle configuration but imply less computational effort.^{43,116,117,119,124–126 }
In striking contrast to macroscopic models, microscopic ones (solved using FEM) can reproduce the experimental dependence of the yield stress on the external magnetic field. This is basically because they can take into consideration the nonlinear magnetization of the particles and the presence of multipoles. Due to the latter, the magnetic field is enhanced in the gaps present between the particles. As a consequence, the particle polar regions close to the gaps reach larger and larger magnetization levels until it starts to saturate. This eventually explains why the interparticle force saturates as well.
Experimental and numerical trends of the yield stress against the external field strength have been also corroborated from an analytical point of view. These theoretical analyses have been applied to a simple system consisting of one isolated chain but introducing the nonlinear magnetic behavior of the particles. In summary, three regions are identified:
where M_{S} is the saturation magnetization of the particles and the proportional constants depend upon the ref. 116 and 127.
As it can be seen, the three reminiscent regimes of the particle magnetization behavior (see Figure 1.2) appear theoretically as well. The frontiers between linear, subquadratic (transition) and saturated regimes are difficult to establish. Indeed, they depend on the constitutive equation and preassembled structure of the particles among others. What is more, the subquadratic dependence is only a particular case because the exponent of H_{ext} continuously changes from 2 to 0 as the particles saturate from their poles to their center.^{128 } Note that, in any regime, the yield stress is always proportional to the volume fraction. This comes from initially supposing isolated noninteracting chains. Assuming a homogeneous distribution of chains, it can be shown that the area per chain w^{2} goes as w^{2}∝ϕ^{−1}, thus τ_{0} = F_{m}/w^{2}∝ϕ.^{129 }
Finally, with the same goal as in the previous analytical derivation, some models based on the MMA (usually neglecting local field computation) for the yield stress of an isolated chain have been proposed as well. In general, the yield stress can be written as:
where M = 3βH_{ext} in the linear regime (see eqn (1.44)) or M = M_{S} in the saturation regime. c_{3} is a proportionality constant that depends on the particular model assumptions: c_{3} = 0.086 for an axisymmetric model,^{127 }c_{3} = 0.114 for a pair or particles,^{130 }c_{3} = 0.123 for a semiinfinite chain^{131 } or c_{3} = 0.137 for an infinite chain.^{132 } More complex models can be also found, for example, considering multipoles for a pair of linear particles,^{133 } or ‘beadrod’ models where the resulting yield stress is dependent on the confinement of the sample.^{131 }
1.4.3 Postyield Regime. Shear Thinning Behavior
When the applied shear stress overcomes the static yield stress, the structure cannot bear it and consequently breaks. At this stage, the sample starts to flow and the shear stress consists of the last four terms in eqn (1.58). This equation can be rearranged as follows:
where the different terms have been grouped depending on whether they are proportional to (first bracket that we will call ‘hydrodynamic’ contribution η_{h}) or proportional to the nonhydrodynamic forces F_{nh} = F_{d} (second bracket or ‘particle’ contribution τ_{p}).
According to the definitions of the stresslets 〈S̃_{h}〉 and 〈S̃_{p}〉 (eqn (1.56)) and the term (eqn (1.57)), it can be seen that neither η_{h} nor τ_{p} depends explicitly on . However, it should be kept in mind that all those terms depend on the particle microstructure and this one is eventually affected by , thus it should be written η_{h} = η_{h}(ϕ,) and τ_{p} = τ_{p}(ϕ,). What is more, taking into account that τ_{p} is proportional to F_{d} and this is, at the same time, proportional to M^{2} (see eqn (1.48)–(1.50)), it can be further written τ_{p} = f(ϕ,)M^{2}. Note that if the dependence on through the microstructure is neglected, CMRFs would follow a plastic Bingham model.
It is well known in the Rheology community that the use of material functions (such as the shear viscosity in a simple shear flow) is more convenient than the use of the stress itself. Simply dividing eqn (1.66) by and bearing in mind the aforementioned dependencies on ϕ, and M^{2} we arrive to:
The ratio M^{2}/ is traditionally introduced through dimensional analysis. It is the socalled Mason number Mn and can be conceived as the ratio between hydrodynamic (drag) and magnetostatic force scales, advection and magneticcontrolled velocity/times scales, etc. In this chapter, the first option is chosen and therefore:
With this and calling η_{∞}(ϕ) the limiting behavior of η_{h}(ϕ,) at large Mn (see below) the CMRF viscosity can be made dimensionless:
where the critical Mason number has been introduced:
In Figure 1.6, the dimensionless viscosity η/η_{∞} for a given ϕ is sketched together with the contributions from the functions η_{h}(ϕ,) and τ_{p}(ϕ,). This figure must be regarded as a qualitative description to roughly catch the viscosity behavior since it is only adapted from SD results for a monolayer of particles with linear polarization behavior.^{17 } Three regions are usually differentiated in Figure 1.6. By analyzing them, it is possible to give the proper physical meaning to the previously introduced functions τ_{p}(ϕ,) and η_{h}(ϕ,).
At low Mn (the borderline being dependent on ϕ), particle motion is expected to be fully governed by the magnetostatic forces since they are much greater than the drag created by . Thus, when the CMRF starts to flow, the initial aggregates induced by the magnetic field will start to rupture and reattach in more relaxed configurations. However, this reconfiguration of the particles takes place almost instantaneously, or, in other words, at a time scale that is much smaller than the time scale related to . Consequently, the microstructure, and with it τ_{p}(ϕ,), are both independent of .^{110 } In addition, according to Figure 1.6, the contribution from η_{h}(ϕ,) (dotted line) is several orders of magnitude smaller than the τ_{p}(ϕ,) contribution and so considered negligible.^{8,17 } These two facts imply that, at low Mn, the CMRF stress level is only a function of the volume fraction, τ_{p}(ϕ,) = τ_{y}(ϕ). Hence, in the magnetostatic regime, the function τ_{y}(ϕ) truly stands for a dynamic yield stress generated by the continuous (and almost instantaneous) rearrangement of the particles.^{89 } From here and remembering that τ_{p}(ϕ,) is proportional to M^{2}, it is straightforward to see that the viscosity will diverge as Mn decreases in the magnetostatic regime: η = τ/ ∼ τ_{y}(ϕ)/ = f(ϕ)M^{2}/∝f(ϕ)/Mn.
On the contrary, in the hydrodynamic regime at large Mn, particles are dragged by the fluid flow without any significant influence from magnetostatics: F_{mag}/ ∝ 1/Mn ∼ 0 (see eqn (1.50) for the definition of F_{mag}). Hence, the suspension behaves as if it consisted of noninteracting spheres. This makes τ_{p}(ϕ,)/ ∼ 0 but also establishes as the only driving parameter in the CMRF dynamics. Bearing in mind that the Stokes flow is linear in its driving parameters, it is expected that 〈S̃_{h}〉 ∝ (regardless the particle microstructure or ϕ). Introducing this in the first bracket of eqn (1.66), it is seen that η_{h} does not depend on in the hydrodynamic regime, thus: η_{h}(ϕ,) = η_{∞}(ϕ) as it was previously anticipated. Because of the fact large Mn values are usually reached by applying large shear rates, η_{∞}(ϕ) is known as the high shear rate viscosity.
Finally, increasing Mn over the transition regime is expected to trigger particles separation from aggregates. As the aggregate size is reduced, they offer less resistance to the flow, thus viscosity decreases as well. Ideally, the typical size of the aggregates continuously goes down until Mn = Mn* when it can be considered that pairs of particles are forming and breaking continuously in the CMRF.^{134 } Note that in this regime, due to the interplay between the particle microstructure and the different stresslets of eqn (1.66), it is not easy to extract any information about η_{h}(ϕ,) or τ_{p}(ϕ,).
To do so, it is necessary to turn to numerical studies that are capable to discern the different contributions to the shear stress (or viscosity) such as SD^{17} or twoway coupling schemes.^{8 } Results from this kind of study were sketched in Figure 1.6. As it can be seen in the transition regime, the particle contribution τ_{p}/(η_{∞}) (dashed line) is roughly proportional to 1/Mn. This means, that its behavior in the magnetostatic regime (that can be predicted analytically) is also applicable at the intermediate Mn.
Regarding the hydrodynamic contribution, a more complicated behavior is observable due to the destruction of the fieldinduced structure. However, the high computational cost required to study this regime (that is, modeling a large number of particles with a complete description of hydrodynamic and magnetostatic interactions) has precluded a deeper and more detailed understanding of it. Currently, what is done in the literature is to extrapolate the behavior of the hydrodynamic contribution away from the hydrodynamic regime, i.e. to assume that η_{h}(ϕ,) = η_{∞}(ϕ) for any Mn. Clearly, this approximation is appropriate in the magnetostatic regime where the hydrodynamic contribution has little effect. However, it may not be valid in the transition regime unless η_{h}(ϕ,) ≪ τ_{p}/(η_{∞}).^{8 }
Under this scenario, eqn (1.69) and (1.70) are greatly simplified since now, only dependencies on ϕ are kept:
Eqn (1.71) and (1.72) show that, under conditions referred to in the previous paragraph, the flow behavior of CMRFs follows the Bingham model and depends only on Mn and ϕ. Furthermore, this last dependence is solely contained in Mn*. According to this, the flow curve of any CMRF has the same shape as the one plotted in Figure 1.6. The only effect of ϕ is to shift the curve (and the frontiers between the three regimes) in the xaxis according to the value of Mn* (see point in Figure 1.6).
Numerous experimental works in magnetic colloids^{107,119,131,135–141 } and, in minor extent, simulation ones^{19,27,33,142 } have tried to collapse viscosity data in a master curve following eqn (1.71). In general, the data do not exactly follow that trend; instead, they can be better fitted as η/η_{∞} = 1 + (Mn*/Mn)^{Δ} with Δ ranging from 2/3 to 1. The origin of this expression is not physical but just a relaxed version of eqn (1.71) to get a better fit. Actually, a perfect agreement is not expected either since the Bingham model arose from the supposition that both η_{∞} and τ_{y} do not depend on and hence involves a certain grade of approximation, especially in the transition regime.
Other experimental sources of error can appear as well. For example, it is impossible to guarantee that nonhydrodynamic forces will be exclusively the magnetostatic ones in any CMRFs. This fact can introduce new factors not accounted for in this section that will probably have a deeper effect in the transition regime (when neither magnetostatics nor hydrodynamics dominate).^{17 } To enhance the fitting in this regime, other macroscopic models such as the Casson one or microstructural viscosity models have been used.^{143 }
1.4.3.1 Critical Mason Number
Eqn (1.71) offers a route to test experimental and simulation data on any MRF subjected to (steady simple shear) flow. However, until now, it has had a very limited predictive utility because nothing has been said about the functional form of Mn*(ϕ). Looking at eqn (1.72) it can be seen that, indeed, both η_{∞}(ϕ) and τ_{y}(ϕ) dependences are needed separately.
The determination of η_{∞}(ϕ) has been extensively addressed in the literature under the context of the calculation of the high shear viscosity of nonBrownian hard sphere suspensions. The solution was provided by Einstein^{1 } for low volume fractions (ϕ < 0.05). In this concentration regime the particles do not interact, even hydrodynamically, and the stresslet S̃_{h} in eqn (1.66) can be identified with the stresslet of an isolated particle (eqn (1.18)), hence:
Here it has been taken into consideration that 〈Ẽ〉 = Ẽ_{∞} (the imposed macroscopic shear rate) and Ẽ_{∞} is given by eqn (1.59). The calculation of η_{∞}(ϕ) for larger concentrations is more complicated. It involves nonnegligible hydrodynamic interactions that require simulation techniques to be solved. As it has already been pointed out, these simulations are highly time/computationalpower consuming and do not always offer results comparable to experimental data due to the further simplifications needed to implement them numerically. Alternatively, several expressions have been derived analytically following mean field theories. In these, each particle is supposed to experience the (averaged) hydrodynamic interactions due to the rest of suspended particles.^{89 } In this chapter, the socalled Quemada expression is used:^{144 }
where ϕ_{max} = 0.64 is the randomly closed packing for spheres where viscosity diverges.
Regarding τ_{y}(ϕ), both theoretical and simulation approaches have been used to predict it depending on the particle concentration. In the two cases, magnetostatic interactions are supposed to be dipolar (eqn (1.48)) while Stokes’ drag law (eqn (1.21)) is used to model particle–fluid interaction.
Similar to those summarized in Section 1.4.2, theoretical models are limited to isolated chains (dilute regime) but now subjected to simple shear flows. Instead of looking for the maximum of the shear stress as a function of the strain (static yield stress), dipolar and drag forces are balanced in the chain to obtain the dynamic yield stress responsible for the chain break. Again, working with an isolated chain implies τ_{y}∝ϕ whereas the high shear rate is simply η_{∞}(ϕ) = η_{c}. Consequently Mn*(ϕ) = c_{4}ϕ where c_{4} is a constant that depends on the particular forces and the specific mechanical stability condition: c_{4} = 5.25,^{131 } c_{4} = 8.49,^{137 } or c_{4} = 8.82.^{145 } At large concentrations, work has been done with oneway coupling schemes resulting also in a linear dependence with the volume fraction, in this case, c_{4} = 0.24.^{129,142,146 }
From an experimental point of view, MRF dominated by magnetostatic dipolar forces (e.g. inverse ferrofluids), regardless of the volume fraction, and CMRFs in the dilute regime (ϕ ≤ 0.05) show an excellent agreement with previous oneway coupling results. On the contrary, CMRFs at higher concentrations, show that Mn*(ϕ) is not linear and almost two orders of magnitude higher than simulations.^{129,142,146 } It is likely that the disagreement can be explained in terms of other nonhydrodynamic colloidal forces and magnetic multipoles present in the experimental samples but not included in the simulations.
Additionally, it has to be borne in mind that this kind of simulation, without coupling between particles and fluid, could introduce another source of error. A priori, only simulations capable of assessing the three stresslets 〈S̃_{h}〉, 〈S̃_{p}〉 and , forming part of the total viscosity, should be used. However, as was previously mentioned, 〈S̃_{h}〉 (through η_{∞}) is the main contributor to the viscosity only at large Mn,^{17 } when it only depends on the volume fraction, and therefore it can be computed analytically from eqn (1.73) or (1.74) depending on the concentration. With regards to , it is properly computed in oneway schemes as nonhydrodynamic forces are properly taken into account as well. Main claims could arise from 〈S̃_{p}〉 that requires the use of twoway methods at the expense of higher computational effort and smaller simulated systems.
As it can be seen, once the dependence Mn*(ϕ) is known, it is possible to predict the steady shear flow behavior of any CMRF (regardless its volume fraction, the applied shear rate and the applied magnetic field) from eqn (1.71), that is, it is a master curve. As a consequence, it could be considered that we have a general constitutive rheological equation for any CMRF. However, it must be taken into consideration that this section is only focused on steady simple shear flows and therefore shear stresses. A true constitutive equation should also describe the rest of the components of the stress tensor (e.g. normal stresses).
1.4.4 Normal Stresses
The appearance of nonisotropic normal stresses, as well as normal stress differences, in simple shear flow together with the shear ratedependent viscosity and the yield stress are classical signatures of nonNewtonian flow behavior. From eqn (1.58) and taking into account the fieldinduced anisotropic structures devised in Section 1.4.1, it is easy to deduce that MR fluids will have also anisotropic normal stress components in the field/gradient σ_{zz}, shear σ_{yy}, and vorticity σ_{xx} directions (see Figure 1.4).
These normal stresses play a very important role in MRF applications. They promote the appearance of macroscopic forces in the field direction acting over the walls confining the sample, an enhanced samplegeometry interaction to avoid wall slip, magnetostriction effects, etc. Unfortunately, the difficulties related to their precise measurement have hindered substantial progress in their understanding and many previous publications in this topic give contradictory results.
To date, experimental work has focused on measuring the normal force (i.e. in the field direction) developed by the sample because it is a magnitude that can be easily accessed experimentally. In all reported cases, a positive normal force is found (i.e. the sample tries to elongate in the field direction) that approximately increases with H_{ext}^{2}. However, consensus does not seem to have been reached regarding the strain/shear rate dependence of the normal force in the pre/postyield regime. In the preyield regime, some authors found that the normal force grows with the strain up to a maximum associated to the yield stress.^{105,147 } Conversely, it has also been reported that the normal force is initially maximum and decreases to a minimum when the yield stress is reached.^{148 }
From a theoretical/simulation point of view, any of the models presented in Section 1.4.2 can be extrapolated to compute the stress in the field direction as it only supposes to redo operations but working with the fielddirection component of the force instead of the sheardirection one. Thus, results share the same properties, σ_{zz}∝μ_{0}M^{2}ϕ, and would reproduce qualitatively the experimental trend within the field. With respect to the dependence on the volume fraction, it cannot even be tested since most of the work deals only with one concentration, which in addition is typically very high (0.3–0.5), questioning the subsequent comparison with these models as they were in the dilute limit.
Furthermore, any of the aforementioned models predict a negative normal force (i.e. the sample squeezes in the field direction) in the nonstrained state because the magnetostatic force between particles in a chain aligned with the field is attractive. It seems that, in order to reproduce experimental trends, it is necessary to include effects related to the finite size of the sample (demagnetization fields^{149 } and Maxwell's stress jump at the boundaries^{148 }), plausible magnetic field gradients in the experimental setup and/or defects in the induced structure.^{52,150,151 }
During flow, very different observations have been reported as well. While some references report a monotonic decrease of the normal force with the shear rate,^{105,152,153 } others describe an initial decrease at low shear rates but followed by an increase that eventually leads to a plateau at high shear rates.^{148 } Simulation studies of normal stresses under flow are scarce as well. Using SD simulations, in a 2D system (areal fraction of 0.4), Bonnecaze et al.^{17 } report (σ_{yy} − σ_{zz})/ < 0 at any Mn. This quantity decreases, in absolute value, with Mn until it vanishes in the hydrodynamic regime.
1.5 Conclusions
Magnetorheological fluids are systems whose microdynamics and macroscopic behavior are dominated mainly by the hydrodynamic and magnetostatic interactions between dispersed particles. Both kinds of forces are multibody and nonlinear which make the analysis of these suspensions inherently complex. Currently, there exist simulation methods able to take into consideration the full physics behind magnetorheological fluids. Nevertheless, due to the large difference between the micro and macroscales for these systems, the direct and combined implementation of such numerical methods is precluded. As a result, most common descriptions to study the rheological behavior are based on oneway schemes (particle and fluid dynamics problems are totally disconnected) and the Mean Magnetization Approximation (dipolar magnetic interactions). These descriptions are able to offer a qualitative agreement with the experimental findings.
Abbreviations
Dyadic product of the vectors t⃗_{n} and r⃗_{S}
 F⃗_{h}
Hydrodynamic force
 F⃗_{S}
Stokes’ drag force
 F⃗_{g}
Gravitational force
 F⃗_{m}
Magnetic force acting on the body
 F⃗_{nh}
Nonhydrodynamic force
 H⃗_{de}
Demagnetization field
 H⃗_{ext}
External magnetic field
 H⃗_{int}
Internal magnetic field
 H⃗_{loc}
Local magnetic field
 H⃗_{sec}
Secondary magnetic field
 R⃗_{αβ}
Vector joining the dipole m⃗_{β} with the dipole m⃗_{α}
 T⃗_{h}
Hydrodynamic torque
 T⃗_{m}
Magnetic torque
 T⃗_{nh}
Nonhydrodynamic torque
 f⃗_{c}
Force exerted by all the immersed particles over the fluid and responsible for the coupling between particle dynamics and fluid flow problems
 f⃗_{ext}
Body force (per unit volume) exerted by any external agent
 r⃗_{S}
Vector joining each point of the particle surface (where v⃗ is evaluated) with the particle center
Vorticity
 B⃗
Magnetic field flux density
 D⃗
Electric field flux density
 Ẽ
Strain rate tensor
 F⃗
Is the total force acting over the particle
 F_{B}
Brownian force scale
 H⃗
Total magnetic field
 Ĩ
Inertia tensor of the particle
 M⃗
Material magnetization
 M_{S}
Saturation magnetization
 P⃗
Material polarization
 R⃗
Vector joining the particle center with the field point (i.e. the place where the magnetic field is evaluated)
 S̃
Stresslet
 S_{p}
Particle surface
 T⃗
Total torque acting over the particle
 T̃
Maxwell's stress tensor
 V_{c}
Fluid unit volume
 V_{e}
Scalar electric potential
 V_{m}
Scalar magnetic potential
 W_{V}
Magnetostatic energy
 ê
Direction of the force
 k_{B}
Boltzmann's constant
 m⃗
Dipole moment
 n̂
Surface normal vector pointing to the fluid
 t̃
Total stress tensor
 u⃗
Translational velocity of the particle
 v⃗
Fluid velocity field
Grand resistance matrix
 δ̃
Secondorder identity tensor
 ε_{0}
Vacuum permittivity
 η_{c}
Viscosity of the continuous phase
 μ_{0}
Vacuum permeability
 μ_{i}
Initial permeability
 μ_{r,c}
Relative permeability of the continuous phase
 μ_{r,p}
Relative permeability of the particulate phase
 μ_{r}(H_{int})
Relative permeability
 ρ_{M} = −∇·M⃗
Magnetization charge density
 ρ_{c}
Fluid density
 ρ_{p}
Particle density
Angular velocity of the particle
 ϱ_{M} = n̂·M⃗
Magnetization charge surface density
 B
Modulus of magnetic field flux density
 H
Modulus of the magnetic field
 N
Number of magnetic monodisperse spherical particles
 S
Arbitrary surface encompassing the body
 T
Absolute temperature
 V
Total volume
 d
Diameter of the magnetic particles
 e
Generalized coordinate
 p
Pressure field
 β
Contrast factor
 λ
Lambda ratio (or coupling factor)
 F⃗_{d}
Magnetic dipolar force
 H⃗_{d}
Magnetic field created by a dipole
 F_{mag}
Magnetic dipolar force scale
 V_{f}
Volume occupied by the continuous phase
 V_{p}
Particle volume
 ϕ
Particle volume fraction
 S̃_{h}
Hydrodynamic stresslet
 S̃_{p}
Particle interaction stresslet
 Ẽ_{∞}
Macroscopically imposed strain rate tensor
Thermodynamic stress

Shear rate
 τ_{zy}
zy component of the total stress tensor
 σ_{xx}
Deviatoric normal stress in the x direction
 σ_{yy}
Deviatoric normal stress in the y direction
 σ_{zz}
Deviatoric normal stress in the z direction
 t_{a}
Aggregation time scale
 τ_{0}
Static yield stress
 τ_{y}
Dynamic yield stress
 γ
Shear strain
 t
Time
 G
Shear modulus
 γ_{0}
Yield strain
 ϕ_{a}
Aggregate internal volume fraction
 μ_{zz}
zz component of the magnetic suspension permeability tensor
 μ_{‖}
Parallel component of the permeability
 μ_{⊥}
Perpendicular component of the permeability
 w^{2}
Area per chain
 η_{h}
Hydrodynamic contribution to the apparent viscosity
 τ_{p}
Particle contribution to the shear stress
 Mn
Mason number
 Mn*
Critical Mason number
 η_{∞}
High shear rate viscosity
 ϕ_{max}
Randomly closed packing for spheres
Acknowledgements
This work was supported by ERDF, FEDER, MICINN AEI PID2019104883GBI00, Junta de Andalucía P18FR2465 and AFQM396UGR20 projects. J. R. Morillas acknowledges FPU14/01576 fellowship and European Union's Horizon 2020 research and innovation program under the Marie SkłodowskaCurie Grant (EFST)H2020MSCAIF2020 (Grant 101030666).