On the Apparent Yield Stress in Non-Brownian Magnetorheological Fluids

TL;DR

This study uses simulations to explore how viscous force laws influence the apparent yield stress and flow behavior in dense magnetorheological fluids, revealing conditions under which a true yield stress emerges.

Contribution

It demonstrates that the form of viscous dissipation critically affects the flow exponent and yield stress presence in magnetorheological fluids, linking microscopic clustering to macroscopic rheology.

Findings

Models with contact damping exhibit a clear yield stress plateau.

Stokes-like drag models show no apparent yield stress.

Increasing volume fraction induces yield stress near jamming.

Abstract

We use simulations to probe the flow properties of dense two-dimensional magnetorheological fluids. Prior results from both experiments and simulations report that the shear stress $\sigma$ scales with strain rate $\dot \gamma$ as $\sigma \sim \dot \gamma^{1-\Delta}$, with values of the exponent ranging between $2/3 <\Delta \le 1$. However it remains unclear what properties of the system select the value of $\Delta$, and in particular under what conditions the system displays a yield stress ($\Delta = 1$). To address these questions, we perform simulations of a minimalistic model system in which particles interact via long ranged magnetic dipole forces, finite ranged elastic repulsion, and viscous damping. We find a surprising dependence of the apparent exponent $\Delta$ on the form of the viscous force law. For experimentally relevant values of the volume fraction $\phi$ and the dimensionless Mason number (which quantifies the competition between viscous and magnetic stresses), models using a Stokes-like drag force show $\Delta \approx 0.75$ and no apparent yield stress. When dissipation occurs at the contact, however, a clear yield stress plateau is evident in the steady state flow curves. In either case, increasing $\phi$ towards the jamming transition suffices to induce a yield stress. We relate these qualitatively distinct flow curves to clustering mechanisms at the particle scale. For Stokes-like drag, the system builds up anisotropic, chain-like clusters as the Mason number tends to zero (vanishing strain rate and/or high field strength). For contact damping, by contrast, there is a second clustering mechanism due to inelastic collisions.

Tables

Table 1. Table 1: Previous work determining the exponent Δ Δ \Delta . These works have been performed under a wide variety of conditions, including variations in the type of particles, carrier fluid, and system geometry. For details about the specific parameters used in each experiment/simulation, readers are referred to the original articles.

Authors	Type	$\Delta$
Marshall et al. Marshall et al. (1989)	experiment ER-fluid	1
Halsey et al. Halsey et al. (1992)	experiment ER-fluid	$0.68$ – $0.93$
Felt et al. Felt et al. (1996)	experiment MR-fluid	$0.74$ – $0.83$
Martin et al. Martin et al. (1994)	experiment ER-fluid	$0.67$
de Gans et al. de Gans et al. (1999)	experiment inverse MR-fluid	$0.8$ – $0.9$
de Gans et al. De Gans et al. (2000)	experiment inverse MR-fluid	$0.94\pm 0.02$
Volkova et al. Volkova et al. (2000)	experiment a) MR-fluid b) inverse MR-fluid	a) $0.86$ – $0.97$ b) $0.74$ – $0.87$
Sherman et al. Sherman et al. (2015)	experiment MR-fluid	1
Bonnecaze and Brady Bonnecaze and Brady (1992)	simulation 2D	1
Melrose Melrose (1992)	simulation 3D	$0.8\pm 0.05$

Equations

$$\frac{\eta }{{\eta }_0}=1+{\left(\frac{{\text{Mn}}^{\ast }}{\text{Mn}}\right)}^{\Delta }$$

$$\sigma (\dot{\gamma })={\sigma }_y+A\dot{\gamma }.$$

$${\ddot{\mathbf{r}}}_i=\frac{1}{m_i}\left({\mathbf{f}}_i^{\text{e}}+{\mathbf{f}}_i^{\text{d}}+{\mathbf{f}}_i^{\text{m}}\right),$$

$$\varphi =\frac{1}{L^2}\sum _{i=1}^N\frac{\pi d_i^2}{4},$$

$$U_{e}(r_{ij})=\left\{\begin{array}[]{cl}\frac{k_{e}}{\alpha}\left(1-r_{ij}/d_{ij}\right)^{\alpha}&\mbox{ for $r_{ij}<d_{ij}$}\\ 0&\mbox{ for $r_{ij}\geq d_{ij}$}\end{array}\right.$$

$${\mathbf{f}}_i^{\text{e}}=-\sum _{j(i)}\frac{d{U}_e(r_{ij})}{d{r}_{ij}}{\hat{\mathbf{r}}}_{ij}=-k_e\sum _{j(i)}{\left(\frac{d_{ij}-r_{ij}}{d_{ij}^2}\right)}^{\alpha -1}{\hat{\mathbf{r}}}_{ij},$$

$${\mathbf{f}}_i^{\text{e}}=-k_e\sum _{j(i)}\frac{d_{ij}-r_{ij}}{d_{ij}^2}{\hat{\mathbf{r}}}_{ij}.$$

$${\mathbf{f}}_i^{\text{RD}}=-k_d\left({\mathbf{v}}_i-{\mathbf{v}}_{\text{RD}}\right),$$

$${\mathbf{f}}_{ij}^{\text{CD}}=-k_d\left({\mathbf{v}}_i-{\mathbf{v}}_j\right).$$

$$\mathbf{B}(\mathbf{r})=\frac{{\mu }_f}{4\pi }\left(\frac{3\mathbf{r}(\mathbf{m}\cdot \mathbf{r})}{r^5}-\frac{\mathbf{m}}{r^3}\right),$$

$$U_m({\mathbf{r}}_{ij})=-\left({\mathbf{m}}_j\cdot {\mathbf{B}}_i\right),$$

$${\mathbf{f}}_{ij}^m=-\nabla U_m(r_{ij}).$$

$${\mathbf{f}}_{ij}^m$$

$$\displaystyle\,\,\,\,\,\,\,+\mathbf{r}_{ij}(\mathbf{m}_{i}\cdot\mathbf{m}_{j})-\frac{5}{r_{ij}^{2}}\mathbf{r}_{ij}(\mathbf{m}_{i}\cdot\mathbf{r}_{ij})(\mathbf{m}_{j}\cdot\mathbf{r}_{ij})\Big{]}\,.$$

$${\mathbf{m}}_i=V_{\mathrm{ci}}\mathbf{M}=V_{\mathrm{ci}}(3\beta \mathbf{H}),$$

$$\beta =\frac{\mu -1}{\mu +2}.$$

$$\sigma ={\sigma }_e+{\sigma }_m+{\sigma }_d+{\sigma }_k.$$

$${\sigma }_{\circ }=-\frac{1}{L^2}\sum _{i<j}{r}_{ijx}{f_{\circ }}_{ijy}$$

$${\sigma }_k=-\frac{1}{L^2}\sum _i{v}_{ix}{v}_{iy}{m}_i$$

