Particle-level Simulation of Magnetorheological Fluids: A Fully-Resolved Solver

1. Introduction

Magnetic particle suspensions, also known as magnetorheological fluids (MRFs), appear in a variety of applications [1,2]. In the traditional fluid engineering field, the magnetorheological effect has been applied to develop mechanical actuators and dampers [3,4]. In the newly emerged bio-engineering and drug delivery fields, there have been strong efforts to synthesize magnetic-based multifunctional particles [5,6]. In addition, in the field of natural resource and environmental engineering [7,8], precious metals or harmful substances dissolved in sea water are captured and recovered using magnetic particles subjected to an externally applied magnetic field.

When a magnetically polarizable particle is subjected to an externally applied magnetic field, they acquire dipole moments and become magnetized [9,10]. A magnetized particle starts interacting with neighboring magnetized particles leading to the formation of chain-like structures or clusters of particles aligned with the magnetic field direction (i.e., particle assembly) [11]. To date, numerous studies have investigated the dynamics of MRFs under magnetic fields. Hayes et al. [12] studied magnetic particles in microchannels by describing reversible self-assembled regularly spaced structures, when particles were exposed to an external magnetic field. From their study, they concluded that magnetic particles can be used in an extensive variety of on-chip applications and unique microfabrication techniques, automating the laboratory procedures. Melle and Martin [13] also developed a chain model for magnetorheological fluids in rotating magnetic fields. Through single-chain simulations as well as through experimental measurements, they showed that the chain shape and orientation depends strongly on the magnetic permeability of the particles ${\mu }_p$. Subsequently, Keaveny and Maxey [14] developed a finite-dipole model, where the magnetization of a particle is represented as a distribution of current density. This was proposed to estimate the magnetic forces between magnetic particles accurately and efficiently such that it can be applicable for systems with thousands of particles. In their model, the induced magnetization of a particle is represented as a localized Gaussian distribution of current that is added as a source term in the Poisson equation for the vector potential of the magnetic field [11]. The procedure yields very accurate solutions to collinear three-body problems. However, the scheme is not as accurate when considering other configurations with a large number of particles, because there is the need to include more information from the far field (e.g., quadrupole moments). Han et al. [15] presented a two-stage computational procedure for the numerical modelling of magnetorheological fluids. At the first stage, the particle dynamics is modelled using the discrete element method (DEM), whereas the hydrodynamic forces on the particles are approximated simply using the Stokes’ law (i.e., the fluid flow was not explicitly resolved) [16]. At the second stage, they deployed a combined approach using lattice Boltzmann method (LBM) and DEM to fully resolve the fluid fields, particle-particle, and particle–fluid hydrodynamic interactions. However, they raised an issue related to the accuracy of the magnetic interaction models while retaining the computational simplicity and efficiency. Subsequently, Ke et al. [17] developed a fully-resolved scheme based on lattice Boltzmann, immersed boundary, and discrete element methods (LBM-IBM-DEM) to simulate the behavior of magnetic particles moving in a fluid subject to an external magnetic field. The numerical results obtained showed that the LBM-IBM-DEM scheme was able to capture the major physical features of magnetic particle’s motion in a fluid. Specifically, they showed that particles first form fragmented chains along the magnetic direction. These chain-like clusters then continue to grow and align, and eventually, they approach an near steady state configuration. Additionally, it was shown that with the increase of the magnetic field a faster particle motion or merging between short chains occurs. Recently, Zhang et al. [18] developed a two-phase numerical simulation method using LBM-IBM-DEM approach to investigate the yielding phenomena during the start-up process of a MRF flowing through a microchannel under a transverse uniform magnetic field. The yielding of the MRF flowing through the microchannel was studied as a proxy to the deformation of the chains composed of magnetic particles. They showed that the yielding of a single-chain at different inlet velocities was regular. However, for a multi-chain system where chains are entangled, the yielding behavior presented an unpredictable regularity. Zhou et al. [19] also studied the motion of magnetic particles in a 3D microchannel flow modulated by the alternating gradient magnetic field. They used the LBM-IBM numerical simulation scheme, and showed that magnetic particles initially agglomerate due to their magnetic dipole force and then move together with the carrier fluid. They also showed that, in an alternating gradient magnetic field, magnetic particles oscillate along the flow direction, disturb the flow field, and increase the overall turbulence intensity. Leps and Hartzell [20] modeled the dynamics of MRFs using DEM method alone leveraging the open source LIGGGHTS [21] software. The algorithm is based on the mutual-dipole model to allow for the use of a large number of magnetic particles with several close neighbors while keeping a good trade-off between model accuracy and computational cost. Using accurate particle size distributions, high heritage contact models, and an uncoupled fluid model, Leps and Hartzell [20] were able to match the experimentally derived yield stress results for MRFs more closely than using mono-disperse particle size distributions. Lastly, Tajfirooz et al. [22] presented an Eulerian–Lagrangian approach for simulating the magneto-Archimedes separation of neutrally buoyant non-magnetic spherical particles within MRFs. A four-way coupled point-particle method [23,24] was employed, where all relevant interactions between an external magnetic field, a magnetic fluid and immersed particles were taken into account. First, the motion of rigid spherical particles in a magnetic liquid was studied in single- and two-particle systems. It was shown that numerical results of single- and two-particle configurations were in good agreement with detailed experimental results on particle position. Subsequently, the magneto-Archimedes separation of particles with different mass densities in many-particle systems interacting with the fluid was also studied. It was concluded that history effects and inter-particle interactions significantly influence the levitation dynamics of particles and have a detrimental impact on the separation performance.

Most of the aforementioned numerical studies around MRFs focus on the formation of magnetorheological structures using the simplified Stokes drag law and the dipole–dipole interaction model, excluding the hydrodynamic interactions between particles and higher order mutual magnetic interactions. The flow characteristics and chain formation features induced by coupled hydrodynamic and magnetic interactions are still missing in the literature. This is mainly due to the lack of proper numerical models that can take into account both inter-particle magnetic and hydrodynamic interactions, in addition to other relevant attributes (e.g., particle type, size, etc.), in a fully coupled algorithm.

In this work, we develop a fully-resolved simulation algorithm using a combination of the finite-volume, immersed boundary and discrete element methods to couple both hydrodynamic and magnetic interactions among magnetic particles suspended in Newtonian fluids. The newly-developed algorithm, so-called FVM-IBM-DEM-MAG solver, is able to describe flows with suspended magnetic particles immersed in a fluid subject to an external magnetic field. The magnetic force exerted on the particles is computed using the gradient of the magnetic field strength, which is obtained from the imposed external magnetic field [17]. The magnetic interactions between the particles are implemented using a mutual dipole model [20] allowing the magnetic fields of other particles to contribute to the magnetization and motion of the particle under consideration. The presented numerical algorithm has several advantages, specifically: (i) it is based on open-source libraries, OpenFOAM and LIGGGHTS, which allows the extension of the algorithm for other applications (e.g., simulation of viscoelastic fluids with suspended magnetic particles); and (ii) it employs a direct particle-level simulation methodology to resolve both hydrodynamic and magnetic interactions allowing accurate predictions of the flow patterns and particle assembly. We focus on simulations of spherical particles suspended in a Newtonian fluid in order to introduce the numerical algorithm and study its feasibility for extension to more complex flows, involving fluids with non-linear rheological behavior, and also particle with different shapes.

