Two-phase volume-averaged predictive theory of dilute ferrofluid spin-up flow in a rotating magnetic field

Abstract

We have developed a parameter-free, two-phase, volume-averaged approach to predictively describe the spin-up flow of dilute, cluster-free ferrofluids excited by low-frequency rotating magnetic fields. Predictive validation of the model was performed through a thorough comparison with local velocity profile measurements, and it demonstrated its ability to capture the spin-up flow dynamics without the need for parameter tuning by carefully delineating the validity domain of the ferrofluid dilutedness conditions. To gain insight into the underlying flow mechanisms, we performed a systematic parametric analysis examining the effects of the induced magnetic field, the dipolar interactions between magnetic nanoparticles and the demagnetizing field. How these mechanisms shape the flow of dilute ferrofluids excited by low-frequency rotating fields in a standard spin-up flow geometry has been addressed using probabilistic nanoparticle orientational dynamics, combining Faxén's laws and the Smoluchowski equation to describe the transport of particle magnetic moments. Our findings revealed that the induced magnetic field is the primary driving force of ferrofluid spin-up flow. The dipole interactions and demagnetizing field, on the other hand, contribute only as secondary phenomena to the overall flow behaviour. Furthermore, we have discussed the potential extension of the two-phase approach, in particular with respect to the formation of chain-like aggregates under the influence of strong magnetic fields. Overall, our study provides valuable insights into the complex dynamics of ferrofluid flow and contributes to a comprehensive understanding of the key mechanisms governing the spin-up flow of dilute ferrofluids excited by low-frequency rotating magnetic fields.

1. Introduction

Since Gilbert's 17th-century seminal treatise ‘De Magnete’, the fascination with magnetism has never waned. The advent of synthetically stabilized magnetic colloidal suspensions, known as ferrofluids, has expanded the frontiers of magnetism. Supported by extensive theoretical studies and experimental validation, ferrofluids have ventured into unexpected territories, opening up new applications, such as micro/nanofluidics (Varma et al. 2016; Liu et al. 2020), microdevice mixing (Yi, Qian & Bau 2002; Boroun & Larachi 2016) and biomedical advances (Kose et al. 2009). This widespread interest stems from the intriguing behaviour of ferrofluids when exposed to a magnetic field. Magnetic particles naturally align themselves in the direction of the field, creating a uniquely magnetized collective state. This alignment allows precise control of the ferrofluid and results in unique flow patterns.

Magnetic fields, both oscillating and rotating, are often used to manipulate ferrofluids. Oscillating fields lead to phenomena such as magnetoviscous effects (Felderhof 2000), while rotating fields induce spin-up flows. These flows arise from confining motionless ferrofluids within cylindrical enclosures, where the magnetic field imparts a torque to the magnetic entities in the fluid, causing them to spin. This angular momentum is then converted to linear momentum, resulting in ferrofluid rotation. The pioneering experimental work of Moskowitz & Rosensweig (1967) was limited to open cylinders, studying only the air–ferrofluid interface. Later research by Rosensweig, Popplewell & Johnston (1990) revealed that understanding bulk spin-up flow requires more than interface measurements due to hydrodynamic interactions between the ferrofluid and the surrounding air. In particular, the ferrofluid corotates with the magnetic field when the interface is concave, while it counter-rotates when the interface is convex. Zaitsev & Shliomis (1969) developed a theory that described the spin-up flow by treating the ferrofluid as a homogeneous magnetic fluid. This theory was based on the continuum approach proposed by Condiff & Dahler (1964) for a pseudohomogeneous fluid flow with molecular rotations as degrees of freedom. It included terms in the Navier–Stokes equation to account for translational motion, as well as the effects of antisymmetric stresses and magnetic forces. It also introduced internal spin diffusion, described by the spin viscosity parameter $\eta '$, which forms the basis of spin diffusion theory. Spin diffusion and antisymmetric stress are the driving forces behind convective flow in ferrofluids.

Zaitsev & Shliomis (1969) estimated spin viscosity $\eta '$ at approximately $10^{-21}\ {\rm kg}\ {\rm m}\ {\rm s}^{-1}$, according to a proportional relationship with the square of the particle radius, $\eta ^{'} \propto R_{p}^2$. Qualitatively, the velocity profile predicted with this value of $\eta '$ corresponds to that obtained experimentally by Rosensweig et al. (1990). However, there are discrepancies in the prediction of diffusive layer thickness and velocity amplitude, which cast doubt on the validity of $\eta '$. Spin diffusion theory, developed primarily to describe bulk ferrofluid flow, has limitations in explaining surface-dominated effects, as acknowledged by Rosensweig et al. (1990). These authors observed that surface effects play a significant role, and spin diffusion theory alone could not account for these phenomena. In response to earlier observations, Glazov (1975) provided a theoretical explanation within the framework of spin diffusion theory, demonstrating that generating spin-up flow with a uniform magnetic field is impossible. However, in the experimental set-up of Rosensweig et al. (1990) the magnetic field was kept uniform and local measurements of the velocity field confirmed the presence of macroscopic flow. Other explanations have explored thermal effects (Shliomis, Lyubimova & Lyubimov 1988; Pshenichnikov, Lebedev & Shliomis 2000) and wall-slip boundary conditions of spinning particles (Kaloni 1992) as potential drivers of the ferrofluid macroscopic flow. However, these hypotheses have yet to be conclusively validated by direct comparison with experimental data.

Local velocity measurements within the spin-up flow geometry are crucial to validate the value of $\eta '$ proposed by Zaitsev & Shliomis (1969) and to determine the direction of bulk ferrofluid rotation. Chaves et al. (2006) performed an experiment using an ultrasound-based technique that allowed them to obtain local measurements of the azimuthal velocity component generated by a rotating magnetic field in a cylindrical device. They found that the ferrofluid bulk rotates with the magnetic field, while the counter-rotation at the interface is due to surface effects. This finding is consistent with spin diffusion theory, although there are some qualitative differences in the velocity profile, such as the velocity amplitude and the thickness of the diffusion layer, which covers approximately $30\,\%$ of the tube radius. Consequently, the spin viscosity value proposed by Zaitsev & Shliomis (1969) would significantly underestimate the presumed hypothetical numerical value required to accurately predict the experimental velocity profile measured by Chaves et al. (2006). Building on this observation, Chaves, Zahn & Rinaldi (2008) performed both an experimental investigation and an asymptotic analysis of the spin diffusion theory using the regular perturbation method. By adjusting the spin viscosity, they were able to achieve quantitative agreement between the results of the asymptotic analysis and the experimental velocity measurements. This analysis revealed that the value of $\eta '$ must be in the range of $10^{-12}$ to $10^{-8}\ {\rm kg}\ {\rm m}\ {\rm s}^{-1}$ to accurately predict the observed velocity profiles. Notably, these $\eta '$ values are orders of magnitude (10 to 12 orders of magnitude) higher than the $\eta '$ value proposed by Zaitsev & Shliomis (1969). In a theoretical study by Finlayson (2013), the influence of the spin viscosity and the Langevin equation on spin-up flow was thoroughly investigated. Flow ceases as $\eta '$ approaches zero, and in the presence of an inhomogeneous magnetic field, an irregular flow occurs, which contradicts the experimental results. Shliomis (2021) conducted a theoretical investigation, expressing the spin viscosity as a function of the particle moment of inertia, which at a very small value could not explain the bulk ferrofluid motion. As a result, the spin diffusion theory was linked to the dissipative effects of rotating particles, resulting in heat dissipation, which generates magnetic field inhomogeneities, thereby inducing ferrofluid flow. However, the interfacial flow was explained by assuming $\eta '=0$. The direction and velocity of ferrofluid rotation were described as functions of contact angle, meniscus height, magnetic field amplitude and frequency. Remarkably, the theoretical results were in good agreement with the experiment performed by Rosensweig et al. (1990). The determination of a physically sound value of $\eta '$ remains a subject of debate since direct measurement of this parameter is not feasible. For the time being, studies are limited to the empirical determination of this parameter, which involves its adjustment from experimental results and reduces spin diffusion theory to a descriptive rather than a predictive method. The lack of fundamental understanding, due to the absence of a predictive model has limited the studies to experimental investigations where different experiments, including open and annular cylindrical settings (Chaves, Torres-Diaz & Rinaldi 2010; Torres-Díaz & Rinaldi 2011) and spherical configurations (Torres-Díaz et al. 2012), have consistently shown macroscopic flow within ferrofluids. Of particular interest is the experiment performed by Torres-Diaz et al. (2014), which focused on very dilute ferrofluids with volumetric concentrations below $1\,\%$.

The literature highlights two main issues: the need to assign a physical value to $\eta '$ to explain ferrofluid bulk flow; and the suitability of the Condiff & Dahler (1964) continuum approach for treating ferrofluid as a pseudohomogeneous fluid, which may not accurately represent the discrete nature of ferrofluid particles. In this framework, the particles can transfer their internal angular momentum by a process similar to molecular diffusion. However, from a phenomenological standpoint, this mechanism is implausible because the particles form a discrete phase and remain distinctly apart from one another. Therefore, it becomes imperative to describe the rotational dynamics of the particles while preserving their discrete nature. Similarly, the magnetization transport equation is defined for the entire ferrofluid, whereas only the particles are magnetized. Consequently, a two-phase approach that distinguishes the magnetic particles from the non-magnetic matrix may provide a more appropriate description. This perspective is consistent with the colloidal suspension theory, as discussed by Batchelor (1977), which separates the particle and liquid contributions in the bulk stress analysis. Saintillan & Shelley (2008a,b) developed a kinetic theory to describe the dynamics of a suspension of self-propelled particles with active stresslets. In their approach, a homogeneous fluid description was used for the linear momentum transport in the suspension, while the particles were characterized by a probability density governed by the Smoluchowski equation, which includes the relevant mechanisms for the active particles. In summary, reconsidering the assignment of a physical value to $\eta '$ for bulk ferrofluid flow and adopting a two-phase approach that distinguishes between the magnetic dispersed phase and the non-magnetic continuous phase may provide a more accurate representation of ferrofluid behaviour. This change in perspective, viewing the ferrofluid as a colloidal suspension rather than a pseudohomogeneous fluid, can elucidate the coupling between the dispersed magnetic particles and the continuous liquid phase.

In this study, we present a two-phase fully predictive model specifically designed for the flow of a dilute clusterless ferromagnetic colloidal suspension under the influence of a rotating magnetic field. The model incorporates the equations for the conservation of mass, linear and angular momentum, and the induced magnetic field at the continuum scale using the average volume theorem (AVT). In a manner similar to Saintillan & Shelley (2008a,b) characterization of active stresslets using the second-order (nematic-order) moment of the Smoluchowski equation, the two-phase equations governing ferrofluid behaviour are then coupled to an adapted Smoluchowski equation tailored to the context of ferromagnetic suspensions. This adaptation focuses on the zeroth- and first-order moments of the Smoluchowski equation, which effectively capture the transport dynamics and the average orientation of the magnetic moments of the individual magnetic nanoparticles within the clusterless ferrofluid. The proposed model is applied in the context of spin-up flow in a cylindrical geometry, aiming to improve our fundamental understanding of the process and to identify the genuine driving mechanism behind the rotating flow in the ferrofluid bulk. To establish the model's validity, it will be compared with the experimental results obtained by Torres-Diaz et al. (2014). This validation process will help to define the range of applicability of the model. Furthermore, a comprehensive parametric study will be performed to quantify the individual contributions of the induced magnetic field, particle transport, dipole–dipole interactions (DDIs) and the demagnetizing field. This analysis will allow us to identify the key parameter(s) responsible for the spin-up flow of the ferrofluid, providing valuable insights into the underlying mechanisms at play. Finally, we will discuss theoretical results that go beyond the scope of the model, highlighting potential avenues for further investigation and development of the two-phase approach. Special attention will be given to non-dilute ferrofluids, where clustering becomes significant and cannot be ignored.

2. Ferrofluid spin-up flow two-phase formulation

Spin-up flow refers to a macroscopic convective rotational motion elicited in an otherwise quiescent magnetic nanofluid (ferrofluid) by a rotating external magnetic field. To illustrate this process, a suspension of magnetic nanoparticles in a Newtonian liquid is confined in a cylindrical enclosure (radius $R$, length $L$) and subjected to a uniform external magnetic field (frequency $f_0$, intensity $H_0$) applied perpendicular to the $z$-axis of the tube (figure 1a). The magnetic inclusions in the suspension have the unique ability to ‘bodily’ align their permanent magnetic moment vector with the constantly moving magnetic field. Hydrodynamic interaction with the surrounding liquid prevents the discrete magnetic inclusions from achieving perfect alignment with the magnetic field. It is this misalignment, measured as a lag angle, that causes the magnetic inclusions to experience a driving torque at the origin of the macroscopic flow (figure 1b).

Figure 1. Schematic of a ferrofluid spin-up flow geometry in a cylindrical enclosure (a) and a zoomed-in view of a single magnetic nanoparticle (b) subjected to a magnetic field. Here $\boldsymbol {m}$ refers to the orientation of the magnetization vector $\boldsymbol {M}$. Spin-up flow is generated by rotating an external constant modulus magnetic field perpendicular to the axis of symmetry of the cylindrical cell containing the ferrofluid (see (5.1)).

In contrast to the standard ferrohydrodynamic model (Zaitsev & Shliomis 1969; Rosensweig 2013), which cannot predict spin-up flows (Finlayson 2013; Torres-Diaz et al. 2014; Shliomis 2021), a comprehensive theoretical representation must account for the two-phase nature of the magnetic suspension. At the very least, this perspective must necessarily include the resolution of the linear and angular momentum equations corresponding to each of the two phases of the magnetic suspension. We propose to tackle the problem by appealing to the AVT, which has been extensively used to describe transport phenomena in two-phase flows and in porous media in particular (Gray 1975; Whitaker 1999; Drew & Passman 2006). The AVT approach models transport phenomena separately in each phase using volume-averaged equations, which are then coupled by jump conditions at interfaces. This makes it possible to analyse the linear and angular momentum transport under the proper distribution of the induced magnetic field within each of the suspension's phases. A notable practical advantage of this method is that it provides average values for unknowns, such as the linear and angular velocities of the magnetic (discrete) phase at each point in the ferrofluid domain. While a macroscopic two-phase description may suffice for non-colloidal suspensions, colloidal ferrofluid particles require the incorporation of Brownian motion into the transport processes, leading to a probabilistic treatment. This coupling can be achieved using the Smoluchowski equation, focusing on its zeroth- and first-order moments, which take into account hydrodynamic conditions and the effect of Brownian motion on the average transport of particles and the orientation of their magnetic moments. Section 4 will go over the connection of the two-phase ferrofluid model with the zeroth- and first-order moments of the Smoluchowski equation.

Consider a two-phase suspension consisting of a liquid, denoted by the subscript ${l}$, and a dispersed magnetic phase (particles), denoted by the subscript ${p}$. Figure 2 is a schematic representation of the ferrofluid system parsed into representative elementary volumes (REV), over which the transport equations of the two phases are averaged. Such two phases (${l}$ and ${p}$) occupy the volumes $V_{l},V_{p}$, respectively, in the REV volume $V$, whose perimeter is bounded by the surface $\varGamma =A_{p-e} \cup A_{{l}-{e}}$. The interface separating the two phases inside the REV is denoted by $A_{{l}-{p}}$ with the normal orientational unit vectors $\boldsymbol {n}_{l}=-\boldsymbol {n}_{p}$. Let $\boldsymbol {x}$ be the macroscopic position vector characterizing the centre of the elementary volume with respect to the macroscopic system, and $\boldsymbol {y}$ the microscopic position vector on the scale of the REV, defined with respect to its centroid. For example, consider a scalar field $\varPhi _{l}$ transported in phase $\rm {l}$, defined with respect to position $\boldsymbol {x}+\boldsymbol {y}$. Using AVT, the mean field $\langle \varPhi _{l} \rangle$ with respect to REV is given by (Gray 1975; Whitaker 1999)

(2.1)\begin{equation} \langle \varPhi_{l}\rangle|_{\boldsymbol{x}} = \frac{1}{V} \int_V \chi_{l} \varPhi_{l}|_{\boldsymbol{x}+\boldsymbol{y}}\,{\rm d} V, \end{equation}

where $\chi _{l}$ denotes the phase indicator describing the spatial distribution of phase ${l}$ in the REV. The latter is defined by

(2.2)\begin{equation} \left. \begin{array}{ll@{}} \chi_{l} = 1 & \boldsymbol{x}+\boldsymbol{y} \in V_{l} \cup A_{{l}-{p}}, \\ 0 & \text{elsewhere}. \end{array} \right\} \end{equation}

Figure 2. Schematic description of a REV of a suspension composed of a liquid and rigid spherical particles.

Supplementary materials (§ 1) available at https://doi.org/10.1017/jfm.2024.32 contain the developments according to the AVT approach leading to the formulation of the mean transport equation for the so-called quantity $\varPhi _{l}$.

2.1. Ferrofluid phase-specific microscale transport equations

To establish average macroscopic transport equations governing the ferrofluid behaviour at each macroscopic position $\boldsymbol{x}$ within the REV, one must first formulate microscopic transport equations for the two suspension phases at each microscopic position $\boldsymbol{y}$ within the REV. These equations involve the conservation of mass, linear and angular momentum and the application of Maxwell's equations to describe the magnetic fields in each ferrofluid phase.

The microscale equation of the incompressible liquid phase (2.3) represents the mass conservation equation of the liquid as a continuum:

(2.3)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (\chi_{l} \boldsymbol{v}_{l}) = 0. \end{equation}

While the liquid phase at the microscale (represented by points $\boldsymbol {y}$ in figure 2) can be described as a continuum using a microscopic mass conservation equation (see (2.3)), this approach does not hold for the discrete magnetic phase. To address this, we establish mass conservation for individual particles within the REV by relying on the mass invariance of a single particle. Specifically, we relate the mass of an individual particle (i) to its density $\rho _{p}$, particle volume $V_{p}$ and phase indicator $\chi ^i_{p}$:

(2.4)\begin{equation} m^i_{p} = \int_{V_{p}} \rho_{p}\,{\rm d} V = \chi_{p}^i \rho_{p} V_{p}, \end{equation}

where $V_{p}$ is the volume of a single particle.

Assuming constant physicochemical properties, the linear momentum conservation equation for the liquid phase at the microscopic scale is

(2.5)\begin{equation} \rho_{l} \left(\frac{\partial \chi_{l} \boldsymbol{v}_{l}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\chi_{l} \boldsymbol{v}_{l} \otimes \boldsymbol{v}_{l})\right) = \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\mathcal{T}}}_{\!l} + \boldsymbol{F}^{e}_{l}, \end{equation}

where ${\boldsymbol{\mathcal{T}}}_{\!l}$ is the viscous-pressure stress tensor containing the symmetric part characterizing shear stress. The term $\boldsymbol {F}^{e}_{l}$ refers to an external volume force.

For a Newtonian and non-polar liquid, the stress tensor, ${\boldsymbol{\mathcal{T}}}_{\!l}$, is given by the following relation:

(2.6)\begin{equation} {\boldsymbol{\mathcal{T}}}_{\!l} =-\chi_{l} P \boldsymbol{\mathsf{I}}+ \mu ({\boldsymbol{\nabla}} \chi_{l} \boldsymbol{v}_{l} + ({\boldsymbol{\nabla}} \chi_{l} \boldsymbol{v}_{l})^{t} -\tfrac{2}{3} \boldsymbol{\nabla} \boldsymbol{\cdot} (\chi_{l} \boldsymbol{v}_{l}) \boldsymbol{\mathsf{I}}) +\mu_b \boldsymbol{\nabla} \boldsymbol{\cdot} (\chi_{l} \boldsymbol{v}_{l}) \boldsymbol{\mathsf{I}}, \end{equation}

where the dynamic viscosity, $\mu$, of the liquid is assumed constant. $\mu _{b}$ is the bulk viscosity associated with expansion (compression) processes.

Taking into account that the bulk viscosity, $\mu _b$, is often negligible in the case of fluids, and based on the liquid continuity (2.3) which describes the incompressibility of the fluid with zero divergence of the velocity field, the stress tensor (2.6), reduces to

(2.7)\begin{equation} {\boldsymbol{\mathcal{T}}}_{\!l} =-\chi_{l} P \boldsymbol{\mathsf{I}}+ \mu ({\boldsymbol{\nabla}} \chi_{l} \boldsymbol{v}_{l} + ({\boldsymbol{\nabla}} \chi_{l} \boldsymbol{v}_{l})^{t}). \end{equation}

By substituting (2.7) in (2.5), the liquid's linear momentum equation becomes

(2.8)\begin{equation} \rho_{l} \left(\frac{\partial \chi_{l} \boldsymbol{v}_{l}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\chi_{l} \boldsymbol{v}_{l} \otimes \boldsymbol{v}_{l} )\right) =-\boldsymbol{\nabla} \chi_l P + \mu \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{\nabla}} \chi_l \boldsymbol{v}_l + \boldsymbol{F}^{e}_{l}. \end{equation}

It is important to emphasize that at the microscopic scale within the REV, where the particles are individually spaced and do not form a continuous medium, the equation governing the angular momentum of the particles does not have a conventional continuum relationship. Therefore, the force balance on a single particle $i$ is used as a microscopic equation as a prerequisite to derive a macroscopic continuum equation for all particles in the REV. The momentum conservation equation for a Brownian particle $i$, also known as the Langevin equation, is expressed as follows:

(2.9)\begin{equation} \rho_{p} \left.\frac{{\rm d} \boldsymbol{v}_{p}^i}{{\rm d} t}\right|_{\boldsymbol{y}=\boldsymbol{y}_i} = \boldsymbol{F}_i^{lp} + \boldsymbol{F}_i^{pp} + \boldsymbol{F}_i^{e}, \end{equation}

where $\boldsymbol {v}_{p}^i(\kern0.7pt \boldsymbol {y}_i)$ is the velocity of the particle $i$ evaluated at its centre of mass $\boldsymbol {y}=\boldsymbol {y}_i$. Here $\boldsymbol {F}_i$ denotes an external force per unit volume. The subscripts $lp,pp$ refer to the liquid–particle and particle–particle interaction forces, respectively.

The conservation equation of angular momentum at the microscopic scale is modelled similarly to the conservation equations of mass and linear momentum of particles by the angular Langevin equation of a particle $i$,

(2.10)\begin{equation} I_{p} \left.\frac{{\rm d} \boldsymbol{w}_{p}^i}{{\rm d} t}\right |_{\boldsymbol{y}=\boldsymbol{y}_i} = \boldsymbol{\varGamma}_i^{lp} + \boldsymbol{\varGamma}_i^{pp} + \boldsymbol{\varGamma}_i^{e}, \end{equation}

where $\boldsymbol {w}_{p}^i(\kern0.7pt \boldsymbol {y}_i)$ is the angular velocity of the particle $i$ evaluated at $\boldsymbol {y}=\boldsymbol {y}_i$; $I_{p}= \frac {2}{5} \rho _{p} R_{p}^2$ is the moment of inertia of the particle $i$ per unit volume of a particle; $\boldsymbol {\varGamma }_i$ denotes the force moment per unit volume of the particle $i$; $\boldsymbol {\varGamma }_i^{e}$ denotes the external moment, per unit volume, exerted on the particle $i$; and $\boldsymbol {\varGamma }_i^{lp}$ and $\boldsymbol {\varGamma }_i^{pp}$ are the force moments of liquid–particle and particle–particle interactions, respectively.

When an external magnetic field is applied to a ferromagnetic suspension, an induced magnetic field is created within its constitutive liquid and nanoparticle phases. While this induced field affects both phases, the magnetization process itself, which involves the alignment magnetic moments, is exclusive to the ferromagnetic particles. Thus, a magnetic torque is exerted only on the nanoparticles. Consequently, the microscopic modelling of the angular momentum considers only the discrete magnetic phase while explicit modelling of the liquid angular momentum is unnecessary. Assuming that the ferrofluid consists of non-conducting discrete magnetic inclusions and a liquid, Maxwell's equations can be simplified to Maxwell–Ampère and Maxwell-flux equations without charge displacement or electric current. At the microscopic level, the Maxwell–Ampère equations for ${l}$ and ${p}$ are as follows:

(2.11)$$\begin{gather} \boldsymbol{\nabla} \times (\chi_{l} \boldsymbol{H}_{l}) = 0 \quad \in V_{l}, \end{gather}$$

(2.12)$$\begin{gather}\boldsymbol{\nabla} \times (\chi_{p} \boldsymbol{H}_{p}) = 0 \quad \in V_{p}, \end{gather}$$

where $\boldsymbol {H}$ is the magnetic field vector expressed in ${A/m}$ and $\mu _0= 4 \times {\rm \pi}10^{-7}\ {\rm N}/{\rm A}^{-2}$ is the vacuum magnetic permeability.

While the Maxwell-flux equations of ${l}$ and ${p}$ read as follows:

(2.13)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} (\chi_{l} \boldsymbol{B}_{l}) = 0 \quad \in V_{l}, \end{gather}$$

(2.14)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} (\chi_{p} \boldsymbol{B}_{p}) = 0 \quad \in V_{p}, \end{gather}$$

where $\boldsymbol {B}$ is the magnetic induction vector which can be formally expressed in terms of the magnetization vector, $\boldsymbol {M}$, by the following relation:

(2.15)\begin{equation} \boldsymbol{B} = \mu_0 (\boldsymbol{H} + \boldsymbol{M}). \end{equation}

Liquids constituting ferrofluids are in general linear, homogeneous and isotropic media, where the magnetization is directly proportional to the magnetic field vector, such as

(2.16)\begin{equation} \boldsymbol{M}_{l} = \varphi \boldsymbol{H}_{l}, \end{equation}

where the dimensionless proportionality constant is called magnetic susceptibility. Liquids are generally weakly magnetized where the value of $\varphi$ is very small. Furthermore, $\varphi$ is positive if the material is paramagnetic and negative if the material is diamagnetic. Using the magnetization expression (2.16), the magnetic induction vector $\boldsymbol {B}_{l}$ in the liquid phase is written as follows:

(2.17)\begin{equation} \boldsymbol{B}_{l} = \mu_0 (1 + \varphi_{l}) \boldsymbol{H}_{l}. \end{equation}

In the case where the liquid is non-magnetic, $|\varphi _{l}| <<1$, the expression for the magnetic induction reduces to

(2.18)\begin{equation} \boldsymbol{B}_{l} = \mu_0 \boldsymbol{H}_{l}. \end{equation}

Compared with the liquid, the particles are superparamagnetic because they are considered monodomain with a permanent magnetic moment $\boldsymbol {M}_{p}=M_{d}\boldsymbol {m}$. The magnetic induction vector of the particles in this case can be written as

(2.19)\begin{equation} \boldsymbol{B}_{p} = \mu_0 (\boldsymbol{H}_{p} + M_{d} \boldsymbol{m}), \end{equation}

where $M_{d}$ refers to the domain magnetization of the particles and $\boldsymbol {m}$ denotes the orientation of particles’ magnetic moments (figure 1b). It should be noted that the model governing the orientation dynamics of the particles’ magnetic moments, $\boldsymbol {m}$, with respect to the magnetic field, $\boldsymbol {H}$, is described in detail in § 4.

Taking into account the expressions of the induction vectors of both phases (2.18), (2.19), the Maxwell-flux equations read as follows:

(2.20)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} (\chi_{l} \boldsymbol{H}_{l}) = 0 \quad \in V_{l}, \end{gather}$$

(2.21)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \chi_{p} (\boldsymbol{H}_{p} + M_{d} \boldsymbol{m}) = 0 \quad \in V_{p}. \end{gather}$$

2.2. Ferrofluid macroscale two-phase transport equations

To establish the macroscale transport equations for a two-phase ferrofluid, a methodology similar to that used in the supplementary materials (§ 2) to derive the hypothetical scalar field transport equation $\varPhi _{l}$ can be applied. This approach involves averaging each transport equation over the REV and using Gray's decomposition (Gray 1975). The derivation of the averaged transport equations is also presented in detail in the supplementary materials.

Applying the AVT to the liquid mass conservation equation at the microscopic scale, (2.3), yields

(2.22)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{l} \langle \boldsymbol{v}_{l} \rangle^{l}) + \frac{1}{V} \int_{A_{{l}-{p}}} \boldsymbol{v}_{l} \boldsymbol{\cdot} \boldsymbol{n}_{l} \,{\rm d} A = 0, \end{equation}

where $\epsilon _{l}$ is the volume fraction of the liquid phase.

The average conservation equation of the mass of the particles, based on the invariability of the total volume of the particles, is written as

(2.23)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{p} \langle \boldsymbol{v}_{p} \rangle^{p}) + \frac{1}{V} \int_{A_{{l}-{p}}} \boldsymbol{v}_{p} \boldsymbol{\cdot} \boldsymbol{n}_{p} \,{\rm d} A = 0 \end{equation}

where $\epsilon _{p}$ is the volume fraction of the magnetic nanoparticles.

Likewise, applying AVT to (2.8), which characterizes the liquid's linear momentum balance equation at the microscopic scale, yields the following macroscopic liquid linear momentum balance equation:

(2.24)\begin{align} & \rho_{l} \left(\frac{\partial \langle \boldsymbol{v}_{l} \rangle}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{l} \langle \boldsymbol{v}_{l} \rangle^{l} \otimes \langle \boldsymbol{v}_{l} \rangle^{l}) + \boldsymbol{\nabla}\boldsymbol{\cdot} (\epsilon_{l} \langle \boldsymbol{\hat{v}}_{l} \otimes \boldsymbol{\hat{v}}_{l}\rangle^{l})\right) \nonumber\\ &\quad =-\epsilon_{l} \boldsymbol{\nabla} \langle P \rangle^{l} + \mu (\epsilon_{l} \nabla^2 \langle \boldsymbol{v}_{l} \rangle^{l} + {\boldsymbol{\nabla}} \langle \boldsymbol{v}_{l} \rangle^{l} \boldsymbol{\cdot} \boldsymbol{\nabla} \epsilon_{l} + \langle \boldsymbol{v}_{l} \rangle^{l} \nabla^2 \epsilon_{l}) + \langle \boldsymbol{F}^{e}_{l}\rangle \nonumber\\ &\qquad +\frac{1}{V} \int_{A_{{l}-{p}}} (-\chi_{l}\hat{P} \boldsymbol{\mathsf{I}} + \mu {\boldsymbol{\nabla}} \chi_{l} \boldsymbol{\hat{v}}_{l})\boldsymbol{\cdot} \boldsymbol{n}_{l} \,{\rm d} A . \end{align}

The Langevin microscale (2.9) and (2.10) once averaged on the REV by the same approach, allow us to express their counterpart for the macroscale concerning the conservation of the linear and angular momentum of the discrete magnetic phase, respectively,

(2.25)\begin{align} & \rho_{p} \left(\frac{\partial \langle \boldsymbol{v}_{p} \rangle}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ( \epsilon_{p} \langle \boldsymbol{v}_{p} \rangle^{p} \otimes \langle \boldsymbol{v}_{p} \rangle^{p}) + \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{p} \langle \boldsymbol{\hat{v}}_{p} \otimes \boldsymbol{\hat{v}}_{p}\rangle^{p})\right) \nonumber\\ &\quad =\langle \boldsymbol{F}^{lp} \rangle + \langle \boldsymbol{F}^{pp} \rangle + \langle \boldsymbol{F}^{e}_{p} \rangle, \end{align}

(2.26)\begin{align} & I_{p} \left(\frac{\partial \langle \boldsymbol{w}_{p} \rangle}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\epsilon_{p} \langle \boldsymbol{v}_{p} \rangle^{p} \otimes \langle \boldsymbol{w}_{p} \rangle^{p}) + \boldsymbol{\nabla}\boldsymbol{\cdot} \epsilon_{p}(\langle \boldsymbol{\hat{v}}_{p} \otimes \boldsymbol{\hat{w}}_{p}\rangle^{p})\right) \nonumber\\ &\quad =\langle \boldsymbol{\varGamma}^{lp} \rangle + \langle \boldsymbol{\varGamma}^{pp} \rangle + \langle \boldsymbol{\varGamma}_{p}^{e} \rangle. \end{align}

The average Ampère–Maxwell equations for the liquid and particles on the REV are written as follows:

(2.27)$$\begin{gather} \boldsymbol{\nabla} \times (\epsilon_{l} \langle \boldsymbol{H}_{l} \rangle^{l})+ \frac{1}{V} \int_{A_{{l}-{p}}} \chi_{l} \boldsymbol{H}_{l} \times \boldsymbol{n}_{l} \,{\rm d} A = 0, \end{gather}$$

(2.28)$$\begin{gather}\boldsymbol{\nabla} \times (\epsilon_{p} \langle \boldsymbol{H}_{p}\rangle^{p}) + \frac{1}{V} \int_{A_{{l}-{p}}} \chi_{p} \boldsymbol{H}_{p} \times \boldsymbol{n}_{p} \,{\rm d} A = 0. \end{gather}$$

Combining (2.27) and (2.28) yields

(2.29)\begin{equation} \boldsymbol{\nabla} \times \langle \boldsymbol{H} \rangle + \frac{1}{V} \int_{A_{{l}-{p}}} (\chi_{p} \boldsymbol{H}_{p} - \chi_{l} \boldsymbol{H}_{l}) \times \boldsymbol{n}_{p}\,{\rm d} A = 0, \end{equation}

where $\langle \boldsymbol {H}\rangle$ is the averaged magnetic field over the suspension,

(2.30)\begin{equation} \langle \boldsymbol{H}\rangle = \epsilon_{l} \langle \boldsymbol{H}_{l} \rangle^{l} + \epsilon_{p} \langle \boldsymbol{H}_{l} \rangle^{p}.\end{equation}

The average Ampère-flux equations, for liquid and particle phases, are written as follows:

(2.31)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{l} \langle \boldsymbol{H}_{l} \rangle^{l}) + \frac{1}{V} \int_{A_{{l}-{p}}} \chi_{l} \boldsymbol{B}_{l} \boldsymbol{\cdot} \boldsymbol{n}_{l} = 0, \end{gather}$$

(2.32)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{l} \langle \boldsymbol{H}_{l} \rangle^{p} + M_{d} \langle \boldsymbol{m} \rangle) + \frac{1}{V} \int_{A_{{l}-{p}}} \chi_{p} \boldsymbol{B}_{p} \boldsymbol{\cdot} \boldsymbol{n}_{p} \,{\rm d} A = 0. \end{gather}$$

Combining (2.31) and (2.32) similarly yields the average Maxwell-flux equation for the suspension:

(2.33)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (\langle \boldsymbol{H} \rangle + M_{d} \langle \boldsymbol{m} \rangle) + \frac{1}{V} \int_{A_{{l}-{p}}} (\chi_{p} \boldsymbol{B}_{p} - \chi_{l} \boldsymbol{B}_{l}) \boldsymbol{\cdot} \boldsymbol{n}_{p} \,{\rm d} A =0. \end{equation}

3. Closure equations

To transform from the microscopic to the macroscopic formulation using AVT, closure relations in the form of additional interfacial transfer terms must be determined and incorporated into (2.24) to (2.26), (2.29) and (2.33) to close the system of equations in the newly derived macroscopic two-phase ferrohydrodynamic formulation. These relations depend on local flow conditions, such as phase velocities, volume fractions and interfacial area. The resulting macroscopic mass, linear and angular momentum equations and Maxwell's equations will then be expressed in terms of the average fields, accounting for the physical phenomena associated with the exchange between the two phases.

3.1. Mass jump conditions

The liquid and particle phases are mutually impermeable. In this case the normal and tangential velocities of the liquid and the particles at the interfaces must be zero,

(3.1) \begin{equation} \left.\begin{array}{l@{}} \boldsymbol{v}_{l} \boldsymbol{\cdot} \boldsymbol{n}_{p} = \boldsymbol{v}_{p}\boldsymbol{\cdot} \boldsymbol{n}_{p} = 0 ,\\ \boldsymbol{v}_{l} \times \boldsymbol{n}_{p} = \boldsymbol{v}_{p} \times \boldsymbol{n}_{p} = 0. \end{array}\right\} \end{equation}

Taking into account the jump condition (3.1) and the incompressibility of the two phases, the mass conservation equations for the liquid and particles can be expressed as follows:

(3.2)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \langle \boldsymbol{v}_{l} \rangle = 0, \end{gather}$$

(3.3)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \langle \boldsymbol{v}_{p} \rangle = 0. \end{gather}$$

3.2. Liquid linear momentum jump condition

In an investigation by Frank et al. (2003), both experimental and theoretical analyses were performed to study the behaviour of colloidal suspensions. These authors have shown that for a dilute suspension (particle volume fraction $\epsilon _{p} < 0.05$), there is no particle migration when the Péclet numbers are in the range $70 \leq Pe \leq 4400$ ($Pe = 6 {\rm \pi}\mu \dot {\gamma } R_{p}^3/K_{B} T$ characterizes the ratio between the characteristic times of particle diffusion under the effect of Brownian motion and shear in the liquid phase). This lack of particle migration can be attributed to weak particle–particle hydrodynamic interactions in the dilute regime, leading to a negligible effect of $(\epsilon _{l}=1-\epsilon _{p})$ gradients on liquid phase hydrodynamics (Frank et al. 2003), as the volume fraction of the particles, $\epsilon _{p}$, remains uniform for dilute suspensions. Furthermore, the non-slip of the liquid on the surface of each particle causes the perturbation of the liquid velocity field, $\boldsymbol {\hat {v}}_{l}$. However, the flow around the particles can be simplified to a Stokes flow $(Re_{p} << 1)$ for dilute colloidal suspensions $(D_{p} < 1\ \mathrm {\mu }{\rm m})$, dwarfing the contribution of the inertial hydrodynamic dispersion term, $\langle \boldsymbol {\hat {v}}_{l} \otimes \boldsymbol {\hat {v}}_{l} \rangle ^{l}$ in (2.24). With this assumption, the linear momentum conservation equation of the liquid (2.24), reduces to

(3.4)\begin{align} & \rho_{l} \left(\frac{\partial \langle \boldsymbol{v}_{l} \rangle}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{l} \langle \boldsymbol{v}_{l} \rangle^{l} \otimes \langle \boldsymbol{v}_{l} \rangle^{l})\right) =-\boldsymbol{\nabla} \langle P\rangle+ \mu \nabla^2 \langle \boldsymbol{v}_{l}\rangle + \langle \boldsymbol{F}^{e}_{l}\rangle \nonumber\\ &\quad + \frac{1}{V} \int_{A_{{l}-{p}}} ( -\chi_{l}\hat{P} \boldsymbol{\mathsf{I}} + \mu {\boldsymbol{\nabla}} \chi_{l} \boldsymbol{\hat{v}}_{l} )\boldsymbol{\cdot} \boldsymbol{n}_{l} \,{\rm d} A. \end{align}

The liquid average linear momentum equation (3.4), which is similar to the Navier–Stokes equation, requires a closure relation that expresses its last right-hand side term as a function of the average field. This specific term refers to surface forces arising from hydrodynamic interactions between liquid and particles. By adopting the Stokes flow approximation (Guazzelli, Morris & Pic 2011; Kim & Karrila 2013) and assuming the absence of contact forces between particles, particularly in dilute regimes, it is possible to obtain a suitable approximation of the exchange of momentum between the two phases. This approximation involves evaluating the forces exerted by the liquid phase on individual particles, then scaling these forces by the particle volume fraction $\epsilon _{p}$,

(3.5)\begin{equation} \frac{1}{V} \int_{A_{{l}-{p}}} (-\chi_{l}\hat{P} \boldsymbol{\mathsf{I}} + \mu {\boldsymbol{\nabla}} \chi_{l} \boldsymbol{\hat{v}}_{l}) \boldsymbol{\cdot} \boldsymbol{n}_{l} \,{\rm d} A = \langle \boldsymbol{F}^{pl} \rangle = \epsilon_{p} \langle \boldsymbol{F}_{p}^{pl} \rangle, \end{equation}

where $\langle \boldsymbol{F}_{p}^{pl} \rangle$ denotes the average force exerted by a particle on a liquid per unit volume and $\langle \boldsymbol {F}^{pl} \rangle$ refers to the average force per unit volume exerted by the particles in the REV, which can be expressed as

(3.6)\begin{equation} \langle \boldsymbol{F}^{pl} \rangle = \frac{1}{V} \int_{V} \boldsymbol{F}^{pl} \,{\rm d} V. \end{equation}

At the microscopic scale, consider a free sphere that is submerged in a flowing liquid and positioned at $\boldsymbol {y}_0$ within the REV confines. To derive insight into the velocity of the liquid, $\boldsymbol {v}_{l}$, at the microscopic position vector $\boldsymbol {y}$, a Taylor series expansion with respect to the position $\boldsymbol {y}_0$ is performed. This expansion can be represented as follows (Graham 2018):

(3.7)\begin{equation} \boldsymbol{v}_{l}(\kern0.08em \boldsymbol{y})=v_{l}(\kern0.08em \boldsymbol{y}_0) + (\kern0.08em \boldsymbol{y}-\boldsymbol{y}_0) \boldsymbol{\cdot} ({\boldsymbol{\varUpsilon}}+{\boldsymbol{\varLambda}}) + O(|\kern0.08em \boldsymbol{y}-\boldsymbol{y}_0|^2). \end{equation}

Equation (3.7) contains the term $({\boldsymbol{\varUpsilon}}+{\boldsymbol{\varLambda}})$, which represents the tensorial gradient of the velocity field. The symmetric portion, ${\boldsymbol{\varUpsilon}}$, of this tensor characterizes the stresslet, defined as the resistance of the rigid particle to deformation (Batchelor 1970). The stresslet accounts for the additional stresses that arise in the liquid phase due to the rigidity of the suspension particles and exerts a significant influence on the flow behaviour of colloidal suspensions (figure 3a). The hydrodynamic force expressed in (3.6) includes three distinct contributions: (i) the first force results from the stresslet as described by ${\boldsymbol{\varUpsilon}}$ (figure 3a); (ii) the second force relates to the torque force represented by ${\boldsymbol{\varLambda}}$ (figure 3b); and (iii) the last is the drag force (figure 3c). When considering a particle immersed in a flowing liquid with non-zero curvature ($\nabla ^2 \boldsymbol {v}_{l} \neq \text {constant}$), Faxén's laws (Jackson 1997; Guazzelli et al. 2011; Kim & Karrila 2013) provide expressions for stresslet $\boldsymbol{\mathsf{S}}$ (Batchelor 1970), angular moment $\boldsymbol {\varGamma }_{\!\!rot}$ and the drag force $\boldsymbol {F}_{drag}$, per unit volume of a particle,

(3.8)$$\begin{gather} \boldsymbol{\mathsf{S}}= 5 \mu {\boldsymbol{\varUpsilon}}, \end{gather}$$

(3.9)$$\begin{gather}\boldsymbol{\varGamma}_{\!\!rot}= 6 \mu \left(\boldsymbol{w}_{p} - \tfrac{1}{2} \boldsymbol{\nabla} \times \boldsymbol{v}_{l}\right), \end{gather}$$

(3.10)$$\begin{gather}\boldsymbol{F}_{drag} =- \frac{9 \mu}{2 R_{p}^2} \left(\boldsymbol{v}_{l}-\boldsymbol{v}_{p} + \frac{R_{p}^2}{6} \nabla^2 \boldsymbol{v}_{l}\right), \end{gather}$$

where $\boldsymbol{\varGamma}_{\!\!rot}$ is the rotational couple describing the lack of synchrony between liquid vorticity and particle angular velocity.

Figure 3. Schematic description of Faxén's laws (3.8) to (3.10), expressing the momentum exchange between the liquid and particle. (a) The arrows pointing in and out describe the stresslet exerted on the liquid by a particle in the case of a straining flow. (b) The rotating arrows illustrate the moment exerted by a spinning particle on the liquid. (c) The arrow describes the drag force exerted by a flowing liquid on a particle.

Using Faxén's laws (3.8) to (3.10), the particle's force exerted on the liquid can be expressed as follows:

(3.11)\begin{equation} \boldsymbol{F}^{pl} = \boldsymbol{F}_{drag} + \boldsymbol{\nabla} \times \boldsymbol{\varGamma}_{\!\!rot} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\mathsf{S}}. \end{equation}

Substituting (3.11) into (3.6) and then into (3.5), the hydrodynamic force describing the momentum exchange in the two phases can be expressed as follows:

(3.12)\begin{align} & \frac{1}{V} \int_{A_{{l}-{p}}} (-\chi_{l}\hat{P} \boldsymbol{\mathsf{I}} + \mu {\boldsymbol{\nabla}} \chi_{l} \boldsymbol{\hat{v}}_{l}) \boldsymbol{\cdot} \boldsymbol{n}_{l}\,{\rm d} A =- \frac{9 \mu \epsilon_{p}}{2 R_{p}^2} \left(\langle \boldsymbol{v}_{l} \rangle - \langle \boldsymbol{v}_{p} \rangle + \frac{R_{p}^2}{6} \nabla^2 \langle \boldsymbol{v}_{l} \rangle\right) \nonumber\\ &\quad + 6 \mu \epsilon_{p} \left(\boldsymbol{\nabla} \times \langle \boldsymbol{w}_{p} \rangle - \frac{1}{2} \boldsymbol{\nabla} \times \boldsymbol{\nabla} \times \langle \boldsymbol{v}_{l} \rangle\right) + \frac{5}{2}\epsilon_{p} \mu \nabla^2 \langle \boldsymbol{v}_{l} \rangle. \end{align}

The final form of the linear momentum conservation equation for the liquid is obtained by replacing the jump condition in (3.4) with (3.12), where

(3.13)\begin{align} & \rho_{l} \left(\frac{\partial \langle \boldsymbol{v}_{l} \rangle}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{l} \langle \boldsymbol{v}_{l} \rangle^{l} \otimes \langle \boldsymbol{v}_{l} \rangle^{l})\right) =-\boldsymbol{\nabla} \langle P \rangle \nonumber\\ &\quad +\mu \left(1 + \frac{5}{2}\epsilon_{p}\right) \nabla^2 \langle \boldsymbol{v}_{l}\rangle - \frac{9 \mu \epsilon_{p}}{2 R_{p}^2} \left(\langle \boldsymbol{v}_{l} \rangle - \langle \boldsymbol{v}_{p} \rangle + \frac{R_{p}^2}{6} \nabla^2 \langle \boldsymbol{v}_{l} \rangle\right) \nonumber\\ &\quad +6 \mu \epsilon_{p} \left(\boldsymbol{\nabla} \times \langle \boldsymbol{w}_{p} \rangle - \frac{1}{2} \boldsymbol{\nabla} \times \boldsymbol{\nabla} \times \langle \boldsymbol{v}_{l} \rangle\right) + \langle\boldsymbol{F}^{e}_{l}\rangle. \end{align}

3.3. Particle linear momentum jump condition

The closure of the average linear particle momentum (2.25) requires the identification of the forces, $\langle \boldsymbol {F}^{lp} \rangle,\ \langle \boldsymbol {F}^{pp} \rangle, \langle \boldsymbol {F}^{e}_{p} \rangle$, as well as the inertial dispersion term, $\epsilon _{p} \langle \boldsymbol {\hat {v}}_{p} \otimes \boldsymbol {\hat {v}}_{p} \rangle ^{p}$, as a function of the average fields. In the dilute regime, the particles are sufficiently separated from each other, so that the particle–particle contact force $\langle \boldsymbol {F}^{pp} \rangle$ can be ignored. Furthermore, given the rigid nature of the particles, it can be inferred that the stresses occurring within them are zero. For the discrete phase, therefore, the contribution of the stresslet or, in other words, the force $\langle \boldsymbol {F}^{lp}\rangle$ is not relevant. Consequently, this whole set of forces is reduced to $\langle \boldsymbol {F}^{lp}\rangle$ in order to consider only the drag force based on the principle of action and reaction (Newton's third law). Given these assumptions, the force $\langle \boldsymbol {F}^{lp}\rangle$ can be written as follows:

(3.14)\begin{align} \langle \boldsymbol{F}^{lp}\rangle &=- \epsilon_{p} \langle \boldsymbol{F}_{drag}\rangle \nonumber\\ &\Rightarrow\langle \boldsymbol{F}^{lp}\rangle = \frac{9 \mu \epsilon_{p}}{2 R_{p}^2} \left(\langle \boldsymbol{v}_{l} \rangle - \langle \boldsymbol{v}_{p} \rangle + \frac{R_{p}^2}{6} \nabla^2 \langle \boldsymbol{v}_{l} \rangle\right). \end{align}

By the same token, the perturbation of the velocity of the particles, $\boldsymbol {\hat {v}}_{p}$, due to mutual particle–particle interactions is negligible, and thus the term for the inertial particle dispersion, $\epsilon _{p} \langle \boldsymbol {\hat {v}}_{p} \otimes \boldsymbol {\hat {v}}_{p} \rangle ^{p}$, is negligible. Considering all previous assumptions, as well as the expression for the drag force (3.14), the equation for the average linear momentum of the particles (2.25), reduces to

(3.15)\begin{align} \rho_{p} \left(\frac{\partial \langle \boldsymbol{v}_{p} \rangle}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{p} \langle \boldsymbol{v}_{p} \rangle^{p} \otimes \langle \boldsymbol{v}_{p} \rangle^{p})\right) = \frac{9 \mu \epsilon_{p}}{2 R_{p}^2} \left(\langle \boldsymbol{v}_{l} \rangle - \langle \boldsymbol{v}_{p} \rangle + \frac{R_{p}^2}{6} \nabla^2 \langle \boldsymbol{v}_{l} \rangle\right) + \langle \boldsymbol{F}^{e}_{p} \rangle. \end{align}

It is important to note that Brownian particles are neutrally buoyant (Graham 2018). Therefore, the force expression, $\langle \boldsymbol {F}^{e}_{p} \rangle$, can only take into account the external magnetic field effect. This expression will be discussed in detail in § 4, as it is related to the average orientation of the magnetic moment of the particles.

3.4. Particles angular momentum jump condition

Similar to the linear momentum (3.15) for particles, the angular momentum equation also requires defining moments such as $\langle \boldsymbol {\varGamma }^{lp} \rangle,\ \langle \boldsymbol {\varGamma }^{pp} \rangle,\ \langle \boldsymbol {\varGamma }_{p}^{e} \rangle$. However, in the case of a dilute suspension, the rotational moment resulting from particle contact, $\langle \boldsymbol {\varGamma }^{pp} \rangle$, and the inertial dispersion of angular momentum, $\epsilon _{p} \langle \boldsymbol {\hat {v}}_{p} \otimes \boldsymbol {\hat {w}}_{p} \rangle ^{p}$, can be considered negligible (Jackson 1997).

The rotational moment exerted by the liquid on the particles in the dilute regime can be related to that of Faxén's laws (3.9), as follows:

(3.16)\begin{equation} \langle \boldsymbol{\varGamma}^{lp} \rangle =- \epsilon_{p} \langle \boldsymbol{\varGamma}_{\!\!rot} \rangle = 6 \mu \epsilon_{p} \left(\tfrac{1}{2} \boldsymbol{\nabla} \times \langle \boldsymbol{v}_{l} \rangle - \langle \boldsymbol{w}_{p} \rangle\right). \end{equation}

Given the above assumptions and the expression for the moment (3.16), the macroscopic conservation of angular momentum (2.26) is written as follows:

(3.17)\begin{equation} I_{p} \left(\frac{\partial \langle \boldsymbol{w}_{p} \rangle}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\epsilon_{p} \langle \boldsymbol{v}_{p} \rangle^{p} \otimes \langle \boldsymbol{w}_{p} \rangle^{p})\right) = 6 \mu \epsilon_{p} \left(\frac{1}{2} \boldsymbol{\nabla} \times \langle \boldsymbol{v}_{l} \rangle - \langle \boldsymbol{w}_{p} \rangle\right) + \langle \boldsymbol{\varGamma}_{p}^{e} \rangle. \end{equation}