$${\text{Mn}}_{\text{RD}}=\frac{3\pi k_d\dot{\gamma }}{8d{\mu }_f{\beta }^2{H}^2}$$

$${\text{Mn}}_{\text{CD}}=Z{\text{Mn}}_{\text{RD}}$$

$$\tilde{\sigma }=\frac{\sigma d^{D-3}}{{\mu }_f\beta H^2}.$$

$${\theta }_H=\frac{1}{NZ}\sum _{ij}\max \left(\theta ({\mathbf{r}}_{ij},\mathbf{H}),\theta (-{\mathbf{r}}_{ij}\mathbf{H})\right).$$

$$\mathbf{B}={\mu }_f\left(\mathbf{H}+\mathbf{M}\right),$$

$${\mathbf{m}}_i=V_{\mathrm{ci}}\mathbf{M}=V_{\mathrm{ci}}(3\beta \mathbf{H}),$$

$$\beta =\frac{\mu -1}{\mu +2}.$$

$$\mu =\frac{{\mu }_i}{{\mu }_f},$$

$${\mathbf{m}}_i=3V_{ci}\beta \left[\mathbf{H}+\frac{1}{{\mu }_f}\sum _{i\ne j}{\mathbf{B}}_{ij}\right],$$

$${\mathcal{O}}_i=\sum _{r_{ij}<r_c}{\mathcal{O}}_{ij}+{\mathcal{O}}_{LR}.$$

$${\mathbf{B}}_j(\mathbf{r})=\frac{{\mu }_0}{4\pi }\left(\frac{3\mathbf{r}({\mathbf{m}}_j\cdot \mathbf{r})}{r^5}-\frac{{\mathbf{m}}_j}{r^3}\right).$$

$${\tilde{\mathbf{m}}}_i=\frac{1}{\pi r_c^2}\sum _{|r_{ij}|<r_c}w(r_{ij}){\mathbf{m}}_j,$$

Full text

Magnetorheological (MR) fluids consist of magnetizable particles suspended in a viscous carrier fluid. An external magnetic field $\mathbf{H}$ induces magnetic moments in the particles, which then rearrange to form chain-like structures, as illustrated in Fig. 1. Chain formation dramatically enhances the stress $\sigma$ needed to maintain a strain rate $\dot{\gamma}$, and by varying $\mathbf{H}$ it is possible to tune the viscosity of the suspension, with applications to damping and switching. An excellent introduction to the fundamental physics and engineering applications of MR fluids can be found in recent review articles by Vicente et al. de Vicente et al. (2011) and Ghaffari et al. Ghaffari et al. (2015) and references therein. Here we numerically study non-Brownian MR fluids in steady shear flow.

Steady state rheology is commonly characterized in terms of the enhancement of the shear viscosity $\eta={\sigma}/{\dot{\gamma}}$ over its value $\eta_{0}$ at zero field. The ratio $\eta/\eta_{0}$ is governed by the dimensionless Mason number $\text{Mn}\propto\dot{\gamma}/H^{2}$ (discussed in detail below), which quantifies the relative strengths of viscous and magnetic stresses in the system; magnetic interactions dominate when Mn tends to zero. In practice, the empirical fitting function

[TABLE]

is often found to give a good description of the viscosity enhancement in MR fluids. Here $\text{Mn}^{\ast}$ is a function of the volume fraction $\phi$; it vanishes as $\phi\rightarrow 0$ and determines the crossover between the Newtonian flow regime $\eta/\eta_{0}\sim 1$ at high Mn and the magnetically dominated regime $\eta/\eta_{0}\sim\text{Mn}^{-\Delta}$ at low Mn. The exponent $\Delta$ controls the rate at which viscosity diverges as the Mason number decreases.

The value $\Delta=1$ is an important reference case, as Eq. (1) is then equivalent to the flow curve of a Bingham plastic

[TABLE]

The Bingham plastic has a nonzero dynamic yield stress $\sigma_{y}$, defined here as the asymptote of the steady state flow curve $\sigma(\dot{\gamma})$ in the limit $\dot{\gamma}\rightarrow 0$ (henceforth “the yield stress”). Experiments are of course performed at small but finite strain rates, hence in practice the yield stress is also identified with an apparent plateau in the flow curve at the lowest accessible rates; i.e. one assumes the plateau persists to asymptotically low strain rates. If instead $\Delta<1$, then the system has no yield stress and the stress vanishes slowly with the strain rate, $\sigma\sim\dot{\gamma}^{1-\Delta}$ as $\dot{\gamma}\rightarrow 0$. Competing theoretical descriptions predict exponents $\Delta=1$ Marshall et al. (1989); Klingenberg and Zukoski IV (1990); Martin and Anderson (1996); de Vicente et al. (2004); de Gans et al. (1999) and $\Delta=2/3$ Halsey et al. (1992); for a discussion of the different models see Bossis et al. (2002); Ruiz-López et al. (2016); Ghaffari et al. (2015). Numerous experimental and numerical studies have measured $\Delta$ values throughout this range in a number of magnetorheological (and electrorheological) systems under varying conditions; a summary of their results is given in table 1. Thermal motion has been shown to give $\Delta<1$ both in experiment de Gans et al. (1999) and simulation Baxter-Drayton and Brady (1996). However there is no effective way to predict a priori whether a given non-Brownian MR fluid will display a yield stress. In the present work we address the presence or absence of a yield stress in a two-dimensional numerical model of athermal MR fluids in which energy is dissipated via viscous forces, in the absence of Coulomb friction.

We model composite particles with a core shell structure as seen in Fig 1a. Similar core-shell structures have been used in experiment to change the surface properties of particles and to lower the effective density of the particles in order to avoid sedimentation You et al. (2007); Choi et al. (2007); Ko et al. (2009); Fang et al. (2010, 2009). The core shell structure is also suitable to model the recent experiments of Cox et al. (2016), which bridge the gap between conventional magnetic suspensions and amorphous magnetic solids.

We consider two types of damping. The first, which we denote reservoir damping (RD) in accord with the terminology of Ref. Vågberg et al. (2014), is a Stokesian drag with respect to the carrier fluid. The second, contact damping (CD), is applied to the relative velocity of particles in contact. Removing friction gives us a cleaner system that helps us to understand the underlying physics. We will ultimately argue that contact damping provides insight into the case of frictional contacts. The use of RD and CD models also allows us to make contact with the extensive literature Durian (1995); Liu and Nagel (1998); O’Hern et al. (2003); Olsson and Teitel (2007); Heussinger and Barrat (2009); Hatano (2009); Tighe et al. (2010); Tighe (2011); Ikeda et al. (2012); Boschan et al. (2016); Baumgarten et al. (2017) on the rheology of yield stress fluids close to the jamming transition.