The remainder of this work is structured as follows. In Section 2, we present the underlying physics and mathematical formulation describing the motion of magnetic particles in a Newtonian fluid. In Section 3, we present the particle-level numerical methodology leading to the FVM-IBM-DEM-MAG solver that couples the continuum and discrete phases in MRFs. In Section 4, we present four case studies with different level of complexities to test the developed algorithm, namely the motion and interaction of two magnetic spheres settling in an incompressible Newtonian fluid under external magnetic field, and the 2D and 3D flow behaviors of random arrays of magnetic spheres immersed in an incompressible Newtonian fluid. Finally, in Section 5, we summarize the main conclusions of this work.

2. Underlying Physics

The magnetorheological fluids (MRFs) considered in this study contain micro-scale magnetic particles with no-Brownian motion suspended in a non-magnetic incompressible Newtonian carrier fluid. MRFs deform and self-organize into mesoscopic structures depending on the internal (e.g., particle concentration) and external stimuli (e.g., temperature, flow, and magnetic fields). Among these stimuli, the application of magnetic fields is shown to provide instant action and contactless control of the mesoscopic physical structures, causing a reversible transition from a fluid-like to a solid-like state. When subjected to an external magnetic field, particle assembly occurs that provides the fluid with the ability to transmit force. In that state, the effective viscosity of the fluid increases to the extent of becoming a viscoelastic solid. The particle assembly promoted by magnetic field can be controlled, i.e., destroyed, deformed, or delayed. To accurately predict the particle assembly and chain formation in MRFs, the coupled interactions between the magnetic field, fluid, and particles must be resolved. The dynamics of MRFs, thus, present a multi-physics problem across different scales. In this section, we present the underlying physics governing the dynamics of MRFs made of rigid micro-scale magnetic spheres suspended in non-magnetic incompressible Newtonian carrier fluids under a static magnetic field.

2.1. Magnetostatic fields

Macroscopic electromagnetic phenomena are described using Maxwell’s fundamental equations [1,9]. In this study, we assume the quantities of interest (e.g., magnetic field strength) do not vary with time, and there is no interaction between electric and magnetic fields. Therefore, we can decouple electrostatic and magnetostatic fields, and consider the problem of magnetostatic field with no free electric currents. The Maxwell’s equations for magnetostatic cases reduce to,

$$\begin{array}{c}\hfill \nabla ·\mathbf{B}=0,\end{array}$$

(1)

$$\begin{array}{c}\hfill \nabla \times \mathbf{H}=0,\end{array}$$

(2)

where $\mathbf{B}$ is the magnetic flux density, and $\mathbf{H}$ is the magnetic field strength. Here ∇ denotes the gradient operator, $\nabla ·$ denotes the divergence operator, and $\nabla \times$ denotes the curl tensor operations. For a linear isotropic domain (matrix) with a constant magnetic permeability, $\mu$, the constitutive equation relating the two field quantities, $\mathbf{B}$ and $\mathbf{H}$, reads as,

$$\begin{array}{c}\hfill \mathbf{B}=\mu \mathbf{H},\end{array}$$

(3)

where