Like the expression of the external force, $\langle \boldsymbol {F}^{e}_{p} \rangle$, the expression of the moment, $\langle \boldsymbol {\varGamma }_{p}^{e} \rangle$, is related to the prevailing magnetic field. The latter will also be discussed in detail in § 4, since it is also related to the average magnetic moment orientation of the particles.

3.5. Maxwell equations jump conditions

In the presence of an external field, a ferromagnetic particle, due to its inherent magnetization, does not have an internal magnetic field identical to the external field. Instead, there is a difference between the external and internal fields known as the demagnetization field. This field is caused by the flux from the magnetization of the particle, which reduces the resulting magnetic field. Specifically, for a particle subjected to a magnetic field propagating in the liquid phase, $\boldsymbol {H}_{l}$, the expression for the magnetic field inside the particle, $\boldsymbol {H}_{p}$, is given by Joseph & Schlömann (1965), Kuznetsov (2018) and Kuznetsov et al. (2022),

(3.18)\begin{equation} \boldsymbol{H}_{p} = \boldsymbol{H}_{l} - \kappa M_{d} \boldsymbol{m}, \end{equation}

where $\kappa$ is the demagnetization factor, which is $1/3$ for a spherical particle (Kuznetsov et al. 2022).

Using the demagnetizing field expression (3.18) as the jump condition, the surface term in the average Ampère–Maxwell equation (2.29) can be written as follows:

(3.19)\begin{equation} \frac{1}{V} \int_{A_{{l}-{p}}} (\chi_{p} \boldsymbol{H}_{p} - \chi_{l} \boldsymbol{H}_{l}) \times \boldsymbol{n}_{p} \,{\rm d} A = \frac{1}{V} \int_{A_{{l}-{p}}} \kappa M_{d} \boldsymbol{n}_{p} \times \chi_{p} \boldsymbol{m} \,{\rm d} A. \end{equation}

Using the following identity (Arfken & Weber 2005):

(3.20)\begin{equation} \frac{1}{V} \int_{A_{{l}-{p}}} \boldsymbol{n}_{p} \times \chi_{p} \boldsymbol{m}\,{\rm d} A = \frac{1}{V} \int_{V} \boldsymbol{\nabla} \times \chi_{p} \boldsymbol{m} \,{\rm d} V \equiv \boldsymbol{\nabla} \times \left\langle \boldsymbol{m} \right \rangle, \end{equation}

the average Ampère–Maxwell (2.29) becomes

(3.21)\begin{equation} \boldsymbol{\nabla} \times (\langle \boldsymbol{H}\rangle + \kappa M_{d} \langle \boldsymbol{m}\rangle) = \boldsymbol{0} .\end{equation}

The superficial term of the average Maxwell-flux equation (2.33) can be modelled using the induction expressions, $\boldsymbol {B}_{l}$ and $\boldsymbol {B}_{p}$, as well as the demagnetizing field (3.18), as follows:

(3.22)\begin{equation} \frac{1}{V} \int_{A_{{l}-{p}}} (\chi_{p} \boldsymbol{B}_{p} - \chi_{l} \boldsymbol{B}_{l}) \boldsymbol{\cdot} \boldsymbol{n}_{p}\,{\rm d} A = (1 - \kappa) M_{d} \frac{1}{V} \int_{A_{{l}-{p}}} \chi_{p} \boldsymbol{m} \boldsymbol{\cdot} \boldsymbol{n}_{p} \,{\rm d} A. \end{equation}

Hence, Gauss’ identity applied to (3.22) gives (Arfken & Weber 2005)

(3.23)\begin{align} & \frac{1}{V} \int_{A_{{l}-{p}}} (\chi_{p} \boldsymbol{B}_{p} - \chi_{l} \boldsymbol{B}_{l}) \boldsymbol{\cdot} \boldsymbol{n}_{p}\,{\rm d} A \nonumber\\ &\quad =(1 - \kappa) M_{d} \frac{1}{V} \int_{V} \boldsymbol{\nabla} \boldsymbol{\cdot} \chi_{p} \boldsymbol{m}\,{\rm d} V \equiv (1 - \kappa) M_{d} \boldsymbol{\nabla} \boldsymbol{\cdot} \langle \boldsymbol{m} \rangle. \end{align}

The final form of the Maxwell-flux equation is obtained by substituting the jump condition (3.23) into (2.33):

(3.24)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (\langle \boldsymbol{H}\rangle + (2 - \kappa) M_{d} \langle \boldsymbol{m}\rangle) = 0. \end{equation}

4. Macroscale orientational dynamics of particles

When exploring the complex dynamics of ferrofluids, it becomes clear that the two-phase approach alone is not sufficient to comprehensively describe the behaviour of these systems subjected to external magnetic fields. A crucial aspect of these dynamics is the average orientation of the magnetic moments of the ferromagnetic particles, which is influenced by magnetic fields as well as by translational and rotational motions resulting from liquid flow or Brownian motion. These phenomena are particularly pronounced in different flow contexts, where regions dominated by inertia show a less sensitive response to the magnetic field, in contrast to regions where liquid shear dominates. To understand this average orientation of magnetic moments, it is essential to separate the characteristic time scales governing flow and Brownian motion. From this perspective, the characteristic time of Brownian motion emerges as the fundamental relaxation time (Graham 2018), imposing a probabilistic approach to accurately represent particle magnetic moments. These considerations lay the foundations for our subsequent investigation in § 4, where the underlying mathematical aspects linking the two-phase model and the Smoluchowski equation are explored. The equations for the concentration and orientation of particle magnetic moments become crucial, not only for understanding the dynamics, but also for the mathematical closure of the two-phase model. It is important to emphasize that the two-phase model and the Smoluchowski equation are complementary, their foundations being distinct, but converging towards a comprehensive description of the behaviour of ferrofluids subjected to external magnetic fields.

Consider a dilute colloidal suspension subjected to an external magnetic field $\boldsymbol {H}_0$. The presence of particles is defined by a probability density, $\varPsi (\boldsymbol {x},\boldsymbol {u},t)$, which depends on the macroscopic position $\boldsymbol {x}$ and an orientation vector $\boldsymbol {u}$ defined on a unit sphere $\varpi$. Integrating $\varPsi$ with respect to all possible orientations on $\varpi$ gives

(4.1)\begin{equation} \int_{\varpi} \varPsi (\boldsymbol{x},\boldsymbol{u},t)\,{\rm d} u = \epsilon_{p}. \end{equation}

The particles’ mass conservation requires that the probability density function, denoted $\varPsi$, satisfies Smoluchowski's equation (Doi & Edwards 1986), formulated as follows:

(4.2)\begin{equation} \frac{\partial \varPsi}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{J}_x \varPsi) + \boldsymbol{\mathcal{L}}_u \boldsymbol{\cdot} (\boldsymbol{J}_u \varPsi)= 0, \end{equation}

where $\boldsymbol {\mathcal {L}}_u =\boldsymbol {u} \times (({\partial }/{\partial \theta }) \boldsymbol {e}_{\theta } + ({1}/{\sin (\theta )}) ({\partial }/{\partial \phi }) \boldsymbol {e}_{\phi })$ is the curl surface gradient operator defined on the unit sphere $\varpi$ (supplementary materials § 3). Here $\boldsymbol {J}_x$ and $\boldsymbol {J}_u$ denote the translational and rotational fluxes, respectively. These fluxes are expressed in the dilute regime as follows:

(4.3)$$\begin{gather} \boldsymbol{J}_x = \langle \boldsymbol{v}_{p}\rangle - D \left(\boldsymbol{\nabla} (\ln (\varPsi)) + \frac{\boldsymbol{\nabla}U}{K_{B} T}\right), \end{gather}$$

(4.4)$$\begin{gather}\boldsymbol{J}_u = \langle \boldsymbol{w}_{p}\rangle - D_u \left(\boldsymbol{\mathcal{L}}_u (\ln (\varPsi)) + \frac{\boldsymbol{\mathcal{L}}_u U}{K_{B} T }\right), \end{gather}$$

where $D=K_{B} T/6 {\rm \pi}\mu R_{p}$ and $D_u=K_{B} T/8 {\rm \pi}\mu R_{p}^3$ denote the translational and rotational diffusion coefficients, respectively; $K_{B}$ and $T$ refer, respectively, to the Boltzmann constant and the absolute temperature; $U$ is the potential energy of magnetic particles; $\alpha = {\mu _0 M_{d} H_0 V_{p}}/{K_{B} T}$ denotes the Langevin parameter.

In our case, ferromagnetic particles are assumed to be single-domain with negligible magnetic anisotropy. Each particle in the suspension has a constant magnetization, but its orientation can be influenced by the magnetic field and thermal energy agitation. The Zeeman potential can be used to model the interactions of the magnetic moments of particles with the external field, where

(4.5)\begin{equation} U_H=- \mu_0 V_{p} \boldsymbol{M} \boldsymbol{\cdot} \langle \boldsymbol{H}\rangle, \end{equation}

where $U_H$ is the Zeeman potential energy of a single particle; $\boldsymbol {M}$ is the particle's magnetization vector, defined as domain magnetization $M_{d}$, which is considered as an intrinsic property of particles, multiplied by the unit vector $\boldsymbol {u}$, which describes the particle's possible orientations with respect to the unit sphere $\varpi$. In this case the expression of the Zeeman potential energy (4.5) is then written as (Jones 2003; Fang 2019; Kuznetsov et al. 2022)

(4.6)\begin{equation} U_H=- K_{B} T \alpha \boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle, \end{equation}

where $\langle \boldsymbol {\breve {H}} \rangle ={\langle \boldsymbol {H}\rangle }/{H_0}$ is the normalized magnetic field by $H_0$, which refers to the magnitude of the external magnetic field $\boldsymbol {H}_0$.

In the case of hydrostatic equilibrium, when the suspension is in a state of quiescence and subjected to a spatially homogeneous, non-rotating magnetic field, the Smoluchowski equation governing the equilibrium probability density, $\varPsi _{eq}$, is formulated under the assumption of the Zeeman potential (4.6):

(4.7)\begin{equation} \nabla_u^2 \varPsi_{eq}- \alpha \boldsymbol{\mathcal{L}}_u \boldsymbol{\cdot} (\varPsi_{eq} \boldsymbol{\mathcal{L}}_u (\boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle)) = 0 .\end{equation}

Equation (4.7), derived from Smoluchowski (4.2) under the assumption of a hydrostatic regime, is identical to that obtained by Doi & Edwards (1986) in the case of a dilute colloidal suspension subjected to an external field of potential $U$ and a purely rotational diffusion motion.

An analytical solution of (4.7) can be obtained (Jones 2003), where

(4.8)\begin{equation} \varPsi_{eq} = \epsilon_{p} \frac{\alpha}{ 4 {\rm \pi}\sinh (\alpha)} \exp (\alpha (\boldsymbol{u}\boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle)). \end{equation}

Equation (4.8) reveals that the equilibrium particle probability density follows the well-known Maxwell–Boltzmann distribution. In this context, the average orientation of particles’ magnetic moments, $\langle \boldsymbol {m}\rangle _{eq}$, can be calculated using the first moment of $\varPsi _{eq}$, such that

(4.9)\begin{equation} \langle \boldsymbol{m}\rangle_{eq} = \int_{\varpi} \boldsymbol{u} \varPsi_{eq}\,{\rm d} u = \epsilon_{p} \mathcal{M} (\alpha)\langle \boldsymbol{\breve{H}}\rangle, \end{equation}

where $\mathcal {M} (\alpha )= \coth (\alpha ) - {1}/{\alpha }$ refers to the Langevin function.

The magnetic moment of a superparamagnetic particle, under the effect of magnetic field, can interact with the moments of neighbouring magnetic particles. These interactions are referred to as DDIs. The strength of the dipolar interactions between ferromagnetic particles depends on several parameters, such as the interparticle distance, which is controlled by the concentration of the suspension, the intensity of the magnetic field and the magnetic moment of the particles, which is determined by the material properties. The DDI potential, $U_{dd}$, of a pair of particles, $i$ and $j$, is given by (Kuznetsov et al. 2022)

(4.10)\begin{equation} U_{dd}= K_{B} T \beta \left(\frac{\boldsymbol{u}_i \boldsymbol{\cdot} \boldsymbol{u}_j}{\vert \boldsymbol{\tilde{r}}_{ij}\vert^3} -3\frac{\boldsymbol{u}_i \boldsymbol{\cdot} \boldsymbol{\tilde{r}}_{ij} \boldsymbol{u}_j \boldsymbol{\cdot} \boldsymbol{\tilde{r}}_{ij}}{\vert \boldsymbol{\tilde{r}}_{ij}\vert^5}\right), \end{equation}

where $\boldsymbol {u}_i$ and $\boldsymbol {u}_j$ denote the particles $i$ and $j$ possible orientations on the unit sphere $\varpi$, respectively, and $\boldsymbol {\tilde {r}}_{ij}={\boldsymbol {r}_{ij}}/{D_{p}}$ refers the vector between the centres of $i$ and $j$ particles, normalized by the particle's diameter $D_{p}$.

The dimensionless number $\beta ={\mu _0(M_{d} V_{p})^2}/{K_{B} T4 {\rm \pi}D_{p}^3}$, in (4.10), is the key parameter controlling the DDI strength, and it represents the ratio of the DDI energy of a pair of particles to the thermal agitation energy. It is clear from (4.10) that the description of DDIs by Smoluchowski's (4.2), by considering the potential $U_{dd}$ (4.10), is impossible using a single-body approach. To account for this phenomenon, a multibody approach is required, the tractability of which is currently impossible to achieve on a continuum scale. However, modelling using approaches such as Langevin dynamics simulations (Berkov, Iskakova & Zubarev 2009; Kuznetsov 2018; Kuznetsov et al. 2022), Brownian dynamics simulations (Soto-Aquino & Rinaldi 2015; Zhao & Rinaldi 2018) or molecular dynamics simulations (Ivanov, Wang & Holm 2004), is possible, but in hydrostatic equilibrium as exemplified in the following studies. For concentrated suspensions ($\beta \geq 3$ and $\alpha >1$), simulations show the formation of chain-like aggregates (Ivanov et al. 2004; Andreu et al. 2012; Faraudo, Andreu & Camacho 2013; Zhao & Rinaldi 2018). These cases are indescribable by the Smoluchowski equation (4.2) in a continuum scale, such as our case. However, if the suspension is sufficiently dilute with low dipole interaction energy and magnetic field, $\beta <1$ and $\alpha <1$, the Weiss mean-field theory allows for the modelling of DDIs on a continuum scale. This theory likewise expands the local magnetic field acting on a particle by an additional term describing the particle's environment (Pshenichnikov, Mekhonoshin & Lebedev 1996; Ivanov & Kuznetsova 2001), where

(4.11)\begin{equation} \boldsymbol{H}^{dd} = \frac{\epsilon_{p} M_{d}}{3} \mathcal{M} (\alpha) \langle \boldsymbol{\breve{H}}\rangle .\end{equation}

Equation (4.11), known as the Weiss mean dipolar field, is valid for $\beta <1$ and a sufficiently diluted suspension $\epsilon _{p} <<1$ (Huke & Lücke 2000; Ivanov et al. 2007). Taking into account the mean dipolar field, the potential $U_{dd}$ is then written as follows:

(4.12)\begin{equation} U_{dd}=- K_{B} T \xi \mathcal{M} (\alpha)\boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle, \end{equation}

where $\xi =8 \epsilon _{p} \beta$ is the initial magnetic susceptibility.

The total potential of the suspension is a linear combination of the Zeeman potential (4.6), and the mean dipolar field potential (4.12), such that

(4.13)\begin{equation} U=- K_{B} T \mathcal{F}(\alpha)\boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle, \end{equation}

where $\mathcal {F} (\alpha )=\alpha +\xi \mathcal {M}(\alpha )$.

The Smoluchowski (4.2) with the potential $U$ (4.13) is then written as

(4.14)\begin{align} & \frac{\partial \varPsi}{\partial t} + \langle \boldsymbol{v}_{p} \rangle \boldsymbol{\cdot} \boldsymbol{\nabla} \varPsi + \langle \boldsymbol{w}_{p} \rangle \boldsymbol{\cdot} \boldsymbol{\mathcal{L}}_u \varPsi= D [\nabla^2 \varPsi - \mathcal{F}(\alpha) \boldsymbol{\nabla} \boldsymbol{\cdot} (\varPsi \boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}} \rangle))] \nonumber\\ &\quad +D_u [\nabla_u^2 \varPsi - \mathcal{F}(\alpha) \boldsymbol{\mathcal{L}}_u \boldsymbol{\cdot} (\varPsi \boldsymbol{\mathcal{L}}_u(\boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle))] . \end{align}

