Chapter 1: Introduction to Magnetorheological Fluids Free

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 micro-configuration, 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 micron-sized magnetizable particles in a Newtonian (nonmagnetic) liquid carrier. In these particular MRFs Brownian motion can be neglected.

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 well-defined 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 fluid-mediated 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 fluid-mediated 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 mass-balance (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 momentum-balance (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 energy-balance 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 well-posed, 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 far-field conditions or periodicity (for unbounded media), a prescribed fluid velocity (the so-called no-slip 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 no-slip 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³/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 non-hydrodynamic 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 second-order 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 second-order tensor would read in index notation a_ib_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:¹ 

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 no-slip 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 many-body 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 mesh-based 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).² 

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.

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 mesh-based 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 mesh-based 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, mesh-free 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 mesh-based 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 space-averaged, are able to reproduce the fluid macroscopic behavior. Nevertheless, these mesh-free 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⃗₀ by its Taylor series:

where and it should be understood that the velocity V⃗₀, the shear rate tensor Ẽ₀ and the vorticity are evaluated on r⃗₀. 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⃗₀, a pure straining flow field Ẽ₀·(r⃗ − r⃗₀) 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:¹ 

Eqn (1.13) is commonly also known as the Stokes’ drag law. Several points can be highlighted from eqn (1.13) and (1.14):

1. 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).

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

3. The flow field and particle movement are supposed to be non-inertial. 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.

4. 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 non-hydrodynamic 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 non-hydrodynamic 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 non-hydrodynamic force, V⃗_α = f⃗_α(u⃗_α,F⃗_nh,α,…).

However, if particle α is also rotating at ω⃗_α and experiencing a non-hydrodynamic 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 ¹²  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.¹³ 

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 two-particle interactions only. Nevertheless, it is stated that during computation, these elements are combined, and thus many-body interactions are finally considered.¹³ 

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_αβ⁻¹. Clearly, this is a long-range force that obliges one to consider force contributions even from particles that are widely separated.¹⁶ 

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 surface-to-surface 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 non-hydrodynamic 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 molecular-like treatment of the problem where only particle–particle interactions (dependent on the relative positions) are considered.

Although molecular-like SD simulations require a smaller computational effort in comparison to mesh-based CFD ones, they are still computationally expensive as a dense matrix (remember that hydrodynamic interactions are long-ranged 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 long-range 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.²⁵  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 two-way 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 many-body problem of particle dynamics is usually solved through molecular-like simulations (such as molecular dynamics or discrete element method) while the flow problem is addressed using mesh-based 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 two-way 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,²⁸  non-spherical particles,^29–31  rotational degrees of freedom,^32,33  polydispersity^34–36  or complex non-hydrodynamic 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 one-way 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.

In Section 1.1 it has been shown that theoretically it is possible to compute particle trajectories in a CMRF as long as the non-hydrodynamic 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 ferri-magnetic 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⁻¹, magnetization reaches values of the order of 1000 kA m⁻¹ (iron), 300 kA m⁻¹ (nickel) and 100 kA m⁻¹ (cobalt).⁴⁰ 

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 non-zero. 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.⁴¹ 

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 so-called internal field H⃗_int) with some phenomenological features:⁴²  (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 non-linear response and assumes the absence of coercivity or anisotropy.

There are several models for the magnetization-field constitutive equation. In the literature, the Fröhlich–Kennelly model is commonly chosen:

where μ_r(H_int) is the scalar and field-dependent 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:⁴¹ 

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 non-linear 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 non-homogenous 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 self-consistent approach. Nevertheless, from first principles, the total field H⃗ at any point can be obtained from Maxwell equations:⁴⁴ 

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 ε₀ and μ₀ 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 quasi-static 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.⁴⁵ 

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 well-known expressions from Lorentz or Kelvin⁴⁶  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 second-order 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.⁴⁷  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:⁴⁸ 

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⃗ = μ₀μ_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 μ₀μ_rH²/2 = BH/2. For non-linear 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,⁴⁹  multipole expansion⁵⁰  or boundary elements³⁸  have been applied in the literature. Nevertheless, the aforementioned solutions are restricted to linear materials (that is, not applicable to ferro- or ferri-magnetic 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 non-linear, 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 mesh-based 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.

Alternatively, one can decide to solve eqn (1.34)–(1.37) up to a first order approximation using the so-called ‘Mean Magnetization Approximation’ (MMA). This is based on a well-known 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:⁵⁵ 

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_pM⃗ 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 non-magnetic particles dispersed in a linear magnetic media also interact, at least, through dipoles. This allows one to consider inverse ferrofluids, dispersions of non-magnetic 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 restoring-like 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 chain-like 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². On the contrary, at high magnetic fields when the material saturates and M² tends to M_S², 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 non-uniformity 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⁵⁶  (as in inverse ferrofluids and magnetic iron oxide-doped latex suspensions), the MMA can reproduce the complete solution. In other cases, the true magnetic interaction will be highly underestimated.

As in any suspension, particles in a CMRF experience other non-hydrodynamic forces besides magnetostatic ones. Most of them are classic interactions in colloidal science, so we refer to well-known 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 surface-to-surface 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⁻⁵ or smaller.⁵⁹  As a result, van der Waals forces can be safely neglected when accounting for other non-hydrodynamic 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 self-assembly 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 double-layer, 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.⁶⁰ 

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_BT/d where k_B is Boltzmann's constant (1.381 × 10⁻²³ J K⁻¹) 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 so-called λ 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 long-ranged hydrodynamics, in the macroscopic sample behavior.⁶¹ 

Different expressions have been proposed in the literature to model contact interactions, although current works in granular rheology usually choose Hertzian repulsion and Coulomb-like 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 non-magnetic 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.

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 macro-descriptions.

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 micro-configuration, 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:¹ 

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:¹ 

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 non-hydrodynamic 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 non-hydrodynamic 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 so-called 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 non-hydrodynamic 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 non-hydrodynamic 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 non-hydrodynamic 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 non-hydrodynamic forces, but also because there is a continuous phase able to mediate the non-hydrodynamic 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 long-range 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 non-hydrodynamic 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 long-range 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 non-hydrodynamic force experienced by the particle in the total stress tensor, the so-called term⁶⁴  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 non-hydrodynamic force acting over it. Although the example proposed here only considers long-ranged non-hydrodynamic forces, eqn (1.57) is valid for hard-sphere 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:

1. 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.

2. In the absence of flow and depending on the particle positions, due to the term, the system can develop shear stresses (non-diagonal 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.

3. 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 shear-thinning behavior.

4. 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.

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 field-induced 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.⁶⁵ 

Numerical studies (using magnetic dipolar forces, the Stokes’ drag approximation and a one-way 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⁷⁵  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.⁸⁰  The latter allows the structure to accommodate by diffusion during the short period when the external field is absent. This approach results in the long-term formation of disperse ‘drops’ consisting of particles that are arranged in a BCT lattice.^81,82 

Undoubtedly, the yield stress τ₀ is one of the most relevant figures of merit in an MRF.⁸³  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 τ₀, 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 τ₀, it is said that the sample leaves the elastic or pre-yield regime, breaks (also called sample failure) and starts to flow or to strain continuously.⁸⁴ 

However, the previous definition together with its related phenomenology is rarely seen experimentally (e.g. Wang et al.⁸⁵ ). 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 τ₀. 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.⁸⁶  As can be seen, the behavior at stresses above the aforementioned interval perfectly agree with the definition of τ₀. On the contrary, the appearance of the Newtonian behavior below τ₀ is controversial as the MRF is not perfectly elastic.

This discussion and related controversy around τ₀ measurement or existence are well-known and not new and highlight the paradox between what flows or not: to properly demonstrate that a material is elastic for any stress below τ₀, 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 τ₀ 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, τ₀ 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.⁸⁹ 

Yield stress ambiguity has given rise to several yield stress definitions. In this chapter, we will deal with two of them: the so-called static and dynamic yield stresses. Static yield stress will be discussed in this section and is identified with τ₀, that is, the lowest stress required for the onset of irreversible flow.⁸⁹  Leading on from this, the dynamic yield stress (τ_y) will be introduced in Section 1.4.3.

Motivated by torque-transfer applications in electro-mechanical 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 τ₀∝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 many-body problem both from an experimental and theoretical point of view.^107,108  Generally speaking, for CMRFs at low-moderate 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 pre-yield regime (i.e. absence of flow), therefore the strain suffered by the sample takes a very long time to relax, so that the only non-isotropic contribution to the total stress tensor is given by eqn (1.57). This depends on the non-hydrodynamic 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 non-hydrodynamic 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 τ₀ = f(γ₀) 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 τ₀ with our definition of static yield stress understood as the onset of flow. The related strain γ₀ will be called here the ‘yield strain’ despite the fact that some references save this name for the limit of the linear elastic behavior.¹⁰⁹  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 non-hydrodynamics 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 field-induced 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 one-way 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 pre-yield 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⁻³ 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.

In the macroscopic approaches the particles are supposed to arrange following the applied field direction in fusiform (cylinders or ellipsoids) or sheet-like 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¹¹¹ ):

where m_z = V(1 − μ_zz⁻¹)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.¹¹² ) 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 +γ²)² 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 gap-spanning 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.¹¹²  Similar to what happened with the first approximation using the one-way coupling scheme, the theoretical predictions are below the experimental data. Clearly, a good match was not expected at intermediate-high 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).¹⁰⁷  Besides, it must be said that these macroscopic models have been also derived for the case of non-linear materials (thus involving a numerical resolution) and grounding the computation of the force on Maxwell's stress tensor.¹¹⁴  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 pre-yield behavior based on the continuum theory at small strains. In this context, the CMRF is seen as a homogeneous body with a field-independent 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.¹¹⁵  Of course, the preferred direction is that along the field-induced 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.

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 (L-th order solution) have been implemented.^49,50 

In the case of non-linear 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 non-linear 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 non-linear 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, sub-quadratic (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 sub-quadratic 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.¹²⁸  Note that, in any regime, the yield stress is always proportional to the volume fraction. This comes from initially supposing isolated non-interacting chains. Assuming a homogeneous distribution of chains, it can be shown that the area per chain w² goes as w²∝ϕ⁻¹, thus τ₀ = F_m/w²∝ϕ.¹²⁹ 

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₃ is a proportionality constant that depends on the particular model assumptions: c₃ = 0.086 for an axisymmetric model,¹²⁷ c₃ = 0.114 for a pair or particles,¹³⁰ c₃ = 0.123 for a semi-infinite chain¹³¹  or c₃ = 0.137 for an infinite chain.¹³²  More complex models can be also found, for example, considering multipoles for a pair of linear particles,¹³³  or ‘bead-rod’ models where the resulting yield stress is dependent on the confinement of the sample.¹³¹ 

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 non-hydrodynamic 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² (see eqn (1.48)–(1.50)), it can be further written τ_p = f(ϕ,)M². 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² we arrive to:

The ratio M²/ is traditionally introduced through dimensional analysis. It is the so-called Mason number Mn and can be conceived as the ratio between hydrodynamic (drag) and magnetostatic force scales, advection and magnetic-controlled 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.¹⁷  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 re-configuration 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 .¹¹⁰  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.⁸⁹  From here and remembering that τ_p(ϕ,) is proportional to M², it is straightforward to see that the viscosity will diverge as Mn decreases in the magnetostatic regime: η = τ/ ∼ τ_y(ϕ)/ = f(ϕ)M²/∝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 non-interacting 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.¹³⁴  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¹⁷ or two-way coupling schemes.⁸  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 field-induced 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/(η_∞).⁸ 

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 x-axis 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 non-hydrodynamic 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).¹⁷  To enhance the fitting in this regime, other macroscopic models such as the Casson one or microstructural viscosity models have been used.¹⁴³ 

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 non-Brownian hard sphere suspensions. The solution was provided by Einstein¹  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 non-negligible hydrodynamic interactions that require simulation techniques to be solved. As it has already been pointed out, these simulations are highly time/computational-power 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.⁸⁹  In this chapter, the so-called Quemada expression is used:¹⁴⁴ 

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₄ϕ where c₄ is a constant that depends on the particular forces and the specific mechanical stability condition: c₄ = 5.25,¹³¹  c₄ = 8.49,¹³⁷  or c₄ = 8.82.¹⁴⁵  At large concentrations, work has been done with one-way coupling schemes resulting also in a linear dependence with the volume fraction, in this case, c₄ = 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 one-way 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 non-hydrodynamic 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,¹⁷  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 one-way schemes as non-hydrodynamic forces are properly taken into account as well. Main claims could arise from 〈S̃_p〉 that requires the use of two-way 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).