Our main finding is that the form of the viscous force law has a dramatic influence on the viscosity enhancement in the magnetically dominated regime, as characterized by the exponent $\Delta$. For reservoir damping we find no evidence of a yield stress over a wide range in $\phi$ and for Mason numbers as low as $10^{-6}$; instead the exponent $\Delta\approx 0.75$ gives an excellent description of the rheology. In sharp contrast, for contact damping there is a clear nonzero yield stress in the same range of Mason numbers, and so $\Delta=1$. We relate this difference to clusters that form in the CD model at intermediate Mn due to inelastic collisions between particles, an effect that is absent in the RD model.

We further investigate the role of finite size effects and volume fraction, both of which we find to promote the emergence of a yield stress.

I Model

The system comprises $N$ spherical particles confined to a plane. The particle distribution is a 50-50 bidisperse mixture with size ratio $1.4:1$. Each particle consists of an elastic non-magnetic outer layer and a hard inner core of a magnetically soft permeable material (Fig. 1a). The diameter of the core is $d_{i}/2$, where $d_{i}$ is the diameter of particle $i$. The mass of each particle $m_{i}$ is directly proportional to its volume $V_{i}$ such that $m_{i}=V_{i}\rho$, where $\rho=1$ is the density of the material. For simplicity we assume the density is constant throughout the particle. We assume the particles are large enough so that thermal motion can be neglected and that there is no static friction.

The motion of each particle is governed by

[TABLE]

where $\mathbf{r}_{i}$ is the position of particle $i$, $\mathbf{f}^{\text{\,e}}_{i}$ is the repulsive contact force, $\mathbf{f}^{\text{\,d}}_{i}$ is a dissipative force caused by the interaction between the particles and the surrounding liquid, and $\mathbf{f}^{\text{\,m}}_{i}$ is the magnetic dipole force. Since the particles are frictionless and do not have permanent dipole moments, there are no torques acting on the particles.

We drive the system by applying a uniform shear strain rate $\dot{\gamma}$ in the $\hat{\mathbf{x}}$ direction using Lees-Edwards boundary conditionsEvans and Morris (1990). The equations of motions were integrated using a Velocity-Verlet scheme modified to better handle dissipative forces Vågberg et al. (2017).

The simulation is controlled by varying three parameters: the shear rate $\dot{\gamma}$, the external magnetic field $\mathbf{H}=H\hat{\mathbf{y}}$ transverse to the flow direction, and the packing fraction

[TABLE]

where $L$ is the side length of the quadratic simulation box.

I.1 Interaction Forces

Overlapping particles repel elastically (Fig. 2a). The elastic contact forces are given by the potential

[TABLE]

where $r_{ij}=\left|\mathbf{r}_{ij}\right|=\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|$ is the distance between particle $i$ and $j$ and $d_{ij}$ is the sum of their radii. The constant $k_{e}=1$ sets the energy scale of the elastic interaction. For the parameter ranges studied here the particle overlaps are small, $d_{ij}-r_{ij}\ll d_{ij},$ so the contact interaction is limited to the outer shell; it is therefore safe to neglect the particle core in the contact potential.

The potential (5) produces an elastic force

[TABLE]

where the sums run over the set of particles $j$ in contact with particle $i$. Using $\alpha=2$ gives the standard harmonic potential with corresponding force

[TABLE]

For the dissipative force $\mathbf{f}^{\text{d}}_{i}$ we use a viscous force proportional to the velocity difference between the particle velocity $\mathbf{v}_{i}$ and a reference velocity. We compare two different viscous dissipations (Fig. 2b and c), by changing the definition of the reference velocity.

With the first viscous force law, which we denote reservoir dissipation (RD), the particle loses energy when moving relative to the carrier fluid. We select the reference as $\mathbf{v}_{\text{RD}}=\dot{\gamma}y_{i}\hat{\mathbf{x}}$, the affine shear velocity. This gives

[TABLE]

where the constant $k_{d}$ allows us to tune the strength of the dissipation.

The second force law is a contact dissipation (CD), wherein the dissipation is proportional to the velocity difference of contacting particles

[TABLE]

To obtain the full dissipative force on particle $i$ one must sum over all particles $j$ in contact with $i$.

We use $k_{d}=1$ for both the RD and CD dissipation. For RD this ensures the dynamics is overdamped for the studied parameter ranges. While the CD-dissipation is overdamped for contacting particles, it is highly sensitive to the average contact number and free particles do not dissipate energy. This mainly affects the behavior of dilute systems at high shear rates, which is not the limit we focus on here.

The RD and CD force laws can be seen as two limiting cases: RD only considers the particle-carrier fluid interaction, while CD only considers the particle-particle interaction. The two force laws have been studied in detail for dense suspensions in the absence of dipole-interactions Vågberg et al. (2014, 2014, 2017). In experimental systems both solvent and particle interactions affect the dissipation and a combination of CD- and RD-dissipation are usually needed to describe the behavior. Simulations are advantageous, in that they allow us to study these effects separately.

The magnetic interaction is modeled using point dipoles positioned at the center of each particle – see Fig. 2d. The dipole moments are induced in the particle core by the external field $\mathbf{H}$. The magnetic flux density $\mathbf{B}$ at a distance $r$ from a dipole $\mathbf{m}$ is given by

[TABLE]

where $\mu_{f}$ is the permeability of the carrier fluid. The potential energy between two dipoles $i$ and $j$ is given by

[TABLE]

which gives the force

[TABLE]

Inserting (10) and (11) into (12) and evaluating gives the force from dipole $i$ acting on dipole $j$,

[TABLE]

The magnitude and direction of the induced dipole-moments are given by

[TABLE]

where $V_{\rm ci}$ is the core-volume of particle $i$, and

[TABLE]

Here $\mu=\mu_{i}/\mu_{f}=1000$ is the relative permeability and $\mu_{i}$ is the permeability of the core of particle $i$. The outer shell is assumed to have the same permeability as the carrier fluid. We consider only direct induction from the external field, ignoring contributions form neighboring dipoles. This is justified by the core-shell structure of the particles, which keeps the magnetic cores separated. We refer to the appendix for a more detailed discussion.

I.2 Stresses

The shear stress $\sigma$ is a sum of four contributing terms

[TABLE]

Each of the first three correspond to one of the forces in (3). The additional term $\sigma_{k}$ is the kinetic stress. In practice only $\sigma_{e}$ and $\sigma_{m}$ are important for the rheology in the magnetically-dominated regime, as $\sigma_{d}$ and $\sigma_{k}$ are orders of magnitude lower and both go to zero in the quasistatic limit $\dot{\gamma}\rightarrow 0$.

The first three stress terms are calculated by substituting $f_{\circ}$ with the corresponding force from equation (3) in the expression

[TABLE]

Here the $x$ and $y$ subscripts indicate the $x$- and $y$-components of respective vector and $L$ is the length of the simulation box. The kinetic stress $\sigma_{k}$ is calculated as