It is impossible to obtain an analytical solution for the Smoluchowski's (4.14) due to its phase-space dependence on position $\boldsymbol {x}$ and orientation $\boldsymbol {u}$. Consequently, a numerical approach is required to solve this equation. Nevertheless, analytical progress has been made by developing the probability density as a series of its moments (Saintillan & Shelley 2008a,b; Ezhilan & Saintillan 2015; Saintillan & Shelley 2015; Theillard & Saintillan 2019). To define the probability distribution of particles as a function of these moments, we introduce a set of tensors $\mathcal {T}^k_{i_1,i_2,\ldots,i_k}$ equivalent to spherical harmonics. The zeroth moment represents the scalar field of particle concentration $\langle \epsilon _{p}(\boldsymbol {x},t)\rangle$. The first-order moment of the probability density, $\varPsi$, is a vector field describing the macroscopic orientation of the magnetic moments of the particles $\langle \boldsymbol {m}(\boldsymbol {x},t)\rangle$. The second-order moment describes the nematic order characterized by the tensor $\langle \boldsymbol {Q} (\boldsymbol {x},t)\rangle$. The tensors, $\mathcal {T}^k_{i_1,i_2,\ldots,i_k}$, are homogeneous polynomials of degree $k$ expressed in terms of the unit vector $\boldsymbol {u}$ and are fully symmetric in terms of the indices $i_1,i_2,\ldots,i_k$. Each tensor $\mathcal {T}^k$ is orthogonal to the others and is normalized by a scalar product defined as

(4.15)\begin{equation} (a,b) = \frac{1}{\varpi} \int_{\varpi} a(\boldsymbol{u}) b(\boldsymbol{u})\,{\rm d} u .\end{equation}

The first three tensors $\mathcal {T}^k$ are given by

(4.16) \begin{equation} \left.\begin{array}{l@{}} \mathcal{T}^0=1 ,\\ \mathcal{T}^1_i= u_i ,\\ \mathcal{T}^1_{i,j}= u_i u_j - \dfrac{\delta_{ij}}{3}. \end{array}\right\} \end{equation}

The particles probability density $\varPsi (\boldsymbol {x},\boldsymbol {u},t)$ admits an exact expansion on the basis of these tensors (Ahmadi, Marchetti & Liverpool 2006), where

(4.17)\begin{equation} \varPsi(\boldsymbol{x},\boldsymbol{u},t) = \sum_{k=0}^{k} a^m_{i1,\ldots,i_k}(\boldsymbol{x},t)\boldsymbol{\cdot} \mathcal{T}^k_{i1,\ldots,i_k}(\boldsymbol{u}), \end{equation}

where $a^m_{i1,\ldots,i_k} (\boldsymbol {x},t )$ denotes $k$th moment of $\varPsi$. The $k$th moment is a tensor that can be calculated using the definition of the scalar product of orthogonal bases (4.15). The application of this scalar product to the expansion of $\varPsi$ (4.17) yields

(4.18)\begin{equation} a^m_{j1,\ldots,j_k}(\boldsymbol{x},t)\frac{1}{\varpi} \int_{\varpi} \mathcal{T}^k_{j1,\ldots,j_k} (\boldsymbol{u}) \mathcal{T}^k_{i1,\ldots,i_k} (\boldsymbol{u})\,{\rm d} u = \frac{1}{\varpi} \int_{\varpi} \mathcal{T}^k_{i1,\ldots,i_k} \varPsi(\boldsymbol{x},\boldsymbol{u},t)\,{\rm d} u .\end{equation}

The three first moments of $\varPsi$ are given by

(4.19)$$\begin{gather} a^0 = \frac{1}{\varpi} \int_{\varpi} \varPsi(\boldsymbol{x},\boldsymbol{u},t)\,{\rm d} u = \frac{\epsilon_p(\boldsymbol{x},t)}{\varpi}, \end{gather}$$

(4.20)$$\begin{gather}\boldsymbol{a}^1 = \frac{3}{\varpi} \int_{\varpi} \boldsymbol{u} \varPsi(\boldsymbol{x},\boldsymbol{u},t)\,{\rm d} u =\frac{\langle \boldsymbol{m}(\boldsymbol{x},t)\rangle}{\varpi}, \end{gather}$$

(4.21)$$\begin{gather}\boldsymbol{\mathsf{a}}^2 = \frac{d(d+2)}{2\varpi} \int_{\varpi} \left(\boldsymbol{u} \otimes \boldsymbol{u} - \frac{\boldsymbol{\mathsf{I}}}{3}\right) \varPsi(\boldsymbol{x},\boldsymbol{u},t)\,{\rm d} u =\frac{\langle \boldsymbol{\mathsf{Q}}(\boldsymbol{x},t)\rangle}{\varpi}, \end{gather}$$

where $d$ denotes the space dimension.

In particular, we focus on the zeroth and first moments, which describe the particle concentration and the macroscopic orientation of the particle magnetic moments. The probability density truncated to first order gives

(4.22)\begin{equation} \varPsi(\boldsymbol{x},\boldsymbol{u},t) = \frac{1}{4 {\rm \pi}}\langle \epsilon_{p}(\boldsymbol{x},t)\rangle + \frac{3}{4 {\rm \pi}} \boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{m} (\boldsymbol{x},t)\rangle .\end{equation}

The zeroth and first moments of the Smoluchowski (4.14), when combined with the probability density expansion expression (4.22), yield a system of two coupled equations governing the macroscopic dynamics of particle transport and their magnetic moment orientation. The supplementary material (§ 4) contains the details of the calculations involving the use of angular operators in the case where the probability density is approximated by (4.22). The zeroth moment of (4.14) gives

(4.23)\begin{equation} \frac{\partial \langle \epsilon_{p}\rangle }{\partial t}+ \langle\boldsymbol{v}_{p} \rangle \boldsymbol{\cdot} \boldsymbol{\nabla} \langle \epsilon_{p}\rangle = D (\nabla^2 \langle \epsilon_{p}\rangle - \mathcal{F}(\alpha) \boldsymbol{\nabla}\boldsymbol{\cdot} (\langle \boldsymbol{m}\rangle \boldsymbol{\cdot} {\boldsymbol{\nabla}} \langle \boldsymbol{\breve{H}}\rangle)) \end{equation}

and that of the first order yields

(4.24)\begin{align} & \frac{\partial \langle \boldsymbol{m}\rangle}{\partial t}+ \langle\boldsymbol{v}_{p}\rangle \boldsymbol{\cdot} {\boldsymbol{\nabla}} \langle \boldsymbol{m}\rangle = D \left(\nabla^2 \langle\boldsymbol{m}\rangle - \frac{\mathcal{F}(\alpha)}{3} \boldsymbol{\nabla} \boldsymbol{\cdot} (\langle \epsilon_{p}\rangle {\boldsymbol{\nabla}} \langle \boldsymbol{\breve{H}}\rangle)\right) \nonumber\\ &\quad +\langle \boldsymbol{w}_{p}\rangle \times \langle \boldsymbol{m}\rangle + \frac{1}{\tau_{B}} (\langle \boldsymbol{m}_0 \rangle - \langle \boldsymbol{m}\rangle), \end{align}

where $\tau _{B} = {1}/{2D_u}$ denotes the Debye relaxation time, and $\langle \boldsymbol {m}_0\rangle$ refers to the average equilibrium orientation which reads as

(4.25)\begin{equation} \langle \boldsymbol{m}_0\rangle = \epsilon_{p} \frac{\mathcal{F}(\alpha)}{3} \langle \boldsymbol{\breve{H}}\rangle. \end{equation}

It is important to note that the equations resulting from the expansion of the particle probability density, in particular (4.23) and (4.24), play an essential role as closure relations for the two-phase model. These equations are not formulated by an averaging over the control volume (REV), but result from the expansion of the probability density. Their presence is crucial for the mathematical closure of the established model. It should also be noted that although the REV averaging and Smoluchowski approaches are distinct, they are complementary and cannot be integrated interchangeably in our context. Equation (4.23) describes particle transport under the effect of magnetophoresis, which is denoted by the term $\boldsymbol {\nabla } \boldsymbol {\cdot } (\langle \boldsymbol {m}\rangle \boldsymbol {\cdot } {\boldsymbol{\nabla}}\langle \boldsymbol {\breve {H}}\rangle )$. Equation (4.24) models the transport of the average orientation of the particles under the effects of the magnetic field, the linear and angular velocity of the particles, and Brownian relaxation, which is characterized by the Debye relaxation time.

Figure 4 depicts the evolution of the average orientation of the particles in hydrostatic equilibrium, $\vert \langle \boldsymbol {m}_{eq}\rangle \vert / \epsilon _{p}$, as a function of $\alpha$. In terms of the effect of DDI on $\vert \langle \boldsymbol {m}_{eq}\rangle \vert / \epsilon _{p}$, the figure shows that for $\xi =0.1$, the predictions of (4.9) and (4.27) are almost identical. In the cases of $\xi =1$ and $2$, there was a significant increase in the orientation at equilibrium when compared with the case without the DDI. In the absence of DDI, $\xi \rightarrow 0$, $\vert \langle \boldsymbol {m}_{eq}\rangle \vert / \epsilon _{p}$ tends to ${\alpha }/{3}$, corresponding to a Taylor series development to order $1$ of the Langevin function (4.9). According to the comparison of ${\alpha }/{3}$ and $\mathcal {M}(\alpha )$ in figure 4, the relaxation term, ${1}/{\tau _{B}} (\langle \boldsymbol {m}_0 \rangle - \langle \boldsymbol {m}\rangle )$, is valid only for a weak magnetic field, where $\alpha < 1$. Indeed, this is due to the truncation of the probability density $\varPsi$ at the first moment. The identification of the equilibrium probability density $\varPsi _{eq}$ (4.8), against the truncated spherical harmonic expansion to first order (4.22), clearly demonstrates that the average equilibrium orientation of the particles is limited to the case $\alpha < 1$, where

(4.26)\begin{equation} \varPsi_{eq} = \frac{\epsilon_{p}}{4 {\rm \pi}} +\frac{3}{4 {\rm \pi}} \frac{\alpha}{3} (\boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle). \end{equation}

Figure 4. The evolution of particle orientation in hydrostatic equilibrium, $\vert \langle \boldsymbol {m}_{eq}\rangle \vert / \epsilon _{p}$, as a function of the Langevin parameter $\alpha$. The continuous curve was calculated using (4.9), while the dashed curves were calculated using (4.25).

Restricting to the case when $\alpha <1$, the first-order Taylor series expansion of the average equilibrium orientation of the particles (4.25) reduces it to

(4.27)\begin{equation} \langle \boldsymbol{m}_0 \rangle = \epsilon_{p} \frac{\mathcal{F}_1(\alpha)}{3} \langle \boldsymbol{\breve{H}}\rangle = \epsilon_{p} \frac{\alpha}{3}\left(1 + \frac{\xi}{3}\right) \langle \boldsymbol{\breve{H}}\rangle. \end{equation}

In conclusion, the average particle orientation (4.24) is expected to be valid when $\alpha <1$. In this case, the effect of magnetophoresis is negligible, resulting in a uniform particle concentration field, which a posteriori justifies our depiction of the two-phase ferrofluid suspension as dilute and clusterless. It is also worth noting that the equilibrium orientation (4.27), is calculated from the probability density in hydrodynamic equilibrium and for the case of a static magnetic field. Our case, on the other hand, corresponds to a rotating magnetic field which is time dependent. Thus, the relaxation term, which involves $\langle \boldsymbol {m}_0 \rangle$, ignores the influence of the magnetic field's time dependence as well as the hydrodynamics of the liquid phase. Therefore, the magnetic field rotation's characteristic time, $\tau _H$, must be greater than that of the Brownian relaxation, $\tau _{B}$. In the case $\tau _H < \tau _{B}$, the orientation of the particles fluctuates significantly to the point where the continuum approach is no longer valid. Consequently, the present theory is necessarily confined solely to $\tau _H > \tau _{B}$ (Fang 2019).

The final point required for the mathematical closure of the ferrofluid two-phase model is the definition of the force and magnetic torque acting on the particles in terms of the average orientation dynamics of the particles’ magnetic moments. The force and torque per unit volume, denoted by $\langle \boldsymbol {F}^{e}_u \rangle$ and $\langle \boldsymbol {\varGamma }_u^{e} \rangle$, acting on a particle whose position is described by the probability density, $\varPsi$, are determined as follows:

(4.28)$$\begin{gather} \langle \boldsymbol{F}^{e}_u \rangle =- \frac{1}{V_{p}} \boldsymbol{\nabla} (U) \varPsi = \frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) \varPsi \boldsymbol{\nabla} (\boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle) , \end{gather}$$

(4.29)$$\begin{gather}\langle \boldsymbol{\varGamma}_u^{e} \rangle =- \frac{1}{V_{p}} \boldsymbol{\mathcal{L}}_u (U) \varPsi = \frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) \varPsi \boldsymbol{\mathcal{L}}_u (\boldsymbol{u} \boldsymbol{\cdot} \langle \boldsymbol{\breve{H}}\rangle) . \end{gather}$$

It is important to note that $\langle \boldsymbol {F}^{e}_u \rangle$ and $\langle \boldsymbol {\varGamma }_u^{e} \rangle$ are derived from the total magnetic potential (4.13), and are multiplied by the particles’ probability density to incorporate the influence of Brownian motion as described by the Smoluchowski equation. This force expresses the interaction between the magnetic field, which attempts to create order, and the disorder caused by Brownian motion.

The magnetic force and torque expressions, $\langle \boldsymbol {F}^{e}_{p} \rangle$ and $\langle \boldsymbol {\varGamma }_{p}^{e} \rangle$, required for the closure of the linear and angular momentum conservation equations of the particles (3.15) and (3.17), are obtained by averaging the expressions (4.28) and (4.29) over the possible orientations $\boldsymbol {u}$ in the unit sphere $\varpi$,

(4.30)$$\begin{gather} \langle \boldsymbol{F}^{e}_{p} \rangle = \frac{1}{\varpi}\int_{\varpi} \langle \boldsymbol{F}^{e}_u \rangle\,{\rm d} u = \frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) \langle \boldsymbol{m} \rangle \boldsymbol{\cdot} {\boldsymbol{\nabla}} \langle \boldsymbol{\breve{H}}\rangle, \end{gather}$$

(4.31)$$\begin{gather}\langle \boldsymbol{\varGamma}_p^{e} \rangle = \frac{1}{\varpi} \int_{\varpi} \langle \boldsymbol{\varGamma}_u^{e} \rangle \,{\rm d} u = \frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) \langle \boldsymbol{m} \rangle \times \langle \boldsymbol{\breve{H}} \rangle, \end{gather}$$

where $du$ denotes the infinitesimal surface element of the unit sphere, $\varpi$.

5. Spin-up flow equations

The derived two-phase model is used for the spin-up flow geometry discussed in § 2. It is important to emphasize that the model includes continuity equations (3.2) and (3.3), linear momentum conservation equations (3.13) and (3.15), particle angular momentum conservation (3.17) and equations governing particle transport and magnetic moment orientation (4.23) and (4.24). In figure 1(a), the Cartesian representation of the external magnetic field $\boldsymbol {H}_0$ is shown as follows:

(5.1)\begin{equation} \boldsymbol{H}_0= H_0 \begin{pmatrix} \breve{H}_{x,0} \\ \breve{H}_{y,0} \\ \breve{H}_{z,0} \end{pmatrix}=H_0 \begin{pmatrix} - \cos (w_0 t) \\ - \sin (w_0 t) \\ 0 \end{pmatrix}, \end{equation}

where its equivalent in cylindrical coordinates is written as

(5.2)\begin{equation} \boldsymbol{H}_0 = H_0 \begin{pmatrix} \breve{H}_{r,0} \\ \breve{H}_{\theta,0} \\ \breve{H}_{z,0} \end{pmatrix}=H_0 \begin{pmatrix} - \cos (w_0 t - \theta) \\ \sin (w_0 t - \theta) \\ 0 \end{pmatrix}. \end{equation}

To make the mathematical formulation of the problem easier, all the mean fields denoted by $\langle \psi \rangle$ will be noted $\psi$. The magnetic field's characteristic rotation time $\tau _H= 1/f_0$ is assumed to be greater than the Brownian relaxation time $\tau _{B}$. In this case, the inertia caused by the rotation of the magnetic field can be ignored in the linear liquid and particles momentum balances (3.13) and (3.15). The volume force in the liquid phase can be ignored because it is entirely due to gravity. The contribution of the liquid shear term $\nabla ^2 \boldsymbol {v}_{l}$ in the drag force of the particles (3.15), which is of the order of $O(R_{p}^2)$, is small in comparison with that of the rotational moment of the particles, which is of the order of $O(1)\gt \gt O(R_{p}^2)$. In this case, we expect that the particle translation contribution is insignificant in comparison with the rotational motion. Given the above assumptions, the linear momentum conservation equations for the liquid and particles (3.13) and (3.15), are written as

(5.3)\begin{gather} \mu \left(1 + \frac{5}{2}\epsilon_{p}\right) \nabla^2 \boldsymbol{v}_{l}- \frac{9 \mu \epsilon_{p}}{2R_{p}^2}\left(\boldsymbol{v}_{l} - \boldsymbol{v}_{p} + \frac{R_{p}^2}{6} \nabla^2 \boldsymbol{v}_{l}\right) \nonumber\\ \quad +\,6 \mu \epsilon_{p} \left(\boldsymbol{\nabla} \times \boldsymbol{w}_{p} - \frac{1}{2} \boldsymbol{\nabla} \times \boldsymbol{\nabla} \times \boldsymbol{v}_{l}\right) - \boldsymbol{\nabla} P = 0, \end{gather}

(5.4)\begin{gather} \frac{9 \mu \epsilon_{p}}{2 R_{p}^2}\left(\boldsymbol{v}_{l} - \boldsymbol{v}_{p} + \frac{R_{p}^2}{6} \nabla^2 \boldsymbol{v}_{l}\right) + \frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) \boldsymbol{m} \boldsymbol{\cdot} {\boldsymbol{\nabla}} \boldsymbol{\breve{H}} = 0 .\end{gather}