$$\begin{array}{c}\hfill \mu =\left\{\begin{array}{cc}{\mu }_p\hfill & \mathrm{in}\mathrm{the}\mathrm{particle}\mathrm{domain},\hfill \\ {\mu }_f\hfill & \mathrm{in}\mathrm{the}\mathrm{fluid}\mathrm{domain},\hfill \end{array}\right.\end{array}$$

(4)

with ${\mu }_p$ and ${\mu }_f$ denoting the particles and base fluid’s magnetic permeability, respectively. Notice that $\mu$ is discontinuous at fluid–particle interfaces, and, therefore, should be evaluated by following a similar interpolation of material properties as the one used in the level set method [11,25]. Hereafter, consider that the total computational domain is represented by $\Omega$ $={\Omega }_s\cup {\Omega }_f$, where ${\Omega }_s$ is solid (“solid particles") domain, and ${\Omega }_f$ is the fluid domain. The total domain, solid and fluid boundaries are represented by $\partial \Omega$, $\partial {\Omega }_s$ and $\partial {\Omega }_f$, respectively.

To solve the first-order differential equations involving the two magnetic field quantities, Eqs. (1) and (2), we first convert them into a second-order differential equation involving only one magnetic field quantity. For that purpose, Eq. (2) admits the existence of a scalar potential, $\varphi$, such that,

$$\begin{array}{c}\hfill \mathbf{H}=-\nabla \varphi ,\end{array}$$

(5)

which can be substituted into Eq (1) with the aid of Eq. (3) to yield the following second-order differential equation,

$$\begin{array}{c}\hfill {\nabla }^2\left(\mu \varphi \right)=0.\end{array}$$

(6)

2.1.1. Magnetic forces and torques

In order to describe the particle motion and the flow around it influenced by a magnetic field, a relationship between the applied magnetic field and the resultant force acting on the particles is needed. This force, known as magnetostatic polarization force [17], on particle i can be evaluated as [17],

$$\begin{array}{c}\hfill {\mathbf{F}}_i^{me}={\int }_{{\Omega }_s}\left({\mu }_f{\chi }_{e}H\nabla H\right)d{\Omega }_s,\end{array}$$

(7)

where ${\chi }_e$ stands for the particle’s magnetic susceptibility given by the Clausius-Mossotti relationship [26],

$$\begin{array}{c}\hfill {\chi }_e=\frac{3{\mu }_p}{3+{\mu }_p}.\end{array}$$

(8)

The torque generated by the magnetostatic polarization force on particle i is computed as,

$$\begin{array}{c}\hfill {\mathbf{T}}_i^{me}={\int }_{{\Omega }_s}\left({\mu }_f{\chi }_{e}H\times H\right)d{\Omega }_s.\end{array}$$

(9)

Another fundamental force in MRFs is the magnetic dipoles evidenced by particles with opposite magnetic point poles [27]. Therefore, in MRFs, particle motion is affected not only by an external magnetic field, but also by other nearby magnetized particles since each particle has a permanent magnetic moment, $\mathbf{m}$. The dipole-dipole interactions between particles i and j results in dipole-dipole inter-particle magnetic force $\left({\mathbf{F}}_{ij}^{d-d}\right)$ and torque $\left({\mathbf{T}}_{ij}^{d-d}\right)$ that are calculated using the dipole-dipole contact model [28,29] as,

$$\begin{array}{c}\hfill {\mathbf{F}}_{ij}^{d-d}=\frac{3}{r^5}\left[\left({\mathbf{m}}_i·{\mathbf{m}}_j\right)\mathbf{r}-\frac{5}{r^2}\left({\mathbf{m}}_i·\mathbf{r}\right)\left({\mathbf{m}}_j·\mathbf{r}\right)\mathbf{r}+\left({\mathbf{m}}_j·\mathbf{r}\right){\mathbf{m}}_i+\left({\mathbf{m}}_i·\mathbf{r}\right){\mathbf{m}}_j\right],\end{array}$$

(10)

and

$$\begin{array}{c}\hfill {\mathbf{T}}_{ij}^{d-d}=-\frac{1}{r^3}\left[\left({\mathbf{m}}_i\times {\mathbf{m}}_j\right)-\frac{3}{r^2}\left({\mathbf{m}}_j·\mathbf{r}\right)\left({\mathbf{m}}_i\times \mathbf{r}\right)\right],\end{array}$$

(11)

where ${\mathbf{m}}_i$ and ${\mathbf{m}}_j$ are the magnetic moment vectors of the two particles, $\mathbf{r}$ is the separation vector between the two particles, and r is the magnitude of the separation vector $\mathbf{r}$. For MRFs consisting of N particles, a direct evaluation of the dipole-dipole interaction alone is $O\left(N^2\right)$ operations. This puts a severe computational constraint on the number of particles that can be simulated with a direct computation of the inter-particle dipole-dipole force. To compute the magnetic moment of each particle, $\mathbf{m}$, the fixed dipole model [15] or the mutual dipole model [20] can be used. In dilute MRFs (i.e., low concentration of magnetic particles), it is often acceptable to use the magnetic moment calculated from the background magnetic field (i.e., fixed dipole model). However, in concentrated MRFs (i.e., high concentration of magnetic particles), the induced magnetic fields from magnetized neighboring particles begin to have a significant effect on the particles magnetic moment vector. Therefore, for accuracy in force calculations, a more complex model (i.e., mutual dipole model) should be leveraged.

2.1.2. Fixed dipole model

When the effect of the extra magnetic field generated by neighboring magnetized particles is negligible on particles dynamics, then it is safe to assume that any particle is theoretically magnetized only by the externally applied magnetic field. Therefore, each particle is considered as a point dipole and the magnetic force between the particles are pairwise only. In this model, the magnetic moment of particle i is given by [15],

$$\begin{array}{c}\hfill {\mathbf{m}}_i=4\pi r^3\frac{\chi -1}{\chi +2}{\mathbf{H}}_0,\end{array}$$

(12)

where $\chi ={\mu }_p/{\mu }_f$ is the relative susceptibility of the particles over the carrier fluid, and ${\mathbf{H}}_0$ is the magnetic field strength of the externally applied uniform magnetic field. Notice that as the carrier fluid is assumed to be non-magnetic, its permeability is the same as that of a vacuum, i.e., ${\mu }_f={\mu }_0=4\pi \times {10}^{-7}$ [Tm/A], where T is Tesla, m is meter, and A is Ampere. This model is accurate when the two particles are far apart, and it loses accuracy when the separation distance of the particles decreases. The accuracy of the model also depends on the relative susceptibility, $\chi$. It has been shown by Keaveny and Maxey [14] that, at $\chi =5$, the fixed dipole model underestimates the maximum attractive force by around $35\%$, whereas overestimates the maximum repulsive force by $50\%$ or more, and the errors increase for larger $\chi$ values.

2.1.3. Mutual dipole model

The mutual dipole model [14] allows for the magnetic fields of the neighboring particles to contribute to the magnetization of the particle under consideration. A particle, thus, is subjected not only to the primary magnetization due to the external magnetic field, but also to a secondary magnetization from the other particles’ magnetic fields. Considering the mutual magnetization of N magnetizable particles with their centres at ${\mathbf{x}}_i$$(i=1,\cdots ,N)$ in a uniform magnetic field with strength ${\mathbf{H}}_0$, the magnetic moment of the particle i, ${\mathbf{m}}_i$, is given by [20],

$$\begin{array}{c}\hfill {\mathbf{m}}_i=4\pi r^3\frac{\chi -1}{\chi +2}\left[{\mathbf{H}}_0+\mathbf{H}\left({\mathbf{x}}_i\right)\right](i=1,\cdots ,N),\end{array}$$

(13)

where $\mathbf{H}\left({\mathbf{x}}_i\right)$ represents the total secondary magnetic field strength generated by other magnetized particles. The total secondary magnetic field strength can be expressed as [20],

$$\begin{array}{c}\hfill {\displaystyle \mathbf{H}\left({\mathbf{x}}_i\right)=\sum _{j,j\ne i}^N{\mathbf{H}}_j({\mathbf{m}}_j,{\mathbf{r}}_{ij})=\sum _{j,j\ne i}^N\frac{1}{4\pi }\frac{3{\widehat{\mathbf{r}}}_{ij}({\mathbf{m}}_j·{\widehat{\mathbf{r}}}_{ij})-{\mathbf{m}}_j}{r_{ij}^3},}\end{array}$$

(14)

with ${\mathbf{r}}_{ij}={\mathbf{x}}_i-{\mathbf{x}}_j$, $r_{ij}=\left|{\mathbf{r}}_{ij}\right|$, and ${\widehat{\mathbf{r}}}_{ij}={\mathbf{r}}_{ij}/r_{ij}$. Once the ${\mathbf{m}}_i$ values are computed for all particles, the inter-particle dipole-dipole force and torque between any two pairs are obtained using Eqs. (10) and (11), respectively.

2.2. Incompressible fluid flow

In the MRFs considered in this study, the carrier fluid is considered to be a non-magnetic incompressible Newtonian fluid. The governing equations for the flow of these fluids consist of the continuity equation,

$$\begin{array}{c}\hfill \nabla ·\mathbf{u}=0\mathrm{in}{\Omega }_f,\end{array}$$

(15)

and the Cauchy momentum equation,

$$\begin{array}{c}\hfill {\rho }_f\left(\frac{\partial }{\partial t}+\mathbf{u}·\nabla \right)\mathbf{u}=-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\mathrm{in}{\Omega }_f.\end{array}$$

(16)

Here ${\rho }_f$ and $\mathbf{u}$ are the fluid density and velocity vector, respectively, t is the time, p is the pressure, and ${\eta }_S$ is the viscosity of the Newtonian fluid. To complete the strong mathematical form describing the flow of MRFs, the following initial and boundary conditions are considered,

$$\left\{\begin{array}{c}\mathbf{u}(\mathbf{x},t=0)={\mathbf{u}}_0\left(\mathbf{x}\right)\mathrm{in}{\Omega }_f,\hfill \\ \mathbf{u}(\mathbf{x},t)={\mathbf{u}}_{\partial \Omega }\mathrm{on}\partial {\Omega }_f,\hfill \\ \mathbf{u}(\mathbf{x},t)={\mathbf{u}}_i\mathrm{on}\partial {\Omega }_s,\hfill \\ \left(-p\mathbf{I}+{\eta }_S\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)\right)·\widehat{\mathbf{n}}={\sigma }_{\partial {\Omega }_s}\mathrm{on}\partial {\Omega }_s.\hfill \end{array}\right.$$

(17)

In Eq. (17), $\widehat{\mathbf{n}}$ is the outward normal unit vector to $\partial {\Omega }_s$, ${\sigma }_{\partial {\Omega }_s}$ is the stress vector acting from the fluid on the solid body surface, and ${\mathbf{u}}_i$ is the (unknown) velocity of the solid-fluid interface. The initial velocity ${\mathbf{u}}_0$ is required to satisfy Eq. (15), and the boundary velocity ${\mathbf{u}}_{\partial \Omega }$ should satisfy the compatibility condition (last equation in Eq. 17) at all times.

The motion of magnetic particles is strongly affected by short-range and long-range hydrodynamic forces (drag, lift, etc.), and the resultant torques, when they are dispersed in a viscous incompressible fluid. The hydrodynamic force acting on the surface of particle i can be obtained using [30,31],

$$\begin{array}{c}\hfill {\mathbf{F}}_i^h={\int }_{{\Omega }_s}\left(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\right)d{\Omega }_s.\end{array}$$

(18)

The resultant hydrodynamic torques on particle i, denoted by ${\mathbf{T}}_i^h$, can be then calculated by taking the cross product between the position vector $\mathbf{r}$ (pointing from the fluid cell centroid to the particle centroid) and the total force from Eq. (18) that reads as,

$$\begin{array}{c}\hfill {\mathbf{T}}_i^h={\int }_{{\Omega }_s}\left\{\mathbf{r}\times \left(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\right)\right\}d{\Omega }_s.\end{array}$$

(19)

The force contribution arising from pressure does not give rise to any torque contribution, due to symmetry of spherical magnetic particles. Thus, normal forces acting perpendicular to the particle surface, such as pressure, do not induce any torque. This is not the case if particle shape departs from the spherical shape (e.g., spheroids). In MRFs, particles also experience the buoyancy force, denoted by ${\mathbf{F}}_i^g$, which is given by the weight of the displaced fluid. The buoyancy force can be calculated as,

$$\begin{array}{c}\hfill {\mathbf{F}}_i^g={\int }_{{\Omega }_s}\left({\rho }_f\mathbf{g}\right)d{\Omega }_s,\end{array}$$

(20)

where $\mathbf{g}$ is the gravitational acceleration vector.

2.3. Particle transient motion

The transient motion of dispersed magnetic particles (i.e., solid phase), can be modeled using the Newton’s second law of motion as,

$$\begin{array}{c}\hfill m_i\frac{d{\mathbf{U}}_i^p}{dt}=\sum _{j=1}^{n_i^c}{\mathbf{F}}_{ij}^c+\sum _{j=1}^{n_i^{cut}}{\mathbf{F}}_{ij}^{d-d}+{\mathbf{F}}_i^{me}+{\mathbf{F}}_i^h+{\mathbf{F}}_i^g,\end{array}$$

(21)

and

$$\begin{array}{c}\hfill I_i\frac{d{\omega }_i^p}{dt}=\sum _{j=1}^{n_i^c}{\mathbf{T}}_{ij}^c+\sum _{j=1}^{n_i^{cut}}{\mathbf{T}}_{ij}^{d-d}+{\mathbf{T}}_i^{me}+{\mathbf{T}}_i^h,\end{array}$$

(22)

for the conservation of linear and angular momentum of the particle i with mass $m_i$ and moment of inertia $I_i$, respectively. Here, ${\mathbf{U}}_i^p$ and ${\omega }_i^p$ denote the translational and angular velocities of particle i, respectively, ${\mathbf{F}}_{ij}^c$ and ${\mathbf{T}}_{ij}^c$ are the contact force and contact torque resulting from the particle-particle and particle-wall interactions (with the number of total contacts, $n_i^c$, for particle i) that can be calculated using different contact models [32,33], ${\mathbf{F}}_{ij}^{d-d}$ and ${\mathbf{T}}_{ij}^{d-d}$ are the dipole-dipole inter-particle magnetic force and torque for a number of $n_i^{cut}$ possible interactions in the admissible cut-off region, respectively, ${\mathbf{F}}_i^{me}$ and ${\mathbf{T}}_i^{me}$ are the magnetostatic polarization force and torque due to the external magnetic field, respectively, ${\mathbf{F}}_i^h$ and ${\mathbf{T}}_i^h$ are the hydrodynamic force and torque acting on particle i, respectively, and ${\mathbf{F}}_i^g$ is the buoyancy force.

We leverage DEM, developed by Cundall and Strack [34] and implemented in LIGGGHTS open-source library [21], to model the transient motion of dispersed magnetic particles described by Eqs. (21) and (22). In DEM, multiple search algorithms are employed to identify contacting pairs of discrete particles [35], and different contact models are developed to integrate various mechanisms and effects such as elasticity, plasticity, viscoelasticity, friction, cohesion, damage, fracture, etc. in the contact points [21]. In this study, we adopted the spring-dashpot contact model that can be extended to other non-linear models depending on the chosen stiffness and damping parameters as function of the particle overlap displacement [33]. In this model, the total contact force between particle i and particle j is calculated using [36],

$${\mathbf{F}}_{ij}^c={\left({\mathbf{F}}_{ij}^c\right)}_n+{\left({\mathbf{F}}_{ij}^c\right)}_t,$$

(23)

where ${\left({\mathbf{F}}_{ij}^c\right)}_n$ is the normal contact force,

$${\left({\mathbf{F}}_{ij}^c\right)}_n=-k_n{\delta }_n\mathbf{n}-{\gamma }_n{\left({\mathbf{U}}_{ij}^p\right)}_n,$$

(24)

and ${\left({\mathbf{F}}_{ij}^c\right)}_t$ is the tangential contact force,

$${\left({\mathbf{F}}_{ij}^c\right)}_t=\min \left(-k_t{\delta }_t-{\gamma }_t{\left({\mathbf{U}}_{ij}^p\right)}_t,{\beta }_s\left|{\left({\mathbf{F}}_{ij}^c\right)}_n\right|\frac{{\delta }_t}{|{\delta }_t|}\right),$$

(25)

with

$${{\delta }_t}^{\left(n\right)}={{\delta }_t}^{(n-1)}+{\left({\mathbf{U}}_{ij}^p\right)}_t\Delta t.$$

(26)

In Eqs. (24), (25) and (26), $\mathbf{n}$ is the unit vector in the normal direction, $k_n$ and $k_t$ are the elastic stiffness for normal and tangential contacts, respectively, ${\gamma }_n$ and ${\gamma }_t$ denote the damping coefficients in normal and tangential directions, respectively, ${\delta }_n$ is normal overlap displacement between two particles, ${\left({\mathbf{U}}_{ij}^p\right)}_n$ and ${\left({\mathbf{U}}_{ij}^p\right)}_t$ are relative velocities in normal and tangential directions of particle i relative to particle j, respectively, with the relative velocity defined as ${\mathbf{U}}_{ij}^p={\mathbf{U}}_i^p-{\mathbf{U}}_j^p$, ${\beta }_s$ is the sliding friction coefficient, ${\delta }_t^{\left(n\right)}$ and ${\delta }_t^{(n-1)}$ are the tangential overlap at the current and previous step, and $\Delta t$ is the time step. The resultant contact torque on particle i due to its contact with particle j, denoted by ${\mathbf{T}}_{ij}^c$, can be then calculated by taking the cross product between the total contact force from Eq. (23) and the position vector leading to,

$${\mathbf{T}}_{ij}^c={\mathbf{F}}_{ij}^c\times ({\mathbf{x}}_c-{\mathbf{x}}_i),$$

(27)

where ${\mathbf{x}}_c$ and ${\mathbf{x}}_i$ are the position of contact point and particle i centroid, respectively.

3. Numerical Methodology

This section presents the numerical formulation for an algorithm using the FVM, IBM and DEM that is able to efficiently handle the rigid body motion of magnetic spherical particles surrounded by a Newtonian fluid. The algorithm considers a fictitious domain formulation, which provides a rigorous basis for the immersed boundary (IB) implementation performed in the open source framework code $CFDEMcoupling$ [31,37,38]. The open source IB solver originally developed by Hager et al. [37] is modified and improved for this study to take into account both hydrodynamic and magnetic interactions between the fluid continuum phase and the particulate disperse phase in a fully coupled manner. Algorithm 1 summarizes the so-called FVM-IBM-DEM-MAG solver describing the solution procedure of the fluid phase and magnetic field equations, the DEM approach to handle the particle’s motion, and the IBM scheme to fully couple the continuum phase with the particulate phase.

Algorithm 1:Fully-resolved FVM-IBM-DEM-MAG algorithm to model magnetorheological fluids

step 1: at time $t=0$ (a) Set initial and boundary conditions (b) Send initial particle position and velocities to CFD solver from DEM solver step 2: at time $t=t+\Delta t$ (a) Compute particle volume fraction (b) Dynamic mesh refinement (c) Calculate loads on particles (hydrodynamic and magnetic external forces, ${\mathbf{F}}_i^h$ and ${\mathbf{F}}_i^{me}$, respectively, and torques, ${\mathbf{T}}_i^h$ and ${\mathbf{T}}_i^{me}$, respectively; particle-particle contact force and torque, ${\mathbf{F}}_{ij}^c$ and ${\mathbf{T}}_{ij}^c$, respectively; magneto dipole-dipole force and torque, ${\mathbf{F}}_{ij}^{d-d}$ and ${\mathbf{T}}_{ij}^{d-d}$, respectively; and buoyancy force ${\mathbf{F}}_i^g$, etc.) given by ${\mathbf{F}}_i^h={\sum }_{c\in {\overline{T}}_h}(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u})\left(c\right)·V\left(c\right){\mathbf{F}}_i^{me}={\sum }_{c\in {\overline{T}}_h}({\mu }_0{\chi }_e\mathbf{H}\nabla \mathbf{H})\left(c\right)·V\left(c\right)$ ${\mathbf{T}}_i^h={\sum }_{c\in {\overline{T}}_h}\left[\mathbf{r}\left(c\right)\times (-\nabla p+{\eta }_S{\nabla }^2\mathbf{u})\left(c\right)\right]·V\left(c\right){\mathbf{T}}_i^{me}={\sum }_{c\in {\overline{T}}_h}({\mu }_0{\chi }_e\mathbf{H}\times \mathbf{H})\left(c\right)·V\left(c\right)$ ${\mathbf{F}}_{ij}^c$ and ${\mathbf{T}}_{ij}^c$ are calculated using the non-linear elastic Hertz-Mindlin contact model. ${\mathbf{F}}_{ij}^{d-d}=\frac{3}{r^5}\left[\left({\mathbf{m}}_i·{\mathbf{m}}_j\right)\mathbf{r}-\frac{5}{r^2}\left({\mathbf{m}}_i·\mathbf{r}\right)\left({\mathbf{m}}_j·\mathbf{r}\right)\mathbf{r}+\left({\mathbf{m}}_j·\mathbf{r}\right){\mathbf{m}}_i+\left({\mathbf{m}}_i·\mathbf{r}\right){\mathbf{m}}_j\right]$ ${\mathbf{T}}_{ij}^{d-d}=-\frac{1}{r^3}\left[\left({\mathbf{m}}_j\times {\mathbf{m}}_i\right)-\frac{3}{r^2}\left({\mathbf{m}}_i·\mathbf{r}\right)\left({\mathbf{m}}_j\times \mathbf{r}\right)\right]$ ${\mathbf{F}}_i^g={\sum }_{c\in {\overline{T}}_h}\left({\rho }_f\mathbf{g}\right)\left(c\right)·V\left(c\right)$ (d) Solve Newton-Euler equations (Velocity-Verlet integration) to obtain new particle position, and linear and angular velocities ($\mathrm{in}{\Omega }_p$) $m_i\frac{d{\mathbf{U}}_i^p}{dt}={\sum }_{j=1}^{n_i^c}{\mathbf{F}}_{ij}^c+{\sum }_{j=1}^{n_i^{cut}}{\mathbf{F}}_{ij}^{d-d}+{\mathbf{F}}_i^{me}+{\mathbf{F}}_i^h+{\mathbf{F}}_i^g{I}_i\frac{d{\omega }_i}{dt}={\sum }_{j=1}^{n_i^c}{\mathbf{T}}_{ij}^c+{\sum }_{j=1}^{n_i^{cut}}{\mathbf{T}}_{ij}^{d-d}+{\mathbf{T}}_i^{me}+{\mathbf{T}}_i^h$ (e) Solve fluid governing equations subjected to an external magnetic field ($\mathrm{in}{\Omega }_f$) ${\nabla }^2\left(\mu \varphi \right)=0$ $\nabla ·\mathbf{u}=0$ ${\rho }_f\left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}·\nabla \mathbf{u}\right)=-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}$ (f) Impose the rigid-body motion of the particles on the fluid velocity field (g) Correct velocity and pressure fields