[TABLE]

where $v_{ix}$ and $v_{iy}$ is the x- and y-component of the non-affine velocity of particle $i$.

I.3 Dimensionless numbers

Much of the observed rheology of MR-fluids can be described using four dimensionless numbers: the Mason number (Mn), the Peclet number (Pe), Lambda ($\lambda$), and the volume fraction $\phi$. The first three characterize the relative strengths of magnetic, viscous, and thermal forces. As we consider non-Brownian particles, the Peclet number (viscous versus thermal forces) and Lambda (magnetic versus thermal forces) play no role in the present results. We are left with the volume fraction and the Mason number, which vanishes when magnetic interactions dominate.

There is some flexibility when selecting the reference forces used to define the Mason number, which has led to competing conventions in the literature Klingenberg et al. (2007). We use microscopic properties to define Mn. Assume there are two particles of the smaller species with diameter $d$ (core diameter $d/2$) placed at a distance $d$ such that their surfaces just touch. The dipole force between these two particles when their dipole moments are aligned is $F_{m}=\frac{3\pi}{8}d^{2}\mu_{f}\beta^{2}H^{2}$. For reservoir damping the typical viscous force is $F_{d}=dk_{d}\dot{\gamma}$, while for contact damping there is an additional dependence on the mean number of contacts per particle, $Z$, such that $F_{d}=Zdk_{d}\dot{\gamma}$. The Mason number $\text{Mn}\equiv F_{d}/F_{m}$ is therefore

[TABLE]

for the RD model and

[TABLE]

for the CD model.

We report shear stresses in the dimensionless form

[TABLE]

where $D$ is the dimensionality of the system. Because the presence or absence of a yield stress is a major focus of the present work, we present most rheological results in the form of a dimensionless flow curve, $\tilde{\sigma}(\text{Mn};\phi)$, as opposed to plotting the viscosity enhancement $\eta/\eta_{0}$. A yield stress is then clearly signaled by a plateau in $\tilde{\sigma}$ at low Mn. When there is no yield stress, the stress vanishes as $\tilde{\sigma}\sim\text{Mn}^{1-\Delta}$.

I.4 Simulation

The length of each simulation in total strain $\gamma$ varies from $\gamma=50$ at $\dot{\gamma}=0.05$ down to $\gamma=4$ at $\dot{\gamma}=10^{-8}$. Simulations are started at high shear rates, and lower shear rate simulations are initialized using starting configurations obtained from the previous higher shear rate. In order to avoid transient effects the first 20% of each run is discarded before calculating time-averaged quantities. For $N\geq 1024$ we perform one simulation for each parameter value, while for $N=256$ and $N=64$ two, respectively five, independent runs are performed to improve statistics. We study the parameter range $0.1\leq\phi\leq 0.86$, $10^{-8}\leq\dot{\gamma}\leq 10^{-1}$, $10^{-4}\leq H\leq 10^{-1}$ and $64\leq N\leq 16384$, which allows us to probe Mason numbers in a window spanning 12 orders of magnitude for $N=256$ and 10 orders of magnitude for $N=4096$. Consequently, we cover a larger parameter space than any of the works referenced in table 1.

For this work we are especially interested in the behavior at low Mason numbers. At $N=4096$ our lowest Mason number is $\text{Mn}=5\times 10^{-7}$, which is significantly lower than the the lowest values accessed by any of the simulations in table 1 and comparable to or slightly lower than the lowest values accessed in experiment Marshall et al. (1989); Felt et al. (1996); Martin et al. (1994); de Gans et al. (1999).

II Bulk rheology

We start by considering the bulk rheological properties of the RD and CD models, with emphasis on the form of their steady state flow curves.

Fig. 3 compares the rheology of the RD and CD models and its dependence on $\dot{\gamma}$ and $\mathbf{H}$ at fixed $\phi=0.5$ and $N=4096$. We first consider rheology of the RD model, shown in the left column of Fig. 3. From top to bottom we plot the same data set as dimensionful flow curves, dimensionless flow curves, and in terms of the viscosity enhancement, respectively. The dimensionless data displays excellent collapse to a master curve that exhibits two flow regimes: a Newtonian regime, $\sigma\sim\text{Mn}$, at high Mason numbers, and a magnetically dominated regime at low Mason numbers. It is clear that the RD model does not exhibit a yield stress over the accessible range of Mn; instead we find $\Delta\approx 0.75$ in the magnetically-dominated regime. The corresponding panels for the CD model (Fig. 3, right column) display a striking difference. There are again two flow regimes, but in this case there is a more gradual crossover to a yield stress in the limit of low Mn, hence $\Delta=1$.

It is natural to ask if the qualitative differences in the flow curves of Fig. 3 are due to finite size effects. To answer this question, we simulate steady state shear flow for a range of system sizes $N=64\ldots 16384$; the corresponding flow curves for the RD and CD models are plotted in Fig. 4a and b, respectively. In both cases, we obtain good data collapse over the entire sampled range of Mn, independent of $N$. We therefore conclude that differences between the RD and CD flow curves are not due to finite size effects.

Boundary effects are closely related to finite size effects. They are also particularly relevant to experiments, of course, as shearing surfaces are necessary to sustain flow. To probe the influence of boundaries on the flow curve, we introduce a wall by fixing the positions of a thin layer of particles intersecting the line $y=0$ (in the center of the cell). The resulting RD flow curves are plotted in Fig. 5. One clearly sees that the system with a wall develops a plateau at low Mn that is absent in the wall-free case for the same system size. This effect is clearly not a material property, but should be borne in mind when interpreting experimental data.

In order to quantify stress fluctuations in flow, we have also sampled the cumulative distribution function (CDF) of shear stress in steady state. In Fig. 6a and b we plot the CDF for a range of strain rates, with the highest chosen to correspond roughly to the “elbow” in the RD flow curve. While the curves shift left with decreasing $\dot{\gamma}$ (as already apparent from the flow curve), their overall shape changes little, indicating that stress fluctuations are insensitive to the strain rate. In Fig. 6c and d we plot the CDF for low Mn and a range of system sizes $N$. There is a slow systematic increase of the median stress (${\rm CDF}=0.5$) with $N$, which is too weak to be seen on the log-log plots of Fig. 4a and b. For small system sizes the flow regularly samples states with negative shear stress; however increasing the system size causes the CDF’s to sharpen, reducing the fraction of negative stress states. For $N=4096$ and $\text{Mn}=10^{-6}$ (Figs. 6c and d), the fraction sampled by the RD model is negligible, while in the CD model it is less than 0.1.

Based on the above results, we conclude that the bulk flow curve of the RD model at $\phi=0.5$ has no apparent yield stress over the experimentally (and numerically) accessible range of Mason numbers. The CD model, by contrast, does have an apparent yield stress.

II.1 Towards Jamming