The particle linear momentum conservation equation (5.4) demonstrates that the drag term is balanced by the external force applied to the particles by the magnetic field. This term also exists in the linear momentum of the liquid (5.3), but with an opposite sign. Because we are interested in the flow of the liquid phase generated by particle rotation, it would be appropriate to substitute the expression of the drag term of (5.4) in (5.3) to eliminate the velocity of the particles and obtain an equation governing only the flow of the liquid phase. Following this modification, (5.3) is written as

(5.5)\begin{align} & \mu \left(1 + \frac{5}{2}\epsilon_{p}\right) \nabla^2 \boldsymbol{v}_{l} + 6 \mu \epsilon_{p} \left(\boldsymbol{\nabla} \times \boldsymbol{w}_{p} - \frac{1}{2} \boldsymbol{\nabla} \times \boldsymbol{\nabla} \times \boldsymbol{v}_{l}\right) \nonumber\\ &\quad +\frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) \boldsymbol{m} \boldsymbol{\cdot} {\boldsymbol{\nabla}} \boldsymbol{\breve{H}} - \boldsymbol{\nabla} P = 0 . \end{align}

In the absence of flow along the $z$ direction, the liquid velocity field is directed only along the azimuthal $\theta$ direction, where it is easy to show through the continuity equation of the liquid (3.2),

(5.6)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}_{l} = 0 \equiv \frac{1}{r} \frac{\partial v_{l,\theta}}{\partial \theta} = 0, \end{equation}

that $v_{l,\theta }$ depends only on the radial coordinate $r$, such that $v_{l,\theta }=v(r)$. In this configuration, the angular velocity of the particles is directed to the $z$-axis, resulting in

(5.7)\begin{equation} \boldsymbol{w}_{p}=\begin{pmatrix} w_{{p},r} \\ w_{{p},\theta} \\ w_{{p},z} \end{pmatrix}=\begin{pmatrix} 0 \\ 0 \\ w(r) \end{pmatrix}. \end{equation}

Taking into account the radial coordinate dependency of the liquid velocity $v$ as well as the expression for the angular velocity of the particles (5.7), the equation for the linear momentum of the liquid (5.5) reduces to

(5.8)\begin{equation} \mu_e \left(\frac{{\rm d}^2v}{{\rm d} r^2} + \frac{1}{r} \frac{{\rm d} v}{{\rm d} r} - \frac{v}{r^2}\right) - 6 \mu \epsilon_{p} \frac{{\rm d}}{{\rm d} r} (w - \varOmega) = 0 \quad 0 < r < R, \end{equation}

where $\mu _e = \mu (1 + \frac {5}{2} \epsilon _{p})$ denotes the Einstein viscosity. Here $\varOmega = \frac {1}{2} ({1}/{r}) ({{\rm d} (rv)}/{{\rm d} r})$ refers to the vorticity of the liquid along the $z$ direction.

Equation (5.8) is subject to two boundary conditions. The first, at $r=0$, is due to the pipe's cylindrical symmetry, such that

(5.9)\begin{equation} v|_{r=0} = 0 \end{equation}

the second, at $r=R$, is due to the liquid's adherence to the wall, where

(5.10)\begin{equation} v|_{r=R} = 0 .\end{equation}

Under the same assumptions, the conservation of particles angular momentum equation (3.17) reduces to

(5.11)\begin{equation} 6 \mu \epsilon_{p} (\varOmega - w) + \frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) (\breve{H}_{\theta} m_r - \breve{H}_{r} m_{\theta}) = 0 \quad 0 \leq r \leq R .\end{equation}

Taking into account the pipe's radial symmetry, the equation describing particle transport (4.23) is written as follows:

(5.12)\begin{equation} \frac{\partial \epsilon_{p}}{\partial t} = D \left(\frac{\partial^2 \epsilon_{p}}{\partial r^2} + \frac{1}{r} \frac{\partial \epsilon_{p}}{\partial r} - \mathcal{F}_1 (\alpha) \left(m_r \frac{\partial^2 \breve{H}_r}{\partial r^2} + \frac{1}{r} \frac{\partial \breve{H}_r}{\partial r} \frac{\partial (r m_r)}{\partial r}\right)\right)\quad 0< r< R . \end{equation}

This equation has two boundary conditions as well as an initial condition. The first boundary condition is due to the pipe's cylindrical symmetry, which assumes that

(5.13)\begin{equation} \left.\frac{\partial \epsilon_{p}}{\partial r}\right|_{r=0} = 0 .\end{equation}

The second condition, at $r=R$, can be obtained by considering the wall's impermeability to particles. This can be obtained by assuming that the translational flux (4.3) is zero in the normal direction of the wall,

(5.14)\begin{equation} (\boldsymbol{J}_x \varPsi) \boldsymbol{\cdot} \boldsymbol{n} = 0, \end{equation}

where $\boldsymbol {n}$ is the unit vector normal to the wall.

The zeroth-order moment of (5.14) gives the particle transport condition at $r=R$, where

(5.15)\begin{equation} \int_{\varpi} (\boldsymbol{J}_x \varPsi \boldsymbol{\cdot} \boldsymbol{n})\,{\rm d} u = 0 \Rightarrow \left. \frac{\partial \epsilon_{p}}{\partial r}\right|_{r=R} = \mathcal{F}_1 (\alpha) \left.\left(m_r \frac{\partial \breve{H}_r}{\partial r}\right)\right|_{r=R} . \end{equation}

The initial condition is obtained by considering the uniformity of the particle distribution:

(5.16)\begin{equation} \epsilon_{p}(r,t=0) = \epsilon_{{p},0}. \end{equation}

Because the magnetic field rotates around the $z$ axis, the equations governing the dynamics of particles orientation only consider the $r$ and $\theta$ directions, resulting in

(5.17)\begin{align} \frac{\partial m_r}{\partial t} &= D \left[\left(\frac{\partial }{\partial r} \left(\frac{1}{r} \frac{\partial (r m_{r})}{\partial r}\right)\right) - \frac{\mathcal{F}_1 (\alpha)}{3} \left(\epsilon_p \frac{\partial}{\partial r} \left(\frac{1}{r} \frac{\partial (r \breve{H}_{r}) }{\partial r}\right) + \frac{\partial \epsilon_p}{\partial r} \frac{\partial \breve{H}_r}{\partial r}\right)\right] \nonumber\\ &\quad + \frac{1}{\tau_{B}} \left(\epsilon_{p} \frac{\mathcal{F}_1 (\alpha)}{3} \breve{H}_r - m_r\right)- w m_{\theta}\qquad 0< r< R , \end{align}

(5.18)\begin{align} \frac{\partial m_{\theta}}{\partial t} &= D \left[\left(\frac{\partial}{\partial r} \left(\frac{1}{r} \frac{\partial (r m_{\theta})}{\partial r}\right)\right) - \frac{\mathcal{F}_1 (\alpha)}{3} \left(\epsilon_p \frac{\partial }{\partial r} \left(\frac{1}{r} \frac{\partial (r \breve{H}_{\theta}) }{\partial r}\right)+\frac{\partial \epsilon_p}{\partial r} \frac{\partial \breve{H}_{\theta}}{\partial r}\right)\right] \nonumber\\ &\quad +\frac{1}{\tau_{B}} \left(\epsilon_{p} \frac{\mathcal{F}_1 (\alpha)}{3} \breve{H}_{\theta} -m_{\theta}\right) + w m_r \qquad 0< r< R. \end{align}

Each of (5.17) and (5.18), which characterize the components of the orientation vector $m_r$ and $m_{\theta }$, requires two boundary conditions and an initial condition. Two first boundary conditions can be obtained by considering the cylindrical symmetry of the pipe, at $r=0$:

(5.19)$$\begin{gather} m_r|_{r=0} = 0, \end{gather}$$

(5.20)$$\begin{gather}m_{\theta}|_{r=0} = 0. \end{gather}$$

Consider the first moment of the particle wall impermeability equation (5.14) at $r=R$ to obtain the two remaining boundary conditions:

(5.21)$$\begin{gather} \left.\frac{\partial m_r}{\partial r}\right|_{r=R} = \frac{\mathcal{F}_1 (\alpha)}{3} \epsilon_{p} \left.\frac{\partial \breve{H}_r}{\partial r}\right|_{r=R}, \end{gather}$$

(5.22)$$\begin{gather}\left.\frac{\partial m_{\theta}}{\partial r}\right|_{r=R} = \frac{\mathcal{F}_1 (\alpha)}{3} \epsilon_{p} \left.\frac{\partial \breve{H}_{\theta}}{\partial r}\right|_{r=R}. \end{gather}$$

The two initial conditions on $m_r$ and $m_{\theta }$ can be obtained by considering an absence of the magnetic field at $t=0$:

(5.23)\begin{equation} m_r(r,t=0)=m_{\theta}(r,t=0) = 0. \end{equation}

All of the ferrofluid transport equations are coupled to the induced magnetic field within the suspension (3.21) and (3.24). By taking into account the pipe's symmetry and the direction of the external field $\boldsymbol {H}_0$, these equations reduce to

(5.24)$$\begin{gather} \frac{\partial \breve{H}_{\theta}}{\partial r} + \frac{\breve{H}_{\theta}}{r} +\frac{1}{3} \frac{M_{d}}{H_0} \left(\frac{\partial m_{\theta}}{\partial r} + \frac{m_{\theta}}{r}\right) = 0 \qquad 0< r< R , \end{gather}$$

(5.25)$$\begin{gather}\frac{\partial \breve{H}_{r}}{\partial r} + \frac{\breve{H}_{r}}{r} + \frac{5}{3} \frac{M_{d}}{H_0} \left(\frac{\partial m_{r}}{\partial r} + \frac{m_{r}}{r}\right) = 0 \qquad 0< r< R. \end{gather}$$

Each of (5.24) and (5.25) requires the definition of a boundary condition at the wall, $r=R$. The integral form of the Maxwell–Ampère equation (3.21)) reads

(5.26)\begin{equation} \int_{A} \left(\boldsymbol{H} + \frac{M_{d}}{3} \boldsymbol{m}\right) \times \boldsymbol{n}\,{\rm d} A = 0. \end{equation}

When the suspension is assumed to be confined by a non-magnetic and non-conductive wall, the integral form of the Maxwell–Ampère equation reveals the necessity of a continuity condition linking the tangential component of $\boldsymbol {H} + \frac {5}{3} M_{d} \boldsymbol {m}$ and that of $\boldsymbol {H}_0$, where

(5.27)\begin{equation} \left(\boldsymbol{H} + \frac{M_{d}}{3} \boldsymbol{m}\right) \times \boldsymbol{n} = \boldsymbol{H}_0 \times \boldsymbol{n}. \end{equation}

For the case of (5.24), the condition (5.27) becomes

(5.28)\begin{equation} \breve{H}_{\theta}|_{r=R} + \frac{1}{3} \frac{M_{d}}{H_0} m_{\theta}|_{r=R} = \breve{H}_{\theta,0}. \end{equation}

Similarly to the Ampère–Maxwell equation, the Maxwell-flux (3.24) in its integral form can be written as follows:

(5.29)\begin{equation} \int_{A} \left(\boldsymbol{H} + \frac{5}{3} M_{d}\boldsymbol{m}\right) \boldsymbol{\cdot} \boldsymbol{n}\,{\rm d} A = 0, \end{equation}

where it is apparent through this representation that a continuity is necessary between the wall normal component of $\boldsymbol {H} + \frac {5}{3} M_{d} \boldsymbol {m}$ with that of $\boldsymbol {H}_0$, such that

(5.30)\begin{equation} \left(\boldsymbol{H} + \tfrac{5}{3} M_{d}\boldsymbol{m}\right) \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{H}_0 \boldsymbol{\cdot} \boldsymbol{n}. \end{equation}

In the case of (5.25), the continuity of the normal component reduces to

(5.31)\begin{equation} \breve{H}_r|_{r=R} + \frac{5}{3} \frac{M_{d}}{H_0} m_r|_{r=R} = \breve{H}_{r,0}. \end{equation}

The solution of the ferrofluid model for the case of a spin up flow geometry can be obtained by numerically solving (5.8), (5.11), (5.12), (5.17), (5.18), (5.24) and (5.25) considering the boundary conditions (5.9), (5.10), (5.13), (5.15), (5.19) to (5.22), (5.28) and (5.31). The stationary solution can be obtained by writing the particles’ concentration field $\epsilon _{p}$, the orientation of the particles’ magnetic moments $\boldsymbol {m}$, and the magnetic field $\boldsymbol {\breve {H}}$ in the following functional form:

(5.32)$$\begin{gather} \epsilon_{p} = Re \lbrace \underline{\epsilon_{p}} \exp (\,j w_0 t)\rbrace, \end{gather}$$

(5.33)$$\begin{gather}\boldsymbol{m} = Re \lbrace \boldsymbol{\underline{m}} \exp (\,j w_0 t)\rbrace, \end{gather}$$

(5.34)$$\begin{gather}\boldsymbol{\breve{H}} = Re \lbrace \boldsymbol{\breve{\underline{H}}} \exp (\,j w_0 t)\rbrace, \end{gather}$$

where $\underline {\epsilon _{p}}$, $\boldsymbol {\underline {m}}$ and $\boldsymbol {\underline {\breve {H}}}$ are the complex amplitudes of $\epsilon _{p}$, $\boldsymbol {m}$ and $\boldsymbol {\breve {H}}$.

6. Results and discussion

This section presents the results of the two-phase ferrofluid model for a cylindrical spin-up flow configuration. To confirm the suitability of the model within this specific framework, the relevance of the experimental conditions with respect to the underlying model assumptions is first examined. Subsequently, the relevant experimental results are confronted with the model in order to allow a well-founded comparison. Finally, a parametric analysis is performed to evaluate the different contributions of the phenomena emphasized by the mathematical model.

6.1. Model confrontation with experimental results for $\alpha <1$ and $w_0 \tau _{B}<<1$

Torres-Diaz et al. (2014) performed experimental measurements of the local azimuthal velocity field using an ultrasonic technique in a cylindrical spin-up flow geometry (tube with radius $R=24.7\ {\rm mm}$ and length $L=63.5\ {\rm mm}$) induced by a uniform, rotating external magnetic field. These experimental results were selected for comparison with the theoretical results because the ferrofluids used have a concentration of less than $1\,\%$, which is consistent with the dilution limit assumption used in model development. Two categories of ferrofluids were examined in the experimental study. The first category included water-based ferrofluids (WBFs) and was further subdivided into three types (WBF-1, WBF-2, WBF-3), which were distinguished by mean particle diameter $(D_{p}=14.3,\ 17.2,\ 28.9\ {\rm nm})$ and particle concentration $(\epsilon _{p}=0.00213, 0.00361,\ 0.00756)$. The second category comprised oil-based ferrofluids, including WBGF-1 $(D_{p}=14.5\ {\rm nm},\ \epsilon _{p}=0.00376)$, GBF-1 $(D_{p}=12. 3\ {\rm nm}, \epsilon _{p}=0.00183)$ and GBF2 $(D_{p}=14.4\ {\rm nm},\ \epsilon _{p}=0.00753)$. The measured domain magnetization for both categories of ferrofluids prepared with cobalt ferrite particles was $M_{d}=425\ {\rm kA}\ {\rm m}^{-1}$.

As explained in § 4, despite the dilute nature of the aforementioned ferrofluids, the equations governing the transport of the particles and the dynamics of their magnetic moment orientation should be valid only under certain conditions. In particular, they are applicable only in the presence of a magnetic field whose strength is small enough ($\alpha <1$) to prevent the formation of cluster chains. Furthermore, the equations retain their relevance when the rotation frequency $f_0$ does not exceed the Brownian relaxation frequency $2D_u$, in which case the relaxation orientation $\boldsymbol {m}_0$ would be restricted to the hydrostatic regime. These conditions can be summed up as the validity ranges of $\alpha <1$ and $w_0 \tau _{B}<<1$ for the proposed model.

Aqueous-based ferrofluids have a significantly shorter Brownian relaxation time, $\tau _{B}$, than oil-based ferrofluids due to the higher viscosity of the latter. Therefore, in the frequency range of $f = 50\unicode{x2013}200\ {\rm Hz}$ used in the experiments of Torres-Diaz et al. (2014) it is only relevant to compare the present model with the aqueous class of ferrofluids due to the fact that the magnitude of $w_0 \tau _{B}$ is less than 1. This condition, even for aqueous ferrofluids, must be accompanied by a restriction on the nanoparticle diameter or the intensity of the magnetic field strength to ensure that $\alpha <1$. Therefore, with the exception of WBF-3 $\alpha > 1$, the theoretical predictions of the developed model are only compared with the experimental results of WBF-1 and WBF-2, which satisfy the assumptions of $\alpha <1$ and $w_0 \tau _{B} << 1$.