At time $t=0$, the fluid and particle initial and boundary conditions are read from the case study input files (step 1(a) in Algorithm 1). Additionally, the DEM solver sends the particle initial position and velocities to the CFD solver (step 1(b) in Algorithm 1). At time $t=t+\Delta t$, the numerical algorithm starts with the location of the magnetic particles, saving the cell ID of the centre position of each particle. This procedure, then, allows to compute the particle volume fraction in each cell (step 2(a) in Algorithm 1). Subsequently, as shown in Figure 1, the algorithm uses dynamic mesh refinement near the particles’ surface ($\partial {\Omega }_s$) to accurately capture the fluid (domain ${\Omega }_f$) forces developed on those regions (step 2(b) in Algorithm 1).

Using the fluid solution from the last time-step in the regions marked by the particle volume fraction, the hydrodynamic, magnetostatic polarization, and buoyancy forces, ${\mathbf{F}}_i^h,{\mathbf{T}}_i^h,{\mathbf{F}}_i^{me},{\mathbf{T}}_i^{me},{\mathbf{F}}_i^g$, that act on each particle’s surface are computed (step 2(c) in Algorithm 1). The hydrodynamic force acting on the surface of particle i, denoted by ${\mathbf{F}}_i^h$ and defined by Eq. (18), can be rewritten as,