The appearance of non-isotropic normal stresses, as well as normal stress differences, in simple shear flow together with the shear rate-dependent viscosity and the yield stress are classical signatures of non-Newtonian flow behavior. From eqn (1.58) and taking into account the field-induced 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 sample-geometry 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². However, consensus does not seem to have been reached regarding the strain/shear rate dependence of the normal force in the pre-/post-yield regime. In the pre-yield 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.¹⁴⁸ 

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 field-direction component of the force instead of the shear-direction one. Thus, results share the same properties, σ_zz∝μ₀M²ϕ, 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 non-strained 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¹⁴⁹  and Maxwell's stress jump at the boundaries¹⁴⁸ ), 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.¹⁴⁸  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.¹⁷  report (σ_yy − σ_zz)/ < 0 at any Mn. This quantity decreases, in absolute value, with Mn until it vanishes in the hydrodynamic regime.

Magnetorheological fluids are systems whose micro-dynamics and macroscopic behavior are dominated mainly by the hydrodynamic and magnetostatic interactions between dispersed particles. Both kinds of forces are multibody and non-linear 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 macro-scales 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 one-way 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.

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

Non-hydrodynamic 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

Non-hydrodynamic 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

δ̃

Second-order identity tensor

ε₀

Vacuum permittivity

η_c

Viscosity of the continuous phase

μ₀

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

τ₀

Static yield stress

τ_y

Dynamic yield stress

γ

Shear strain

t

Time

G

Shear modulus

γ₀

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²

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

This work was supported by ERDF, FEDER, MICINN AEI PID2019-104883GB-I00, Junta de Andalucía P18-FR-2465 and A-FQM-396-UGR20 projects. J. R. Morillas acknowledges FPU14/01576 fellowship and European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant (EFST)-H2020-MSCA-IF-2020 (Grant 101030666).