Figure 5 shows a comparison of the liquid velocity predictions of the current model and the experimental results of Torres-Diaz et al. (2014) for WBF-1 and WBF-2. The liquid bulk corotates with the magnetic field in all cases considered. Furthermore, both theoretical and experimental results predict the same thickness of the velocity field's diffusive layer. This thickness, which measures the distance between the wall and the extremum point of the velocity profile, is approximately $40\,\%$ of the tube's radius $R$. This value remains constant regardless of the type of ferrofluid, WBF-1 or WBF-2, the magnetic field strength, $H_0$, or the rotation frequency, $f_0$. This result was also observed experimentally for the case of a concentrated WBF EMG705, $\epsilon _{p}> 1\,\%$, carried out by Chaves et al. (2008) in the same configuration considered by Torres-Diaz et al. (2014). This implies that the parameters of the external magnetic field and the type of ferrofluid do not affect the thickness of the liquid linear momentum diffusion layer. As mentioned above, in spin diffusion theory, the thickness of the velocity field diffusion layer is determined by the spin viscosity $\eta '$, since when $\eta ^{'} \rightarrow 0$, spin diffusion theory predicts no macroscopic flow in the ferrofluid. The value of $\eta '$, that scales as $\mu R_{p}^2$, predicts a very thin diffusion layer, of the same order of magnitude as the radius of a particle, which contradicts the experimental measurements of Chaves et al. (2008) and Torres-Diaz et al. (2014). According to Chaves et al. (2008), spin diffusion theory predicts the same diffusion layer only when $\eta '$ is artificially and arbitrarily inflated by 8 to 10 orders of magnitude compared with the prescription of Zaitsev & Shliomis (1969). The magnitude of such a substantial correction can only indicate a physical deficiency that is difficult to justify within the standard pseudohomogeneous ferrohydrodynamic model, and in any case calls into question its predictive character. It is important to note that despite the fact that the current model is independent of $\eta '$, it still predicts the thickness of the diffusion layer, as shown in figure 5. Therefore, it can be assumed that the origin of the macroscopic flow in the spin-up flow geometry does not need to invoke the spin viscosity concept and is independent of $\eta '$. Figure 5(a) shows the influence of the magnetic field strength on the velocity profile for WBF-1 over the range of the considered magnetic field strength $\mu _0 H_0=3.4\unicode{x2013}6.8\ {\rm mT}$, corresponding to $\alpha \lesssim 1$. In this figure, the predictions of the developed model and the experimental results clearly show that increasing the magnetic field intensity increases the liquid velocity. Regarding the effect of the rotational frequency of the magnetic field, figure 5(b) shows the velocity profile of WBF-1 for $f_0$ ranging between $75$ and $150\ {\rm Hz}$, $w_0 \tau _{B} << 1$. The experiment and the developed model show that increasing $f_0$ causes an increase in the liquid velocity. Figure 5(c) shows that for the case of ferrofluid WBF-2, the experiment and the current model quantitatively predict the same behaviour for the cases $\alpha =0.62,\ f=175\ {\rm Hz}$ and $\alpha =0.77,\ f=125\ {\rm Hz}$.

Figure 5. Comparison of experimental and theoretical predictions for the local velocity field along the azimuthal direction for the cases of ferrofluids WBF-1 and WBF-2: (a) WBF-1 – $f_0=150\ {\rm Hz}$; (b) WBF-1 – $\mu _0 H_0 = 5.6\ {\rm mT}$; (c) WBF-2.

Figure 5 shows that for all cases considered, the model is in global agreement with the experimental results. It is important to note that the measurement uncertainties recorded by Torres-Diaz et al. (2014) for the local velocity represent approximately $40\,\%$ of the mean field. For the WBF-1 cases, a maximum mean error of $15\,\%$ was recorded, with the predictions always remaining within the measurement uncertainties, except for the case $\alpha =1.07$, which can be considered outside the model validity range. As for the WBF-2 cases, a maximum mean error of $\alpha =1.07$ was recorded, which is significantly lower than that of WBF-1 and always remains within the measurement uncertainties.

The liquid in corotation with the magnetic field (figure 5) causes shear at the tube surface at $r = R$, thereby slowing down its rotation. This is tantamount to a torque exerted on the tube wall, which is another important hydrodynamic parameter of the ferrofluid in the spin-up flow geometry. Using the linear momentum equation for the liquid and the angular momentum equation for the discrete phase (5.8) and (5.11), the torque can be calculated from the total force exerted by the liquid at $r = R$ as follows:

(6.1)\begin{equation} L_w = 2 {\rm \pi}R^2 L \left(\frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) (\breve{H}_{\theta} m_{r} - \breve{H}_{r} m_{\theta})|_{r=R} - 2 \mu_e \varOmega|_{r=R}\right). \end{equation}

Figure 6 shows a comparison of the experimental torque results of Torres-Diaz et al. (2014) and those of the theoretical expression (6.1), for the cases WBF-1 and WBF-2. The experimental and model results show that the torque increases with increasing magnetic field strength and rotation frequency $f_0$, in the range where $\alpha < 1$ and $w_0 \tau _{B} << 1$. Figure 6 also shows that the evolution of $L_w$ appears to be linear with $f_0$ for WBF-1 and WBF-2. However, the evolution with $H_0$ is clearly nonlinear. Regarding the accuracy of the predictions of the theoretical model in comparison with the experimental results, a precise evaluation of the error is unfortunately not possible due to the unreported uncertainty estimation by the experimentalists. However, the model seems to agree with the experimental results in terms of the behaviour and the order of magnitude of the torque.

Figure 6. Comparison of Torres-Diaz et al. (2014) experimental results and model prediction with respect to the torque $L_w$ (6.1), for WBF-1 and WBF-2 ferrofluids: (a) WBF-1; (b) WBF-2.

6.2. Parametric study of spin-up flow physics under $\alpha <1$ and $w_0 \tau _{B} << 1$ conditions

The ability of the two-phase model to capture the basic mechanisms governing the ferrofluid spin-up flow has been demonstrated by comparing its predictions with the experiment of Torres-Diaz et al. (2014). However, to gain a comprehensive understanding of the process within the valid range of the model ($\alpha <1$ and $w_0 \tau _{B} << 1$), it is crucial to quantitatively assess the impact of each of the phenomena inventoried in the model. Therefore, a careful parametric study is essential to identify the influential factors that govern the process. The following simulation study will be performed on the same configuration used by Torres-Diaz et al. (2014), where the dimensions are $R=24.7$ and $L=63.5\ {\rm mm}$. The ferrofluid used will be water-based, consisting of cobalt ferrite nanoparticles, with a domain magnetization $M_{d}= 425\ {\rm kA}\ {\rm m}$. The primary objectives of the parametric study are to investigate the effects of the induced magnetic field, to evaluate the role of the magnetic Kelvin body force in nanoparticle transport, to analyse the effects of DDIs, and to investigate the contribution of the demagnetizing field. The systematic study of these parameters provides a more complete understanding of the dynamics underlying ferrofluid flow.

Figures 7 and 8 depict the dynamic evolution of the induced magnetic field, $\boldsymbol {\breve {H}}$, and the orientation of magnetic moments of particles, $\boldsymbol {m}$, in the $(r,\theta )$ plane. Due to the perpendicular orientation of the external field with respect to the $z$-axis, the components of $\boldsymbol {\breve {H}}$ and $\boldsymbol {m}$ along $z$ vanish consistently according to (5.24) and (5.25). Figure 7 highlights the influence of the intensity of the magnetic field $H_0$, which is encapsulated in the Langevin parameter $\alpha$, on the radial profiles of the azimuthal component of the induced magnetic field and the azimuthal and radial components of the nanoparticles magnetic moment orientation for a given rotational frequency of the magnetic field. Conversely, figure 8 demonstrates the impact of the rotational frequency of the external magnetic field on these same variables. Despite the uniform nature of the external magnetic field, the induced field exhibits localized gradients in all scenarios presented in figures 7 and 8. The magnitude of the induced $\breve {H}_r$ and $\breve {H}_{\theta }$ fields varies monotonically, exhibiting zero values at $r = 0$ and reaching their maximum at $r = R$. These findings underscore the significance of considering the induced field in the spin-up flow process.

Figure 7. Effect of an external magnetic field on the evolution of the induced local magnetic field and the particle magnetic moment orientations in the case of a WBF ($\epsilon _{{p},0}=0.001$, $f_0=100\ {\rm Hz}$, $D_{p}=15\ {\rm nm}$, $M_{d}=425\ {\rm kA}\ {\rm m}^{-1}$).

Figure 8. Effect of magnetic field frequency $f_0$, on the evolution of the radial profiles of (a) the (dimensionless) azimuthally induced magnetic field component and the orientation of the nanoparticle magnetic moment in terms of the (b) azimuthal and (c) radial components in the case of a WBF ($\epsilon _{{p},0}=0.001$, $D_{p}=15\ {\rm nm}$, $M_{d}=425\ {\rm kA}\ {\rm m}^{-1}$, $H_0=2.18\ {\rm kA}\ {\rm m}^{-1}$). In panel (c), the radial projection of the nanoparticle orientation via the Langevin equation (4.25) using the induced radial magnetic field is compared with $m_r$.

Figures 7(a) and 7(b) show the effect of the external magnetic field strength on $\breve {H}_{\theta }$ and $m_{\theta }$ for $\alpha =0.1, 0.5$ and $1$ when $f_0=100\ {\rm Hz}$. These figures show that the induced magnetic field is the cause for the magnetic moments of the particles to be oriented in the direction of $\theta$, because in the absence of $\breve {H}_{\theta }$ the particles are mostly oriented in the direction of $r$. This means that, despite the fact that $\breve {H}_{\theta }$ is omnipresent in the external field, disregarding the induced field does not yield a component $m_{\theta }$ due to the uniformity of the external field. These figures also show that the amplitudes of $\breve {H}_{\theta }$ and $m_{\theta }$ are affected by the external magnetic field intensity. Figure 7(b) shows that an increase in $\alpha$ causes an increase in $m_{\theta }$. However, for $\breve {H}_{\theta }$, figure 7(a) shows that the induced field is almost insensitive to the external magnetic field strength for $\alpha <1$, except in the region close to the wall, where $\breve {H}_{\theta }$ decreases with increasing $\alpha$ due to the effect of the demagnetizing field. Furthermore, $\breve {H}_{\theta }$ also becomes insensitive to $\alpha$ in the absence of the demagnetizing field, demonstrating that the azimuthally induced magnetic field component is generated by the rotation of the magnetic field. Figure 7(c) illustrates the impact of $\alpha$ on the evolution of $m_r$ and $\breve {H}_{r}$. It compares $m_r$ with $\epsilon _{p} ({\alpha }/{3}) \breve {H}_r$ to quantify the deviation from the Langevin relation (4.25). The Langevin relation describes the orientation of particle moments in a ferrofluid under hydrostatic equilibrium when subjected to a static magnetic field of equal strength instead of a rotating field. As shown in this figure, $m_r$ follows the same trend as $\breve {H}_{r}$, where its evolution is nearly linear. Furthermore, the difference caused by $f_0$ between $m_r$ and $\epsilon _{p} ({\alpha }/{3}) \breve {H}_r$ is vanishingly small. As a result, we can conclude that $m_r$ obeys Langevin's law with a small deviation caused by the hydrodynamic interaction between the liquid and the particles.

In terms of frequency influence, figure 8 illustrates the effect of $f_0$ on $\boldsymbol {\breve {H}}$ and $\boldsymbol {m}$ when $\alpha =0.5$ for $f_0$ values of $f_0=50,100$ and $150\ {\rm Hz}$. The rotation of the magnetic field induces a field along $\theta$, and the particles, being ferromagnetic, align with this field, resulting in a component $m_{\theta }$. Figures 8(a) and 8(b) show that increasing $f_0$ leads to an increase in the amplitudes of $\breve {H}_{\theta }$ and $m_{\theta }$. When there is no rotation of the magnetic field $(\,f_0 \rightarrow 0)$, $m_{\theta }$ tends to zero. Under the condition $w_0 \tau _{B} <<1$, the influence of $f_0$ on the radial profiles of $m_r$ and $\epsilon _{p} ({\alpha }/{3}) \breve {H}_r$ is negligible (figure 8c). It is worth noting that according to the current model's assumption of a uniform external field throughout the tube, the suppression of the induced field does not generate any flow in the ferrofluid. This can be seen from (5.8), where ${\boldsymbol{\nabla }} \boldsymbol {\breve {H}} \rightarrow 0$ transforms into a Poisson equation, ${\boldsymbol{\nabla }} \boldsymbol {v} \rightarrow 0$, with homogeneous boundary conditions, resulting in a trivial solution. This leads to the conclusion that the induced field serves as the driving force for the bulk motion of the ferrofluid, considering the continuity of its components at the tube surface $r = R$ with the external field $\boldsymbol {\breve {H}}_0$. As mentioned earlier, the rotation of the magnetic field at a given frequency causes the particles to rotate in the direction of the field, resulting in a component $m_{\theta }$ in the orientation of their magnetic moments. However, due to the viscous interaction between the particles and the liquid, the rotation of the magnetic moments of the particles lags behind the rotation of the magnetic field. The following expression can be used to determine the local lag angle between the direction of the magnetic moments of the particles and the magnetic field:

(6.2)\begin{equation} \varTheta = \arccos \left(\frac{\boldsymbol{\breve{H}} \boldsymbol{\cdot} \boldsymbol{m}}{|\boldsymbol{H}|\,|\boldsymbol{m}|}\right) .\end{equation}

Figure 9 provides insight into the influence of $\alpha$ (figure 9a) and $f_0$ (figure 9b) on the lag angle as estimated by (6.2). These figures show that the lag angle is uneven within the pipe in all cases considered. The maximum value of the lag angle occurs in the centre of the pipe where the induced field is the weakest. Conversely, the lag angle reaches its minimum value at the pipe surface, $r = R$, where the induced field is the strongest. Figure 9(a) illustrates the effect of $\alpha$ on $\varTheta$, showing an increase in the amplitude of $\varTheta$ with higher magnetic field strength. Regarding the effect of $f_0$, figure 9(b) shows that an increase in $f_0$ leads to a larger lag angle. These observations are consistent with expectations since higher values of $\alpha$ and $f_0$ result in increased viscous friction between the liquid and the particles. Under the conditions $\alpha <1$ and $w_0 \tau _{B} <<1$, the magnitude of $\varTheta$ remains small despite the radial gradients of the lag angle. Consequently, the effect of $\varTheta$ on the spin-up flow process is negligible. Another factor considered in this model is the magnetic force acting on the particles, commonly referred to as the Kelvin body force. It is worth noting that this force drives magnetophoresis, which causes particle transport within the tube, as described by the zeroth moment of the Smoluchowski (4.14). The components of the Kelvin body force are expressed along $r$ and $\theta$:

(6.3)\begin{equation} \boldsymbol{F}_{p}= \frac{K_{B} T}{V_{p}} \mathcal{F}_1 (\alpha) \begin{pmatrix} m_{r} \dfrac{{\rm d} \breve{H}_{r}}{{\rm d} r}-\dfrac{\breve{m}_{\theta} \breve{H}_{\theta}}{r} \\ m_{r} \dfrac{{\rm d} \breve{H}_{\theta}}{{\rm d} r} + \dfrac{\breve{m}_{\theta} \breve{H}_{r}}{r} \\ 0 \end{pmatrix}. \end{equation}

Figure 9. The effect of magnetic field intensity and rotation frequency on the lag angle in the case of a WBF ($\epsilon _{{p},0}=0.001$, $D_{p}=15\ {\rm nm}$, $M_{d}=425\ {\rm kA}\ {\rm m}^{-1}$): (a) $f_0=100\ {\rm Hz}$; (b) $\alpha =0.8$.

Figures 10(a) and 10(b) illustrate the evolution of the Kelvin force local components $F_{{p},r}$ and $F_{{p},\theta }$ for $\alpha$ values of $0.1$, $0.5$ and $1$. Notably, both figures clearly demonstrate the absence of magnetic field forces acting on the particles, indicating no translation of the particles within the tube. Consequently, particle transport within the tube is negligible, as depicted in figure 10(c), where no gradient of $\epsilon _{p}$ is observed for the three considered $\alpha$ values ($\alpha =0.1, 0.5, 1$). Therefore, both the particle transport and the Kelvin body force can be considered negligible in the spin-up flow process. However, unlike the Kelvin body force, it is important to emphasize that the torque exerted on the particles is not negligible. In the presence of the induced field, this torque introduces an asynchrony between the particle spin, $w$, and the liquid's vorticity, $\varOmega$, resulting in an antisymmetric stress. This antisymmetric stress is then converted into a linear momentum that drives the ferrofluid's movement (Zaitsev & Shliomis 1969; Rosensweig 2013). It is critical to note that ignoring the induced field results in the absence of asynchrony between the spin and vorticity of the fluid, effectively suppressing macroscopic flow within the ferrofluid. The magnetic torque can be calculated based on the net particle spin, $w-\varOmega$. By employing (5.11), the net particle spin, $w^{\ast }=(w-\varOmega )/w_0$, normalized with respect to the rotational pulsation of the magnetic field, $w_0$, can be expressed as

(6.4)\begin{equation} w^{{\ast}} = \frac{K_{B} T}{6 \mu \epsilon_{p} w_0 V_{p}} \mathcal{F}_1 (\alpha) (\breve{H}_{\theta} m_r - \breve{H}_{r} m_{\theta}).\end{equation}

Figure 10. Effect of magnetic field intensity on the Kelvin body force and particle transport, for $f_0=100\ {\rm Hz}$, in the case of WBF ($\epsilon _{{p},0}=0.001$, $D_{p}=15\ {\rm nm}$, $M_{d}=425\ {\rm kA}\ {\rm m}^{-1}$).