$$\begin{array}{c}\hfill {\int }_{{\Omega }_s}\left(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\right)d{\Omega }_s={\int }_{\Omega }\left(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\right){\delta }_{\Omega }d\Omega ,\end{array}$$

(28)

where, $\mathbf{x}$ is an arbitrary region within the domain $\Omega$, and ${\delta }_{\Omega }=1$ if $\mathbf{x}\in {\Omega }_s$, otherwise ${\delta }_{\Omega }=0$. Assuming that $T_h$ is a decomposition of $\Omega$ consisting of computational cells c, we can approximate Eq. (28) as,

$$\begin{array}{c}\hfill {\int }_{\Omega }\left(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\right){\delta }_{\Omega }d\Omega =\sum _{c\in T_h}{\int }_{V\left(c\right)}\left(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\right){\delta }_{\Omega }dV\left(c\right),\end{array}$$

(29)

where $V\left(c\right)$ is the volume of cell c. Notice that for notation purposes we use the parentheses $\left(c\right)$ to evaluate a function on cell c. Numerical integration of Eq. (29) leads to the final form of the hydrodynamic forces acting on the particle,

$$\begin{array}{c}\hfill {\mathbf{F}}_i^h=\sum _{c\in {\overline{T}}_h}\left(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\right)\left(c\right)·V\left(c\right),\end{array}$$

(30)

where ${\overline{T}}_h$ is the set of all cells covered, in full or in part, by a magnetic particle. The resultant hydrodynamic torque on particle i, denoted by ${\mathbf{T}}_i^h$ and defined by Eq. (19), can be then approximated by taking the cross product between the position vector $\mathbf{r}$ and the total force from Eq. (30) that reads as,

$$\begin{array}{c}\hfill {\mathbf{T}}_i^h=\sum _{c\in {\overline{T}}_h}\left\{\mathbf{r}\left(c\right)\times \left(-\nabla p+{\eta }_S{\nabla }^2\mathbf{u}\right)\left(c\right)\right\}·V\left(c\right).\end{array}$$

(31)

Similarly, the magnetostatic polarization force and torque, defined by Eqs. (7) and (9), are approximated numerically as,

$$\begin{array}{c}\hfill {\mathbf{F}}_i^{me}\sum _{c\in {\overline{T}}_h}({\mu }_0{\chi }_e\mathbf{H}\nabla \mathbf{H})\left(c\right)·V\left(c\right),\end{array}$$

(32)

and

$$\begin{array}{c}\hfill {\mathbf{T}}_i^{me}=\sum _{c\in {\overline{T}}_h}({\mu }_0{\chi }_e\mathbf{H}\times \mathbf{H})\left(c\right)·V\left(c\right).\end{array}$$

(33)