We now consider the role of packing fraction $\phi$ in the bulk rheology. Intuitively, one expects the stress required to sustain steady flow to increase with $\phi$. Moreover, soft sphere packings in the absence of an applied field (i.e. $H=0$) are known to develop a yield stress at a critical volume fraction $\phi_{c}$ (the jamming point)Liu and Nagel (1998); O’Hern et al. (2003); Olsson and Teitel (2007). The precise value of $\phi_{c}$ depends and particle size distribution Koeze et al. (2016) as well as the protocol used to generate the packings Vågberg et al. (2011). For sheared systems in the quasistic limit and $H=0$ both the CD and RD model have been shown to jam at the same packing fraction $\phi_{c}\approx 0.8433$ Olsson and Teitel (2012); Vågberg et al. (2014).

It is therefore reasonable to ask what happens when the volume fraction is increased towards $\phi_{c}$ in the presence of a magnetic field $H>0$. We start by looking at the RD model.

In the top row of Fig. 7, panels (a-e), we plot RD flow curves for $\phi=0.4$, 0.6, 0.7, 0.82, and 0.84 at varying strain rate and field strength. For $\phi\leq 0.7$ we do not observe a plateau, although fitting a power law $\tilde{\sigma}\sim\text{Mn}^{1-\Delta}$ to the low Mn data reveals an effective exponent $\Delta$ approaching 1 as $\phi$ increases. For $\phi=0.82$ there is an unambiguous plateau at low Mn. Data above the plateau no longer collapse with Mn, which is an indication that critical effects near jamming have begun to play a significant role; at the same time, flow curves at $\phi=0.82$ and $H=0$ do not show a yield stress Vågberg et al. (2014). For $\phi=0.84$ the dynamics is completely dominated by the proximity to the jamming transition and data collapse with Mn is wholly absent. There are also strong finite size effects (as expected near a critical point), as seen in Fig. 4c. The flow curves at high Mn are no longer Newtonian but shear thinning – also a signature of the approach to jamming. For comparison, in the bottom row of Fig. 7, panels (f-j), we plot flow curves for the CD model for the same volume fractions; in all cases there is a plateau at low Mn, and we observe identical trends regarding data collapse with Mn.

In order to compare stresses at low Mn directly, we plot the stress over a range of volume fractions for constant $\text{Mn}=2\times 10^{-6}$ in both drag models – see Fig. 8. The stresses display an approximately exponential growth with $\phi$ over a wide range of volume fractions, before increasing more rapidly close to jamming. To test whether the flow curve has approached a plateau, we numerically evaluate the logarithmic derivative $q\equiv{\rm d}\,\ln{\tilde{\sigma}}/{\rm d}\,\ln{\dot{\gamma}}$ and plot the stress only when $q<0.2$ (filled symbols). For comparison we also plot the unfiltered stress (open symbols). It is apparent that the CD model always reaches a plateau (apart from a small number of outliers), while the RD model only shows a clear plateau at sufficiently high volume fractions. The particular value of $\phi$ where the plateau appears has some dependence on system size (compare panels (a) and (b)). While the stress in the RD model always exceeds that in the CD model, the two curves grow closer with decreasing Mason number. This is suggestive of convergence to a common asymptote, and therefore indirect evidence that the RD flow curves display a plateau at asymptotically low strain rates.

The data of Fig. 7 demonstrates that a plateau in the flow curve (i.e. an apparent yield stress) emerges in the RD model at sufficiently high volume fractions. We speculate that the plateau is present for all $\phi$ where the particles form a percolating cluster, which at lower $\phi$ values occurs for smaller Mason numbers than those accessed here. This hypothesis cannot be tested directly using present methods, but in the following Section we provide supporting evidence based on the evolution of microstructural measures with Mn and $\phi$.

III Microscopic structure

What microscopic features of the system correlate with (changes in) the bulk rheology? To gain insight into the qualitative differences between the RD and CD models apparent in the flow curves, we now seek to characterize the microstructural evolution of MR fluids in steady shear flow as a function of the Mason number and volume fraction. Our goals are twofold. First, at sufficiently low Mason numbers one expects MR fluids to quasistatically sample states that minimize the sum of the elastic and magnetic potential energies, with viscous forces playing a negligible role. Hence we will seek evidence that our simulations are approaching, if not definitively reaching, this asymptotic regime. Second, the qualitatively different flow curves in the RD and CD models should be reflected in their microstructure. Therefore we seek evidence of qualitative differences, in general, and competing clustering mechanisms in particular.

Snapshots of the system are presented in Fig. 9 for both the RD (top row) and CD (bottom row) model and several values of Mn (increasing from left to right). In the RD model there is an apparently homogeneous and isotropic microstructure in the Newtonian regime at high Mn. Chains gradually emerge as the Mason number is lowered and magnetic interactions increasingly dominate. The microstructural evolution in the CD model is comparatively complex. There is anisotropy even in the Newtonian regime. More strikingly, large clusters appear at intermediate Mn. These clusters are more compact than the chains that eventually form at low Mn, and which resemble those seen in the RD model. In the remainder of this Section we quantify the above observations.

III.1 Coordination

At asymptotically low Mason numbers, particles must follow quasistatic trajectories that track minima of the (magnetic and elastic) potential energy as parameterized by the strain coordinate; viscous dissipation can only play a subdominant role. Therefore the $\dot{\gamma}\rightarrow 0$ (and hence $\text{Mn}\rightarrow 0$ at fixed $H$ and $\phi$) limit of the flow curve $\sigma(\dot{\gamma};H,\phi)$, i.e. the “true” yield stress, must be the same in both the RD and CD models. To obtain evidence of the approach to this asymptotic limit, we now study the evolution of the mean contact number $Z$ at low Mn. $Z$ plays an important role in determining whether a network (e.g. the contact network of a soft sphere packing) can elastically support a load. Here we present evidence that microstructure is indeed independent of the damping mechanism in the limit of vanishing strain rate.

In the absence of a magnetic field, a packing jams (develops a shear modulus and yield stress) when it satisfies Maxwell’s Maxwell (1864) counting argument $Z\geq Z_{\rm iso}=2D+{\cal O}(1/L^{D-1})$, where $Z$ is the mean number of contacts per particle calculated after removing non-load bearing “rattlers” and $D$ is the spatial dimension. The correction term accounts for boundary effects. For several reasons, one expects magnetic interactions to generate elastically rigid states with mean contact numbers $Z<Z_{\rm iso}$. First, magnetic interactions enhance boundary effects due to clusters’ anisotropic shape Tighe and Vlugt (2011); Goodrich et al. (2012); Dagois-Bohy et al. (2012); Goodrich et al. (2014). They also introduce long range, potentially tensile forces between particles. The connectivity of the contact network still provides a useful characterization of the flow, however, because the tail of the magnetic interaction potential falls off rapidly with distance, so that the strongest magnetic forces are between nearest neighbors. Finally, when chains are present at low Mason numbers, to minimize the potential energy the particles will arrange such that nearest neighbor magnetic forces are nearly always tensile. Tensile forces increase the likelihood of a structure containing states of self stress, which reduce the number of contacts needed to render a structure rigid. Maxwell’s original counting argument can be extended to correctly count states of self stress as described by Calladine Calladine (1978), a procedure which has also been adopted for studying dense sphere packings Lois et al. (2008); Wyart et al. (2005); Head (2007).