Figure 11 shows the effect of $\alpha$ and $f_0$ on the net dimensionless spin of the particles, $w^{\ast }$. It is evident that the particle spin in the centre of the tube tends to zero for all cases shown in this figure. This behaviour can be explained by referring to (6.4) which states that $w^{\ast } \rightarrow 0$ when $r = 0$, due to the rotational symmetry of the tube and the cancellation of the induced field at the centre. However, this does not imply that the particle spin is completely zero in this region, but rather indicates a synchronization between the particle spin and the vorticity of the liquid. The maximum value of $w^{\ast }$, which represents the peak asynchrony between $w$ and $\varOmega$, is observed near the wall where the induced field is at its maximum.

Figure 11. Effect of frequency and magnetic field on the net dimensionless particles spin $w^{\ast }$ in the case of a WBF ($\epsilon _{p}=0.001$, $D_{p}=15\ {\rm nm}$, $M_{d}=425\ {\rm kA}\ {\rm m}^{-1}$): (a) $f_0=100\ {\rm Hz}$; (b) $\alpha =0.5$.

To solve the conservation of angular momentum equation (5.11), no a priori knowledge of the wall boundary condition is required, as can be seen from the range $0 \leq r < R$ in figure 11 where $w^{\ast }$ is plotted. Mathematically, this equation does not require boundary conditions, as it is an algebraic equation that effectively bypasses the knowledge of any wall effect on $w$. From a phenomenological perspective, predicting the interactions between spinning particles and the wall at the continuum scale is challenging. Spin diffusion theory (Zaitsev & Shliomis 1969; Rosensweig 2013) considers a non-slip condition for $w$, that generates spin diffusion through the spin viscosity $\eta '$. Thus, Kaloni (1992) suggests describing $w$ at the wall with a non-slip condition that relates $w$ to the vorticity $\varOmega$ using a correction factor. In the present model, the antisymmetric stress caused by $w-\varOmega \neq 0$ is modelled by Faxén's laws (3.16), where the factor $6 \mu \epsilon _{p}$ represents the vortex viscosity in spin diffusion theory (Zaitsev & Shliomis 1969; Rosensweig 2013). Leach et al. (2009) performed an experimental study of the translational and rotational motion of a microsphere near a wall, highlighting the need for corrections to Faxén's laws to account for particle–wall interactions. However, these corrections depend on the particle–wall distance, which is difficult to predict. Despite uncertainties regarding particle slip or non-slip on walls, the present model accurately predicts the hydrodynamic behaviour of the ferrofluid, as demonstrated in § 6.1. In addition, since the ferrofluid is assumed to be dilute, the microrheological behaviour has a negligible influence on the hydrodynamic behaviour, allowing the wall effect to be omitted in the angular momentum balance. Lukaszewicz (1999) has also explored boundary condition issues in the context of polar fluid theory as a promising approach to the treatment of fluid flow with microrotational degrees of freedom. This approach deserves to be revisited in the future in the context of two-phase ferrofluid treatments of spin-wall boundary conditions, in order to assess its contribution to the understanding of wall–particle interaction boundary conditions, particularly in intense magnetic fields where chain-like aggregates can form. Indeed, the hydrodynamic interactions of these chain-like aggregates with the wall can be significant, depending on their size.

The parameterization of the magnetic field frequency on $w^{\ast }$ radial profiles results in a collapse to a single master curve (figure 11b) for $f_0=50, 100$ and $150\ {\rm Hz}$. This observation indicates that $w_0$ serves as a relevant scaling coefficient for $w-\varOmega$, implying that $6 \epsilon _{p} w_0 R$ acts as a velocity scaling factor in the linear momentum equation (5.8), as shown in figure 12(a). The convergence of the velocity profiles to a single master curve for the three frequency values considered (figure 12a) demonstrates the importance of $6 \epsilon _{p} w_0 R$ as a velocity scaling factor in terms of frequency influence. The influence of $\alpha$ on $w^{\ast }$ is illustrated in figure 11(a), which clearly shows that an increase in the magnetic field leads to an increase in the antisymmetric stress. This effect is further emphasized in figure 12(b), which illustrates the effect of $\alpha$ on the velocity profile. As $\alpha$ increases, the increase in antisymmetric stress leads to an increase in linear momentum, resulting in an increase in velocity. For $\alpha =0.1$, a tiny fraction close to zero of the rotational pulsation of the external field, $w_0$, is converted into antisymmetric stress, resulting in a velocity field close to zero. For $\alpha =0.5$ and $1$, approximately $4\,\%$ and $14\,\%$, respectively, of $w_0$ is converted to antisymmetric stress. This conversion manifests itself as an increase in velocity (figure 12b). However, the scaling factor $6 \epsilon _{p} w_0 R$ fails to unify the velocity profiles for the three $\alpha$ values considered (figure 12b). This discrepancy arises because the scaling factor does not account for $\alpha$. In addition, unlike frequency, the nonlinear dependence between $v$ and $\alpha$ makes it difficult to identify a scaling factor that characterizes their nonlinear relationship.

Figure 12. Effect of frequency and magnetic field on the velocity profile in the case of a WBF ($\epsilon _{p}=0.001$, $D_{p}=15\ {\rm nm}$, $M_{d}=425\ {\rm kA}\ {\rm m}^{-1}$): (a) $\alpha =0.5$; (b) $f_0=100\ {\rm Hz}$.

Figure 13 depicts the influence of DDIs on the liquid velocity profile. The comparison between velocity profiles with and without DDI is performed for a WBF to highlight its quantitative effect on the spin-up flow process. It is clear from all cases shown in the figure that the DDI contributes to the advective flow of the liquid rather than serving as a fundamental parameter to trigger the spin-up flow process. This observation can be attributed to the fact that the model consistently predicts macroscopic flow even in the absence of DDIs. The inclusion of DDIs in our model is achieved through mean field theory, which accounts for the influence of the dipolar interaction energy through the initial magnetic susceptibility $\xi$, as explained in § 4. Notably, this parameter is influenced by the mean particle distance, which is determined by the ferrofluid concentration $\epsilon _{p}$ and the average particle diameter $D_{p}$. Figure 13(a) illustrates the effect of particle concentration on the liquid velocity profiles in the presence and absence of DDI for particles with a diameter of $D_{p}=10\ {\rm nm}$. For $\epsilon _{p}=0.001$ and 0.002, the effect of DDI is negligible, except for a slight increase in velocity observed for $\epsilon _{p}=0.005$. The weak influence of DDI for particles with a diameter of $10\ {\rm nm}$ can be attributed to the fact that as the particle size decreases, the average distance between the particles increases, resulting in a weakening of the dipolar interaction energy. This weakening of the DDI is evident in the low initial magnetic susceptibility values of $\xi =0.0096,0.019$ and $0.048$, corresponding to concentration values of $\epsilon _{p}=0.001,0.002$ and $0.005$, respectively. Therefore, for ferrofluids with small particle diameters, DDIs are typically negligible, except in the case of highly concentrated suspensions $\epsilon _{p} > 1\,\%$, which fall outside the validity range of the present model. The DDI has a greater effect when the particle diameter is $D_{p}= 20\ {\rm nm}$ (figure 13b). As the particle concentration increases, the increase in the velocity field with DDI is more pronounced. Due to the increase in particle diameter, the initial magnetic susceptibility values increase, $\xi =0.077,0.15,0.39$, compared with the case where $D_{p}= 10\ {\rm nm}$. Therefore, as the particle diameter increases, the average distance decreases, which increases the dipolar interaction energy.

Figure 13. Effect of DDI on the velocity profile for $\alpha =0.8$ and $f_0=100\ {\rm Hz}$ in the case of water based ferrofluid ($M_{d}=425\ {\rm kA}\ {\rm m}^-$): (a) $D_{p}=10\ {\rm nm}$; (b) $D_{p}=20\ {\rm nm}$.

The last effect to be studied on the spin-up flow process is the demagnetizing field in terms of particle diameter and concentration (figure 14). It should be noted that the demagnetizing field is controlled by the geometrical parameter $\kappa$, which tends to $1/3$ for spherical particles. In the literature, the demagnetizing field is often pointed out as the physical phenomenon behind the spin-up flow process (Chaves et al. 2008; Torres-Díaz et al. 2012; Finlayson 2013). However, figure 14 demonstrates that even in the absence of the demagnetizing field $\kappa \rightarrow 0$, the present model still predicts a macroscopic flow in rotating magnetic fields. This indicates that the demagnetizing field, together with DDIs, contributes to the spin-up flow process rather than being the main driving parameter. The demagnetizing field tends to slow down the spin-up flow, leading to a reduction in the liquid velocity. This can be attributed to the decrease in the induced magnetic field caused by the magnetic moments at the particle surface (Joseph & Schlömann 1965). Regarding the influence of the concentration of the ferrofluid suspension on the demagnetizing field, the velocity profiles show an insensitivity to this phenomenon for a diameter $D_{p}=10\ {\rm nm}$ (figure 14a,c,e). However, for $D_{p}=20\ {\rm nm}$ (figure 14b,d,f), where the velocity field is reduced, the effect of the demagnetizing field becomes more pronounced. This reduction in the velocity profile is particularly significant for high $\epsilon _{p}$. The analysis of figure 14 reveals that the demagnetizing field is a volume-dependent process, where the reduction of the induced field leading to a reduction of the flow velocity is proportional to the particle volume. Therefore, it can be concluded that the demagnetizing field has a significant effect on particles with larger diameters.

Figure 14. Effect of particles’ demagnetization on the velocity profile for, $f_0=100\ {\rm Hz}$, in the case of WBF ($M_{d}=425\ {\rm kA}\ {\rm m}^{-1}$): (a) $\epsilon _{p}=0.001$ and $D_{p}=10\ {\rm nm}$; (b) $\epsilon _{p}=0.001$ and $D_{p}=20\ {\rm nm}$; (c) $\epsilon _{p}=0.005$ and $D_{p}=10\ {\rm nm}$; (d) $\epsilon _{p}=0.005$ and $D_{p}=20\ {\rm nm}$; (e) $\epsilon _{p}=0.01$ and $D_{p}=10\ {\rm nm}$; (f) $\epsilon _{p}=0.01$ and $D_{p}=20\ {\rm nm}$.

6.3. Cases when $\alpha >1$ and $w_0 \tau _{B} \sim 1$

As mentioned above, the validity of the model lies in the range where $\alpha <1$ and $w_0 \tau _{B} <<1$. However, in order to gain a comprehensive understanding of the physical phenomena involved, it is also necessary to investigate the behaviour of the model outside its range of validity. To achieve this, a comparison should be made between the velocity profiles predicted by the present model and those obtained from the experiment of Torres-Diaz et al. (2014) for cases where $\alpha >1$ and $w_0 \tau _{B} \sim 1$. The experiment performed by Torres-Diaz et al. (2014) on the WBF-1 ferrofluid showed that the velocity field increases with increasing $\alpha$, which is consistent with the theoretical predictions of our model. However, for $\alpha > 1$, a saturation of the velocity increase was observed, suggesting that the velocity tends to reach a limit. It has been observed that when the magnetic field is intense, the ferromagnetic suspension can form chain-like aggregates (Ivanov et al. 2004; Andreu et al. 2012; Faraudo et al. 2013; Zhao & Rinaldi 2018), which reduce the Brownian relaxation time and the magnetic torque exerted on the chains due to their increased length. Consequently, this effect is believed to be the cause of the velocity saturation observed for $\alpha > 1$ in the experiment of Torres-Diaz et al. (2014).

When chain-like aggregates form, Faxén's laws, (3.8) to (3.10), are no longer valid because they are based on the single-particle view. More precisely, (3.9) includes the term $6 \mu \epsilon _{p}$, which expresses the antisymmetric stress resulting from the asynchrony of the spin and the vorticity of the liquid, which, as discussed in § 6.2, is essential for the generation of the spin-up flow. This relation is similar to the vortex viscosity in spin diffusion theory. Consequently, a modification of this relation to account for chain formation is required to describe the macroscopic flow for $\alpha > 1$. Thus, it can be concluded that this model does not account for the effect of chain formation because its constitutive equations are based on the single-particle (or single-body) view. Accordingly, for $\alpha > 1$, figure 15(a) shows an overestimation of the experimental results by the theoretical predictions. However, it is important to note that the two-phase view used in the present model may provide a potential opportunity for investigation and development in cases where $\alpha > 1$. By studying the liquid and the discrete magnetic inclusions separately, it may be possible to model cause and effect separately. For example, this approach provides a particle transport equation with negligible effects for $\alpha < 1$, but it can be improved for $\alpha > 1$ by using a multibody approach to predict chain formation. Regarding the effect of frequency $f_0$ in cases where $w_0 \tau _{B} \sim 1$, figure 15(b) shows that the model overestimates the experimental results, as expected. This is due to the fact that the orientation of the magnetic moments at equilibrium has been modelled in the case of hydrostatic equilibrium. At high frequencies, however, the inertia of the fluid can affect the relaxation, and hydrostatic equilibrium is no longer valid.

Figure 15. Comparison of the present model with the experimental results of Torres-Diaz et al. (2014) for $\alpha >1$ and $w_0 \tau _{B} \sim 1$: (a) WBF1; (b) WGBF.

7. Concluding remarks

We have developed a two-phase, parameter-free, volume-averaged macroscopic approach to predictively describe the spin-up flow of dilute, cluster-free ferrofluids excited by low-frequency rotating magnetic fields. In contrast to the standard ferrohydrodynamic model, which treats the ferrofluid as a pseudohomogeneous phase, the proposed approach takes a different perspective, recognizing from the outset the inherent heterogeneity of two-phase colloidal magnetic suspensions by emphasizing the discrete nature of nanoparticles. Two fundamental aspects of the particle angular momentum equation have therefore been re-examined. First, spin diffusion, associated with a continuum-type molecular diffusion process, is excluded because the particles that make up the ferrofluid are very far apart, even on a microscopic scale. Second, the boundary condition of the standard ferrohydrodynamic model for the angular velocity of particles at the wall, which assumes a continuum-like non-slip condition on questionable grounds, was reconsidered. The two-phase approach revealed that the wall region must be treated with a slip condition in which nanoparticles can rotate close to the wall. This is due to the size of the particles, significantly larger than at the molecular scale, which phenomenologically invalidates the application of a no-slip boundary condition. The induced magnetic fields in the ferrofluid have a pronounced effect on the orientation of the magnetic moments and, consequently, have been treated inclusively by coupling using the Maxwell–Ampère and Maxwell-flux equations, taking into account the continuity of the normal and tangential components of the external field at the wall.

To account for the Brownian motion of nanoparticles, we combined macroscale two-phase transport equations with the Smoluchowski equation. This allowed us to derive two key equations: one describing the transport of the particle concentration field and another describing the transport of the average orientation of their magnetic moments. Notably, the resulting average orientation transport equation differs from the classical magnetization transport equation in the standard ferrohydrodynamic model, which completely neglects the need for a particle transport equation. In the case of dilute clusterless spin-up flow $\alpha < 1$, we observed limited spatial gradients in the particle concentration fields due to weak magnetophoresis resulting from the Kelvin body magnetic force. It is important to note that the mean orientation transport equation is only valid when $\alpha < 1$, since it relies on the truncated first-order expansion of the probability density. However, attempting to increase the order of the probability density expansion beyond the first moment is impractical for predicting spin-up flow when $\alpha > 1$. In such cases, the literature clearly shows the formation of cluster chains in the suspension. Therefore, outside the limit $\alpha < 1$, neither the Faxén nor the Smoluchowski equations, formulated under the assumption of suspensions consisting of single nanoparticles, are applicable for predicting the spin-up flow.

The proposed model was then used to simulate dilute ferrofluid flow in a cylindrical spin-up flow geometry induced by a rotating magnetic field. We first compared the model predictions with the experimental results of Torres-Diaz et al. (2014). The comparison showed that the theoretical predictions are in good agreement with the experimental results for the case of a weak magnetic field, where $\alpha <1$, and a rotation frequency much lower than the Brownian motion, where $w_0 \tau _{B} <<1$. An in-depth parametric study within the validity domain of the model was conducted to quantify and compare the relative weights of the different mechanisms highlighted by the model. This comparative study revealed that the induced field is the mechanism responsible for the process spin-up flow. The magnetic field gradients, combined with the torque generated by the asynchrony of the spinning particles and the liquid vorticity, produce an antisymmetric stress that is converted into linear momentum. In contrast, the DDIs and the demagnetizing field seem to play secondary roles rather than being the primary mechanisms prompting the spin-up flow. It has also been shown that the influence of the magnetic force is negligible compared with the magnetic torque acting on the particles. As a result, there are no significant particle transport or concentration gradients induced by magnetophoresis.

The parametric study demonstrated the importance of phase separation by using a two-phase approach at the continuum scale. When the magnetic field is rotating, the particles’ motion is triggered in the first place by the magnetic field excitation, since the liquid is not magnetic. On the other hand, the transition from a rotational to a linear motion for the liquid phase is a purely hydrodynamic process, obviously related to the effect of the induced magnetic field. A comparison of the experimental results outside the validity range of the model, where $\alpha > 1$ and $w_0 \tau _{B} \sim 1$, has also been carried out. The differences between the model predictions and the experimental results were also discussed. Two main points have been raised: the formation of chain-like aggregates in strong magnetic fields and the restriction of the relaxation to the hydrostatic regime. The importance of using a multiphase approach to describe the flow for $\alpha > 1$ has been emphasized in this discussion. This can be achieved by adopting a multibody approach using the particle transport equation, which is derived from the zeroth-order moment of the Smoluchowski equation, and by extending Faxén's laws to incorporate chain formations.