The buoyancy force, defined by Eq. (20), can be also approximated numerically by integrating the fluid density $\left({\rho }_f\right)$ over the volume of the solid region in the mesh, i.e., $V\left(c\right)$ with $c\in {\overline{T}}_h$, to obtain the total displaced fluid mass, i.e., ${\rho }_{f}V\left(c\right)$ with $c\in {\overline{T}}_h$. Next, by multiplying the fluid mass by the gravitational acceleration vector $\left(\mathbf{g}\right)$, the buoyancy force can be calculated as,

$$\begin{array}{c}\hfill {\mathbf{F}}_i^g=\sum _{c\in {\overline{T}}_h}\left({\rho }_f\mathbf{g}\right)\left(c\right)·V\left(c\right).\end{array}$$

(34)

As the next step in the FVM-IBM-DEM-MAG algorithm, the resulting forces and torques for each particle are returned to the DEM solver. Additionally, if collision between particles or particle-wall are detected, the collision force and torque, ${\mathbf{F}}_{ij}^c$ and ${\mathbf{T}}_{ij}^c$, are calculated using Eqs. (23)–(27). Finally, the dipole-dipole magnetic force and torque, ${\mathbf{F}}_{ij}^{d-d}$ and ${\mathbf{T}}_{ij}^{d-d}$, are calculated using Eqs. (10) and (11) with either the fixed dipole model for dilute suspensions or the mutual dipole model for non-dilute suspensions to retrieve the particle dipole moment.

A data exchange model is also used to run a DEM script, which computes the particles’ positions, translational and angular velocities (Eqs. (21)–(22)), using Velocity-Verlet integration [39] (step 2(d) in Algorithm 1). The particles’ new positions and velocities are then transferred to the CFD solver. The CFD solver proceeds with the PISO (Pressure-Implicit with Splitting of Operators) algorithm [40] (step 2(e) in Algorithm 1), which solves the magneto-static potential equation, Eq. (6), and fluid flow governing equations, Eqs. (15)–(17). An intermediate velocity field $\widehat{\mathbf{u}}$ is first obtained by solving the momentum balance equations, Eq. (16), and then an intermediate pressure $\tilde{p}$ is obtained from the continuity equation, Eq. (15), which results in a Poisson equation for the pressure correction.

The next step is to correct the intermediate velocity field $\widehat{\mathbf{u}}$ in the particle region by imposing the rigid body velocity provided by the DEM calculation (step 2(f) in Algorithm 1). This correction is equivalent to adding a body force per unit volume defined as,

$$\begin{array}{c}\hfill \mathbf{f}=\rho \frac{\partial }{\partial t}(\tilde{\mathbf{u}}-\widehat{\mathbf{u}}),\end{array}$$

(35)

in the momentum balance equations, Eq. (16), to obtain a corrected velocity field $\tilde{\mathbf{u}}$. Here $\tilde{\mathbf{u}}={\mathbf{U}}_i^p+{\omega }_i\times \mathbf{r}$ is defined only for the cells within the solid body. The translational and angular velocities, ${\mathbf{U}}_i^p$ and ${\omega }_i$, respectively, were previously computed in step 2(d).

The previous step introduces a discontinuity in the velocity field at the interface, giving rise to a non-zero divergence in that location. Hence, the velocity field $\tilde{\mathbf{u}}$ and the pressure field $\tilde{p}$ need to be corrected (step 2(g) in Algorithm 1). For that purpose, $\tilde{\mathbf{u}}$ is projected onto a divergence-free velocity space, $\overline{\mathbf{u}}$, by using a scalar field $\psi$, as:

$$\begin{array}{c}\hfill \overline{\mathbf{u}}=\tilde{\mathbf{u}}-\nabla \psi ,\end{array}$$

(36)

where $\psi$ is obtained by solving the following Poisson equation,

$$\begin{array}{c}\hfill {\nabla }^2\psi =\nabla ·\tilde{\mathbf{u}}.\end{array}$$

(37)

Then $\overline{\mathbf{u}}$ is calculated by Eq. (36). The last step is equivalent to adding a pressure force $-\rho \frac{\nabla \psi }{\Delta t}$ in the momentum conservation equations, which requires the pressure field to be corrected by,

$$\begin{array}{c}\hfill p=\tilde{p}+\frac{\psi }{\Delta t}.\end{array}$$

(38)

This new FVM-IBM-DEM-MAG solver is implemented within the $CFDEMcoupling$ [41] framework.

4. Results and discussion

This section presents the validation of the proposed FVM-IBM-DEM-MAG solver against several benchmark case studies. The first case study is devoted to the sedimentation of two sphere’s in a rectangular duct containing a Newtonian fluid, mimicking the so-called drafting-kissing-tumbling (DKT) phenomenon. We start by turning off the external magnetic field to verify the solver’s capabilities to simulate the motion and interaction of the two settling spheres. Subsequently, in the second case study, we activate both the external magnetic field and the dipole-dipole force (and resultant torque) between the spheres for the DKT problem. This case study allows us to test the implementation of the magnetic force acting on the particles induced by the external magnetic field and the nearby magnetized particles. The magnetic potential gradient is applied in the vertical and horizontal directions to verify the ability of the algorithm to predict particle chaining in both directions. The third case study tests the robustness of the FVM-IBM-DEM-MAG solver by computing multi-particle chaining with 260 and 390 spheres whose centers are located in a 2D plane. Finally, the fourth case study describes the multi-particle chaining when particles are randomly distributed in a 3D domain.

4.1. DKT phenomenon under zero magnetic field

The objective of this test case is to simulate the motion and the interaction of two equal rigid spheres settling in a duct as shown in Figure 2. The spherical particles are placed vertically with a distance equal to four particle’s radius. The leading sphere (i.e., the one in below) is slightly off-centred to avoid the symmetric solution. In this case, we expect the simulations to reproduce the well-documented DKT phenomenon, which has been observed in laboratory experiments [42] and modeled through numerical simulations using different computational methods [17,43,44,45]. This benchmark is specifically selected to test the accuracy and effectiveness of the FVM-IBM-DEM-MAG algorithm, when the magnetic field is set to zero [31,37].

The computational domain is $\Omega =\left[0,1\right]\times \left[0,1\right]\times \left[0,4\right]$ cm³. The diameter of the spheres is $d=1/6$ cm. The initial positions of the spheres centers are $(0.5,0.5,3.5)$ and $(0.5,0.49,3.16)$, and the fluid and spheres are initially at rest. On the boundary of the channel, no-slip fluid velocity is imposed. The fluid density is ${\rho }_f=1$ g/cm³, the sphere’s density is ${\rho }_s=1.14$ g/cm³, and the fluid kinematic viscosity is $\nu =0.01$ cm²/s [30]. For the potential inter-particle contacts and particle-wall contacts, the coefficient of normal restitution, coefficient of friction, Poisson’s ratio and Young’s modulus are considered to be $0.97$, $0.10$, $0.45$, and $2\times {10}^9$ Pa, respectively.

The numerical experiments were performed using two hexahedral meshes with initial configuration M1: $40\times 40\times 160$ and M2: $60\times 60\times 240$ grid cells. In addition, dynamic mesh capability ($dynamicRefineFvMesh$) [46] is used to refine the mesh near the solid-fluid interface at each time-step. In this work, the maxRefinement parameter (a property of the dynamic mesh method defining the maximum number of layers of refinement that a cell can experience) is equal to two layers. The simulation time-step is set to $\Delta t={10}^{-4}$ s corresponding to an average Courant number of $0.1$. The total computational elapsed time for the simulations was 1h52m and 6h16m for M1 and M2, respectively, executed on a 3.00-GHz 48 cores Intel Xeon Gold 6248R CPU processor with 128 GB of RAM.