We now empirically determine the scaling of $Z(\text{Mn})$ at low Mason number, including its asymptote $Z_{0}$ as Mn tends to zero. The contact number is a “bare” $Z$ with no correction for rattlers. Recalling that $Z_{\rm iso}\approx 4$ in large systems with no magnetic interactions, in Fig. 10a (crosses) we plot $4-Z$ as a function of Mason number the RD model with $\phi=0.5$ and $N=4096$. While in the Newtonian regime at high Mason number the contact number is insensitive to Mn, the quantity $4-Z$ decreases ($Z$ increases) as chains form in the magnetically dominated regime. There is an apparent leveling off at the lowest simulated values of Mn, suggesting that $Z$ asymptotes to a value below 4. In order to estimate this value, we plot $Z_{0}-Z$ (Fig. 10b, circles) and adjust the value of $Z_{0}$ to find the cleanest power law at low Mn. For $Z_{0}=3.78$ we find a power law $Z_{0}-Z\sim\text{Mn}^{a}$ with exponent $a\approx 0.37$. Interestingly, a similar scaling relation $Z_{\rm iso}-Z\sim\dot{\gamma}^{0.38}$ has been observed in hard sphere suspensions with no magnetic interactions Lerner et al. (2012). In Fig. 10b we plot the same quantities for the CD model, finding nearly identical values for the extrapolated asymptote $Z_{0}=3.78$ and exponent $a\approx 0.41$. We note that the small difference in the exponent $a$ seems to be entirely due to the $Z$ factor in the definitions of Mn, which differs between the RD and CD model. If we fit both data sets using the same definition of Mn the exponent $a$ is the same for both models within statistical error. We have verified that both $a$ and $Z_{0}$ are independent of $N$ for sufficiently large system sizes, and that their values vary little over a wide range of volume fractions (not shown). Between $\phi=0.3$ and $0.7$ the value of $Z_{0}$ trends from $Z_{0}\approx 3.78$ to $3.85$ and eventually approaches $Z_{0}\approx 4$ as $\phi\rightarrow\phi_{c}$ for both RD and CD models.

To further verify that the microstructure in both models is statistically indistinguishable in the zero Mason number limit, we now investigate the distribution of local contact numbers. In Fig. 10c we plot the fraction of contacts $f_{z}$ having $z$ contacts, for $z=0\ldots 7$, in both the RD and CD models. At large Mason numbers, $f_{z}$ differs strongly between the two models, both in its magnitude and its trend with Mn. However, at low Mn each fraction $f_{z}$ approaches a constant value. To within the accuracy of our measurements, the asymptotes of each $f_{z}$ are equal in the RD and CD models.

To summarize our results on contact number, we have seen that for two types of damping, the flow samples states with the same mean value $Z_{0}$ of the contact number, as well as the same contact number frequencies $\{f_{z}\}_{z=0\ldots 7}$. This provides strong evidence that steady shear flows in the RD and CD models sample the same ensemble of states as Mn tends to zero. However it is also clear that the asymptotically low-Mn regime is at the limit of the lowest Mason numbers we can practically access numerically.

III.2 Cluster statistics

From the snapshots in Fig. 9 it is apparent that the build-up of clusters proceeds differently in the RD and CD models. Here we present evidence that, whereas clustering in the RD model is driven solely by magnetic interactions, inelastic collisions between particles provide a second, unrelated clustering mechanism in the CD model. Clustering due to inelastic collisions is well known in granular gases: particles exit a collision with a lower relative velocity, and hence tend to stay closer togetherGoldhirsch and Zanetti (1993); McNamara and Young (1996). In the CD model, and unlike the RD model, dissipation indeed occurs via collisions. Moreover, due to the model’s overdamped dynamics, particles remain in contact after colliding; i.e. their relative velocity is zero (see Fig 11).

We now seek to quantify the degree of clustering in the RD and CD models. If, as hypothesized above, inelastic clustering is present only in the CD model, one should find differences in, e.g., the time-averaged size $C_{\rm max}$ of the largest cluster in the system. We consider a particle to belong to a cluster if it has a non-zero overlap with any other particle belonging to that cluster. A size $C_{\rm max}=N$ indicates that every particle participates in one cluster. In the left panel of Fig. 12, we plot $C_{\rm max}/N$ in the RD model as a function of Mason number. Note, first, that the data collapse with Mn. Second, there are no clusters of significant size at high values of Mn, when the rheology is Newtonian; however, there is a sharp rise in cluster size below $\text{Mn}\sim 10^{-3}$, coinciding with the magnetically-dominated regime in the flow curve (c.f. Fig. 3). We conclude that “clusters” in the RD model correspond to chains supported by magnetic interactions.

As with the flow curves, the clustering data for the CD model (Fig. 12b) are comparatively complex. First, there is a degree of clustering even in the Newtonian regime. Second, the data do not collapse with Mason number. This clearly indicates the presence of a clustering mechanism independent of magnetic interactions, which we identify with inelastic collisions. Finally, for sufficiently low Mn all particles participate in a single cluster, as in the RD model.

$C_{\rm max}$ also shows qualitatively different dependence on the volume fraction $\phi$ in the two models. In Fig. 12c and d we plot $C_{\rm max}/N$ as a function of $\phi$ at high field strength $H=0.1$. It is clear that the clustering in the RD model shows a much stronger $\phi$ dependence than in the CD model. This $\phi$ dependence is consistent with our previous observations that the Mn needed to reach the plateau in $\sigma$ decreases as $\phi$ is lowered, and that this shift is stronger in the RD model. Here we also include data for $\phi=0.1$ and $0.2$. At these low values of $\phi$, the Mason number needed to reach the yield stress plateau is currently inaccessible in simulation. However there is an increase in $C_{\rm max}$ at low Mn, suggesting that a plateau does emerge at lower Mn. Another way of visualizing the $\phi$ dependence over a wider range of Mn is shown in Fig. 13a and b, where we plot contours of $C_{\rm max}/N$ over the same range of $\phi$ and Mn for the RD and CD model, respectively. Differences are most easily seen by considering, e.g., the $C_{\rm max}/N=0.9$ contour. In the CD model this contour is nearly independent of $\phi$, up to some maximum $\phi$ close to $\phi_{c}$. This suggests that large clusters appear in the CD model at a characteristic Mason number that is independent of $\phi$. In the RD model, by contrast, the value of Mn where clusters appear is an increasing function of $\phi$.