Figure 3 shows the z-component of the spheres centers, $z_i^p$, and the z-component of the spheres translation velocities, ${\left(U_z\right)}_i^p$, as function of time for calculations using M1 and M2 meshes. Additionally, the results obtained by Glowinski et al. [30] using two levels of mesh refinement, $h_{\Omega }=1/60$ and $h_{\Omega }=1/80$, are included for comparison purposes. Our results for M1 and M2 meshes obtained using the newly-developed algorithm (Algorithm 1) are in good agreement indicating that the results are mesh independent. As can be observed, the particle on top (following particle) is first carried by the wake generated by the particle on the bottom (leading particle) forcing to the so-called drafting phenomenon ($0\le t<0.14$ s). Then, the following particle velocity increases, the distance between the two particle’s centres decreases, and ultimately a contact forms between them forcing to the so-called kissing phenomenon ($0.14<t<0.35$ s). Since the vertical configuration is unstable and particles cannot stay attached [47], the particles start tumbling and are found side by side, which is known as tumbling phenomenon ($0.35<t<0.5$ s). Subsequently, the following particle passes ahead of the leading particle causing the deviation of the leading particle from the middle of the channel influenced by the fluid’s back-flows along the wall. Ultimately, the particle stagnate against the wall ($t\approx 0.65$ s) [48].

When comparing our results with the results computed by Glowinski et al. [30] with $h_{\Omega }=1/60$ and $h_{\Omega }=1/80$, it can be seen that they both predicted similar physical behaviors but with small discrepancy on timing. It must be noted that the kissing, drafting, and tumbling (DKT) benchmark case study is a non-smooth case involving several symmetry breaking. The exact agreement between different numerical algorithms after the kissing phenomenon is difficult to achieve, in part due to the lack of achieving mesh-independent results, or the use of different inter-particle contact models that influences the particles’ position drastically. To show the completeness of our solution, Figure 4 presents particles location and the contour distribution for the longitudinal fluid velocity, $u_z$ (cm/s), obtained at the middle plane $x=0.5$ cm for times $t=0.01,0.30,0.35,0.45,0.50$, and $0.65$ s obtained with M2. One can distinctly observe that the drafting ($t=0.3$ s), kissing ($t=0.35$ s), and tumbling ($t=0.45$ s) phenomena are indeed taking place. Next test-cases explore how the DKT benchmark case study changes when particles are magnetized under a constant magnetic field.

4.2. DKT phenomenon under magnetic field

This computational experiments examines the effect of the application of an external magnetic field on the DKT benchmark case study described in Section 4.1. We apply the external magnetic field both vertically (see Figure 5(a)) and horizontally (see Figure 5(b)), and explore how the magnetic field affects the sedimentation of the two magnetic spheres, i.e., the DKT phenomenon [17]. In both cases, the applied magnetic potential gradient field, $\nabla \varphi$, is set to 50 A/m. In addition, the fixed dipole model (see Eq. (12)) is employed to magnetize the particles with a relative susceptibility of $\chi =2000$ [49].

Figure 6 shows the z-component of the spheres centers, $z_i^p$, and of the z-component of the spheres translation velocities, ${\left(U_z\right)}_i^p$, as function of time for the calculations using mesh M2. In addition, the results of the DKT benchmark case study under no magnetic field (obtained in Section 4.1) are also shown for comparison purposes. When the external magnetic field is applied in the vertical direction, the two magnetic particles are attracted together forming a string ($t\approx 0.5$s). The string last until they contact the bottom wall of the domain. On the other turn, when the external magnetic field is induced in the horizontal direction, the spherical particles do not approach each other, but instead they tumble side-by-side. During the rest of the sedimentation process, the wake generated by the leading particle leads to a faster settling of the following particle ($t>0.5$s), also see Figure 8 for an illustrative representation of this phenomenon.

Figure 7 presents particles settling under vertical magnetic field. The particle location and the contour distribution for the longitudinal fluid velocity, $u_z$(cm/s), obtained at the midplane $x=0.5$ cm for times $t=0.01,0.30,0.35,0.45,0.50$, and $0.65$ s are shown. As can be seen, the particles experience longer drafting period, and form a tight string that does not get separated in the rest of the sedimentation process.

Figure 8 shows the settling of the spherical particles under an horizontal magnetic field. The particle location and the contour distribution for the longitudinal fluid velocity, $u_z$ (cm/s) obtained at the midplane $x=0.5$ cm for times $t=0.01,0.30,0.35,0.45,0.50$, and $0.65$ s are shown in Figure 8. In this case, the direction of the particles sedimentation is transverse to the magnetic field direction, and hence, particles form a repulsive magnetic force [17]. For that reason, the particles, instead of approaching and contacting each other, just tumble as a non-kissing pair ($t\approx 0.50$ s). Shortly after the tumble, the particles approach the vertical walls, where the external magnetic field is applied ($t\approx 0.65$ s).

Figure 8. The change in drafting-kissing-tumbling (DKT) benchmark case study under horizontal magnetic field simulated using Algorithm 1. The positions of spheres at $t=0.01,0.30,0.35,0.45,0.50$ and $0.65$ s, and the contour of the longitudinal ($z-$component) fluid velocity, $u_z$ (cm/s), at the midplane $x=0.5$ cm are shown.

Figure 8. The change in drafting-kissing-tumbling (DKT) benchmark case study under horizontal magnetic field simulated using Algorithm 1. The positions of spheres at $t=0.01,0.30,0.35,0.45,0.50$ and $0.65$ s, and the contour of the longitudinal ($z-$component) fluid velocity, $u_z$ (cm/s), at the midplane $x=0.5$ cm are shown.

4.3. Multi-particle chaining under magnetic field: 2D

In this test case, we analyze the motion of a random array of magnetic particles whose centers are located at the midplane of a rectangular box filled with a Newtonian fluid under the influence of external magnetic fields [15,17,49]. Two computational domains are employed as ${\Omega }_h=\left[0,4\right]\times \left[0,1\right]\times \left[0,1\right]$ cm³ and ${\Omega }_v=\left[0,1\right]\times \left[0,4\right]\times \left[0,1\right]$ cm³ (see Figure 9). The initial positions of the spheres centers, with diameter $d=1/16$ cm and density of ${\rho }_s=1.01$ g/cm³, are randomly generated and constrained such that the minimum distance between particles and between the particles and walls is equal to $1.5d$. The spheres move under the action of gravity, hydrodynamic forces, mutual dipole-dipole forces, and the applied external magnetic force [17].

Two area fractions of spheres were tested, $20\%$ and $30\%$, corresponding to 260 and 390 spheres under the effect of both a vertical and horizontal magnetic fields with magnetic potential gradient of $\nabla \varphi =50$ A/m. The fluid and the spheres are initially at rest. On the channel walls, the no-slip boundary condition is applied for the fluid velocity. A cyclic boundary condition is applied on the other boundaries. In addition, the fluid density and kinematic viscosity are set to ${\rho }_f=1$ g/cm³ and $\nu =0.01$ cm³/s, respectively. The dipole-dipole magnetic forces and torques are calculated using the mutual dipole model, see Eq. (13), with a relative susceptibility of $\chi =2000$[49]. For the inter-particle contacts and particle-wall contacts, the coefficient of normal restitution, coefficient of friction, Poisson’s ratio and Young’s modulus are considered to be $0.90$, $0.33$, $0.33$, and $7\times {10}^8$ Pa, respectively [17].

The calculations were performed in an hexahedral mesh with initial configuration $128\times 32\times 32$ grid cells for the horizontal domain (M_h) and $32\times 128\times 32$ grid cells for the vertical domain (M_v). Again, the dynamic mesh refinement was employed in the calculations with two levels of refinement. The time-step used in the simulations is $\Delta t={10}^{-4}$ s corresponding to a maximum Courant number of $0.1$. The total computational elapsed time for the simulations was approximately 2h05m and 2h35m for the $20\%$ and $30\%$ particle’s area fractions executed on a 3.00-GHz 48 cores Intel Xeon Gold 6248R CPU processor with 128 GB of RAM.

Figure 10 and Figure 11 show the snapshots of 260 and 390 particles moving in the rectangular channel under the effect of gravity and an external magnetic field applied in the vertical direction. At the initial instants of the simulations ($t\le 0.2$ s), short fragmented chains or clusters of particles are formed in the $y-$direction (the same as the applied external magnetic field direction). Subsequently, at later instants ($0.4\le t\le 0.6$ s), the short chains start to merge together and form long chains, i.e., they form mesoscale structures made of magnetic particles with shapes and orientations comparable to the results presented by Han et al. [15], Ke et al. [17] and Ly et al. [49].

Figure 12 and Figure 13 show the snapshots of 260 and 390 particles moving in the rectangular channel under the action of gravity and of an external magnetic field applied in the horizontal direction. Again, at the initial instants of the simulations ($t\le 0.2$ s), short fragmented chains or clusters of particles are formed in the $x-$direction (the same as the applied external magnetic field direction). Then, at later instants ($0.4\le t\le 0.6$ s), the short chains start to merge together and form long horizontally aligned chains, i.e., they form mesoscale structures made of magnetic particles with distinct shapes and orientations.

Figure 10 to Figure 13 also show the presence of isolated magnetic particles and a number of chains with shorter lengths. Predominantly, these chains are linear as head-to-tail aggregation of magnetic dipoles, but as claimed by Ke et al. [17], Mohebi et al. [50] and Fermigier and Gast [51], it is also observed thick particle clusters due to the lateral merging of the linear chains.

4.4. Multi-particle chaining under magnetic field: 3D

In this subsection, we analyze the robustness of the proposed FVM-IBM-DEM-MAG solver by studying the chain formation in MRFs within a three-dimensional (3D) domain. A random array of magnetic spheres is placed in a rectangular box filled with a Newtonian fluid under the influence of gravity and external magnetic fields applied in different directions [15]. The computational domain employed was $\Omega =\left[0,2\right]\times \left[0,2\right]\times \left[0,1\right]$ cm³ (see Figure 14). The diameter of the spheres is $d=1/16$ cm. The initial positions of the spheres centers are randomly generated with a restriction such that the minimum distance between particles and between the particles and walls is higher than $1.5d$. The spheres move under the action of gravity, hydrodynamic forces, mutual dipole-dipole forces, and the applied external magnetic force. The sphere volume fraction was fixed at $1.85\%$, corresponding to 580 spheres. We considered an external magnetic field with magnetic gradient potential $\nabla \varphi =50$ A/m applied vertically or horizontally. The fluid and the spheres are initially at rest. On the channel walls, the no-slip boundary condition is imposed for the fluid velocity. A cyclic boundary condition is applied on the other boundaries. The fluid and particle densities are ${\rho }_f=1$ g/cm³ and ${\rho }_s=1.01$ g/cm³, respectively. The fluid kinematic viscosity is $\nu =0.01$ cm²/s. The dipole-dipole magnetic forces and torques are calculated using the mutual dipole model, see Eq. (13), with a relative susceptibility of $\chi =2000$[49]. For the inter-particle contacts and particle-wall contacts, the coefficient of normal restitution, coefficient of friction, Poisson’s ratio and Young’s modulus are considered to be $0.90$, $0.33$, $0.33$, and $7\times {10}^8$ Pa, respectively [17].

The calculations were performed in an hexahedral mesh with initial configuration of $64\times 64\times 32$ grid cells. Again, the dynamic mesh refinement was employed in the calculations with maxRefinement = 2. The time-step used in the simulations is $\Delta t={10}^{-4}$ s, corresponding to a maximum Courant number of $0.1$. The total computational elapsed time for the simulations was approximately 18h12m executed on a 3.00-GHz 48 cores Intel Xeon Gold 6248R CPU processor with 128 GB of RAM.

Figure 15 and Figure 16 depict the evolution of the particles at six time instants for the two directions of the imposed external magnetic potential field. It can be seen that with the application of the magnetic field, the particles become magnetized and acquire a magnetic dipole moment [15], which promotes the particles to aggregate and form short fragmented chains ($t\le 1$ s). As time advances, these short chains merge together and form longer chains (i.e., mesoscopic structures) that align in the direction of the applied magnetic field [15].

5. Conclusions

A numerical formulation for fully-resolved simulation of magnetorheological fluids (MRF) consisting of solid magnetic particles suspended in a Newtonian carrier fluid was presented. The implementation was carried out by extending the open-source $CFDEMcoupling$ framework with a force calculation at the particles surface due to the applied external magnetic field, and with the implementation of the fixed and mutual dipole-dipole magnetic models to account for the magnetic interactions between the particles. The overall algorithm procedure solves a second-order differential equation for the magnetic potential field, followed by the flow equations, including the continuity and momentum balance equations, and an immersed boundary algorithm to model the flow around discrete magnetic particles present in the flow domain. This approach guarantees a tight coupling between the dynamics of the fluid and the magnetic solid discrete phase. The coupling is provided by the calculation of the net hydrodynamic and magnetic forces and torques exerted by the fluid on the solid particles. The algorithm subsequently uses the discrete element method to model the particle motion, comprising linear and rotational motions, as well as the particles magnetic moment, which in turn provides new boundary conditions for the fluid domain.

The accuracy and robustness of the proposed FVM-IBM-DEM-MAG algorithm were evaluated using four benchmark studies. First, for the sedimentation of two spheres in a rectangular duct containing a Newtonian fluid without the presence of an external magnetic field (mimicking the drafting-kissing-tumbling, DKT, phenomena), the particles velocity and location were compared with numerical data available in the literature and a good agreement was obtained. The velocity contour profiles of the particles falling through the Newtonian fluid distinctly showed several symmetry breaking physical aspects of the non-smooth DKT phenomenon. We also demonstrated the capability of the algorithm to predict the dynamics of two magnetic particles falling under the action of gravity and an external magnetic field, i.e., the simulation of the DKT benchmark case study but with activating the magnetic forces calculations. For a vertical magnetic field, the particles experience a longer drafting period and form a tight string which does not separate during the rest of the sedimentation process. For a horizontal magnetic field, the particles just tumble as a non-kissing pair and approach the vertical walls of the domain, where the external magnetic field is applied. The FVM-IBM-DEM-MAG solver was also used to study the multi-particle chaining when particles are placed randomly on a 2D-plane. Two area fractions of spheres were tested, $20\%$ and $30\%$, corresponding to 260 and 390 spheres under the effect of gravity and a vertical or horizontal magnetic fields. The snapshots of the particles locations showed that, at the initial instants, short fragmented chains or clusters of particles are formed. With time advancing, the short chains merge together and form longer column-like chains always aligned with the direction of the externally imposed magnetic field. Finally, the robustness of the FVM-IBM-DEM-MAG solver was tested in a 3D domain, where an array of 580 randomly distributed magnetic particles were subjected to gravity and a horizontal or vertical magnetic field. Again, the snapshots of the particles location demonstrated the formation of long column-like chains in the direction of the applied magnetic field.

In summary, the results presented in this study show that the newly developed code can accurately predict the flow patterns and particle assembly in MRF for a number of benchmark problems.