In the snapshots of Fig. 9, it is also evident that the orientation of the emergent chains differs between RD and CD flows. To characterize chain orientation, we study $\theta_{H}$, defined as the average contact angle measured counter-clockwise relative the magnetic field axis (the $\hat{y}$-axis),

[TABLE]

The sum runs over all bonds with a positive overlap. $\theta(\mathbf{u},\mathbf{v})$ is the angle between the vectors $\mathbf{u}$ and $\mathbf{v}$ measured counterclockwise from $\mathbf{v}$ such that $-\pi/2<\theta(\mathbf{u},\mathbf{v})<\pi/2$, giving $0<\theta_{H}<\pi/2$. In Fig. 14 we plot $\sin{2\theta_{H}}$ as a function of Mn for three values of $\phi$. Chains emerge in both models for sufficiently low Mn, indicated by $\sin{2\theta_{H}}\approx 0$. Likewise, at high Mn there is a positive bias, indicating that contacts tend to be rotated in a positive sense with respect to $\bf H$ – as one would expect for collisions due to rapid shear flow. The height of the plateau at high Mn shows stronger $\phi$-dependence in the RD model than in the CD model.

There is a dramatic difference in how the two models cross over between the plateaus at high and low Mn. Whereas $\sin{2\theta_{H}}$ has a sigmoidal shape in the CD model, in the RD model the curve overshoots its low-Mn asymptote. In this intermediate range of Mn, the two models approach their asymptotic values from opposite “directions”: chains in the CD model are rotated counter-clockwise with respect to $\bf H$, while chains in the RD model have a clockwise rotation.

One expects the clusters promoted by inelastic collisions to have a different character from the chain-like structures formed due to magnetic interactions – they should be comparatively compact and isotropic (see Fig. 9). We find the clearest signature of this difference is found by plotting mean number of triangles $\Delta_{s}$ formed by small particles in contact. For a given cluster size, one expects $\Delta_{s}$ to be larger for a compact cluster than for an anisotropic, chain-like structure. $\Delta_{s}$ is plotted in Fig. 15 as a function of Mason number. While $\Delta_{s}$ increases monotonically with decreasing Mn in the RD model, its evolution is non-monotonic in the CD model. There is a peak at intermediate Mn, which we associate with the more compact collisional clusters, followed by a decrease as those clusters are converted to chains.

The data for cluster size, contact angle, and mean triangle number suggest the following picture. In the RD model chain-like clusters build up monotonically as Mn is lowered. In the CD model, in contrast, isotropic clusters form “earlier” (at higher Mn) due to inelastic collisions. As Mn is further lowered and magnetic interactions grow dominant, these compact clusters are reshaped into chains. All relevant observables approach the same asymptotic value in the two models, but they may do so from opposite sides (e.g. $\theta_{H}$ an $\Delta_{s}$). This provides some insight into how the two models’ flow curves can display qualitative differences even as they approach the same asymptote. It also provides indirect support of the hypothesis suggested in the previous section, namely that flow curves approach a finite yield stress plateau at inaccessible values of Mn. Of course one might instead infer that the common asymptote of the RD nor CD flow curves is at zero stress, i.e. that neither has a true yield stress. However this interpretation is disfavored by Occam’s Razor, as all simulated values of $\phi$ show a plateau in the CD flow curve.

IV Conclusions

We have studied the steady state rheology of MR fluids interacting via magnetic, elastic, and two distinct viscous forces. Performing numerical simulations that meet or exceed the lowest values of the Mason number accessed experimentally, we have shown that for moderate volume fractions only systems with contact damping (CD) show a clear plateau in their flow curve. Systems with reservoir damping (RD), by contrast, appear to follow a power law $\sigma\sim\dot{\gamma}^{1-\Delta}$ with $\Delta<1$ – which, if extrapolated to zero strain rate, would imply the absence of a dynamic yield stress. We have argued, instead, that viscous forces must play a subdominant role at asymptotically low $\dot{\gamma}$, and hence either both models possess a yield stress or neither does. The fact that both models display a plateau in their flow curves at sufficiently high volume fractions strongly suggests it is the former: both models possess a dynamic yield stress, with the plateau in the RD flow curve appearing outside the accessible window of Mn for moderate $\phi$. This interpretation is supported by statistical measures of the microstructure, which approach the same asymptote in each model – albeit at the edge of our numerically accessible window in Mn. Cluster statistics suggest that the difference in bulk rheology is related to cluster formation due to inelastic collisions in the CD model, which are absent in RD systems. Despite this conclusion, the clear qualitative difference between the RD and CD flow curves evidenced in our simulations is significant for at least two reasons. First, it persists over a wide interval in Mn including, as previously noted, the lowest values of Mn accessed experimentally. Second, the difference is clearest for moderate values of $\phi$ far below the jamming transition, which are typical of the MR fluids used in applications.

Our work raises several (computationally expensive) questions that might profitably be addressed in future work. One, of course, is whether the speculated crossover to a plateau is in fact seen in RD flow curves at volume fractions around $0.5$ or lower. We have focused on higher $\phi$ values in part to make the connection to jamming, but also because a yield stress, if present, should be more readily apparent. In practice, $\phi$ values around 0.1 are common in experiments and applications. In this dilute limit, chains form, break, and re-form slowly. Hence transients are long and it becomes necessary to simulate for comparatively (and impractically) long total strains.

A second question concerns the role of dimensionality. Inelastic collisions are also present in the CD model (and absent in RD) in higher dimensions, which would suggest that qualitative differences persist. However simulations are needed to determine details such as the apparent value of $\Delta$ and the Mn-interval over which effects are observed.

Third, one can ask about the origins of the exponent $\Delta\approx 0.75$ in the RD flow curves. We note that the critical exponent $\beta$ in directed percolation (DP), which characterizes the mass of the percolating cluster, has a value $\beta\approx 0.276$ in 1+1 dimensions Essam et al. (1988). It is tempting to think there might be a connection to $1-\Delta$ in MR flows, with the applied field defining the time-like dimension. However such a connection is purely speculative.

Finally, one can ask about the role of Coulomb friction, which presumably plays a role in the laboratory. Insofar as Coulomb friction renders collisions between particles inelastic, we expect that shear flows in the CD model more closely resemble systems with friction. While the results and equations in this paper are presented in the context of MR fluids, the model and the findings are more general and we expect that they can be generalized to electrorheological fluids or other similar dipolar systems.

V Acknowledgments

This work was sponsored by NWO Exacte Wetenschappen (Physical Sciences) for the use of supercomputer facilities, with financial support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Netherlands Organization for Scientific Research, NWO). For the final stages of the work D.V. also recevied support from the European Research Council under the European Unions Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement No. 306845. D.V. acknowledges simulating discussions with H.P. Svea.

VI Appendix: Simulation details

VI.1 Dipole moments

There are many effects to consider when modeling the dipole moments $\mathbf{d}$ induced by the external field $\mathbf{H}$. For simplicity we assume our materials are ideal so that we do not need to consider saturation effects at high field strengths. We also assume the magnitude and direction of the induced dipole-moments are given by

[TABLE]

where $\mathbf{H}$ is the the applied magnetic field and $\mathbf{M}$ the magnetization. The dipole moment induced by the external field in a single particle is

[TABLE]

where $V_{\rm ci}$ is the core-volume of particle $i$, and

[TABLE]

The relative permeability of the particles is

[TABLE]

and $\mu_{i}$ is the permeability of the core of particle $i$. The outer shell is assumed to have the same permeability as the carrier fluid.

When there are multiple particles the fields from the induced dipoles interact, giving a total dipole moment of

[TABLE]

where $\mathbf{B}_{ij}$ is given by (10). This is an implicit relation, since $\mathbf{B}$ itself depends on $\mathbf{m}$. Eq. (27) can be solved by iteratively evaluating the expression until it converges Vesely (1977). However we find that for the parameter range investigated here, the correction due to this iterative scheme is negligible, except at the highest field strength we consider. Since our goal for this paper is to reach the lowest Mason numbers possible, and the Mason number is more sensitive to changes in the field strength than to the shear rate, we chose to ignore this effect for all values of $\mathbf{H}$. Consequently all the data presented here are generated using the much faster single particle relation for $\mathbf{m}$ given in equation (24). A major reason why the self-interaction is so low in our system is the core-shell structure of the particles, which prevents the magnetic cores from directly touching each other and ensures the point dipoles remain separated. Note that since $V_{\rm ci}$ and $\beta$ always appear together, these parameters can be varied without changing the result as long as their product stays constant, meaning our results can be mapped to a model where $V_{ci}=V_{i}$ by lowering the value of $\beta$ accordingly.

VI.2 Long range interactions

The dipole-dipole potential between two particles decays as $1/r^{3}$. The interaction is therefore long-ranged in 3D and decays too slowly to be easily truncated in 2D, and care must be taken to correctly include the influence of distant particles. There are several methods to do this, of which the lattice-based Ewald summation Ewald (1921) and cutoff-based reaction field methods Barker and Watts (1973) are the most common – see e.g. Allen and Tildesley (1989); Fennell and Gezelter (2006) for comparisons of different methods. We use a cutoff-based method because it is more computationally efficient (computational complexity $O(N)$) and easier to generalize when changing the geometry of the simulation cell and applying external deformations such as shearing.

While cutoff-based methods are commonly used to simulate MR systems, the long range correction terms used (if any) are rarely published. We therefore include the correction terms employed here. We consider only dipole-dipole interactions; free point charges are not treated. The expressions stated in this Appendix are for 2D systems.

We introduce a cutoff distance $r_{c}$ and evaluate all pair interactions at close distances $r_{ij}<r_{c}$ directly. Evaluating each pair interaction at longer distances quickly becomes computationally expensive. Instead we assume the space outside the sphere given by $r_{c}$ is filled with a uniformly polarized continuous phase. It is then possible to analytically integrate over the continuous phase to obtain the long range correction.

For each observable $\mathcal{O}$ dependent on the dipole potential, it is necessary to calculate a correction term $\mathcal{O}_{LR}=\int_{r_{c}}^{\infty}\tilde{\mathcal{O}}dV$ by integrating the corresponding observable density function $\tilde{\mathcal{O}}$ over $r>r_{c}$. The observable for a single particle $i$ is then given by

[TABLE]

We now show how this is applied to the dipole-dipole potential energy. The magnetic flux density $\mathbf{B}_{j}$ from a dipole $\mathbf{m}_{j}$ at a distance $\mathbf{r}$ is given by

[TABLE]

At short distances $\mathbf{r}$ the local field can be calculated by summing over all particles $j$ located within a sphere of radius $r_{c}$. At longer distances we integrate over the uniformly polarized continuous phase to obtain the long range contribution to the magnetic flux density. In order to perform the integration the discrete particle dipole moment $\mathbf{m}_{i}$ is replaced with the an average dipole moment density $\tilde{\mathbf{m}}$. There are several ways to approximate $\tilde{\mathbf{m}}$; we use

[TABLE]

where we estimate the density of the whole space using the local density. An alternative would be to use the system average or the asymptotic value at infinity (if known) to estimate $\tilde{\mathbf{m}}$. Here we have introduced a weight factor $w(r_{i})$ used to taper the interaction as the cut off distance $r_{c}$ is approached. This prevents discontinuous jumps in measured quantities when particles move in or out of the cutoff sphere. We use a simple linear taper function

[TABLE]

Inserting $\tilde{\mathbf{m}}_{i}$ into (29) and integrating over all $r>r_{c}$ yields the correction term

[TABLE]

The correction to the magnetic potential energy for a given particle $i$ then follows as

[TABLE]

We note that this is an approximation. For a more careful calculation the correction term should be integrated over all space where $w(r)\neq 1$ including the weight function $\mathcal{O}_{LR}=\int_{0}^{\infty}(1-w(r))\tilde{\mathcal{O}}dV$.

It is straight forward to repeat the above procedure for other observables. For the force one obtains

[TABLE]

as expected from symmetry. For the pressure one finds

[TABLE]

and correspondingly for the stress

[TABLE]

While this expression works for isotropic distributions of dipole moments, in our specific case all the dipoles are aligned with the $y$-axis and the correction term is identically zero. We solve this by introducing a second correction term

[TABLE]

where the coefficient

[TABLE]

is a measure of the anisotropy of the packing. This correction term approximates the $\phi$ and Mn dependence over the parameter range we study. However it still assumes that all the dipoles are aligned with the $y$-axis, and it becomes increasingly inaccurate at $\phi<0.3$.

Figure 16 shows the effect of the above mentioned correction terms. In our simulations we use $r_{c}=15r_{m}$ for $\phi>0.3$ and $r_{c}=60r_{m}$ for dilute systems with $0.1<\phi\leq 0.3$. Here $r_{m}$ is the radius of the magnetic core of the larger particles.

In general the the need for corrections is lower for isotropic packings, i.e packings with high Mn or high $\phi$, and their contribution is often insignificant at the $r_{c}$ we use. At the other end in dilute low Mn packings the corrections play an important role as they can reduce the $r_{c}$ needed during simulation. In figure 17 we see the flow curve $\tilde{\sigma}$ vs Mn for $\phi=0.5$ with and without corrections. It is clear from the figure that the corrections are only important at the lowest Mn. The use of the stress correction term shifts the onset of the yield stress plateau to higher Mn, making the plateau easier to observe. However our main conclusions are not sensitive to the use of the correction term; most significantly, our observations regarding the presence or absence of a yield stress plateau at low Mn are also supported by looking at the raw stress without the correction term.