Ferrofluids and magnetically guided superparamagnetic particles in flows: A review of simulations and modeling

TL;DR

This review discusses the mathematical modeling of ferrofluids and superparamagnetic particles in flows, highlighting their behavior, applications, and recent advances in simulations relevant to biomedical and engineering fields.

Contribution

It provides a comprehensive overview of modeling approaches for superparamagnetic particles and ferrofluid flows, emphasizing recent developments and applications.

Findings

Models for nanoparticle motion in drug targeting

Analysis of ferrofluid drop dynamics

Insights into ferrofluid film breakage

Abstract

Ferrofluids are typically suspensions of magnetite nanoparticles, and behave as a homogeneous continuum. The production of nanoparticles with a narrow size distribution and the achievement of colloidal stability are important technological issues. The ability of the ferrofluid to respond to an external magnetic field in a controllable manner has made it emerge as a smart material in a variety of applications, such as seals, lubricants, electronics cooling, shock absorbers and adaptive optics. Magnetic nanoparticle suspensions have also gained attraction recently in a range of biomedical applications, such as cell separation, hyperthermia, MRI, drug targeting and cancer diagnosis. In this review, we provide an introduction to mathematical modeling of three problems: motion of superparamagnetic nanoparticles in magnetic drug targeting, the motion of a ferrofluid drop consisting of chemically bound nanoparticles without a carrier fluid, and the breakage of a thin film of a ferrofluid.

Equations

$$\rho V\frac{d^2\mathbf{x}}{d{t}^2}=-D(\frac{d\mathbf{x}}{dt}-{\mathbf{u}}_{\mathbf{b}})+\frac{3\chi V{\mu }_0}{3+\chi }\nabla \left(\frac{1}{2}|\nabla \varphi {|}^2\right)+\frac{\sqrt{2DkT}}{\sqrt{\mathbf{d}t}}\mathbf{N}(0,1),$$

$$dX$$

$$f(X)$$

$$g$$

$$\frac{\partial \rho }{\partial t}+\mathbf{u}\cdot \nabla \rho$$

$$\rho \frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}$$

$${\mathbf{\Pi }}_M=-\frac{{\mu }_0}{2}(\mathbf{H}\cdot \mathbf{H})\mathbf{I}+\mu \mathbf{H}\mathbf{H},$$

$$\nabla \cdot \mathbf{B}$$

$$\mathbf{B}(\mathbf{x},t)$$

$$\nabla \cdot (\mu \nabla \psi )=0.$$

$$\mathbf{M}(\mathbf{H})=M_{s}L\left(\frac{{\mu }_{0}m|\mathbf{H}|}{k_{B}T}\right)\frac{\mathbf{H}}{|\mathbf{H}|},$$

$$[[\mathbf{u}]]=0,\mathbf{n}\cdot [[\mathbf{B}]]=0,\mathbf{n}\times [[\mathbf{H}]]=0,[[{{\mathbf{t}}_1}^T\cdot \mathbf{\Sigma }\cdot \mathbf{n}]]=0,[[{{\mathbf{t}}_2}^T\cdot \mathbf{\Sigma }\cdot \mathbf{n}]]=0.$$

$$[[{\mathbf{n}}^T\cdot \mathbf{\Sigma }\cdot \mathbf{n}]]=\sigma k,$$

$$C(r,z,t)$$

$$\frac{\partial C}{\partial t}+\mathbf{u}\cdot \nabla C$$

$$\rho \frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}$$

$${\mathbf{F}}_s=\sigma k{\delta }_s\mathbf{n},$$

$$\frac{1}{r}\frac{\partial }{\partial r}\left(\mu r\frac{\partial \psi }{\partial r}\right)+\frac{\partial }{\partial z}\left(\mu \frac{\partial \psi }{\partial z}\right)=0.$$

$$\tilde{\text{x}}=\frac{\text{x}}{L_c}$$

$$\tilde{\text{u}}=\frac{\tau }{L_c}\text{u}$$

$${\tilde{\mathbf{F}}}_s=\frac{{\mathbf{F}}_s}{\sigma /L_c^2}$$

$$\tilde{p}=\frac{p}{\sigma /L_c}$$

$$\tau =\frac{{\eta }_d{L}_c}{\sigma },$$

$$L_c=\frac{{\eta }_d^2}{{\rho }_d\sigma },$$

$$\text{mbox}{Oh}^{-2}\left(\rho \frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)$$

$$\mathbf{H}=(\frac{\partial \psi }{\partial r},\frac{\partial \psi }{\partial z},0),$$

$${\mathbf{\Pi}}_{M}=\mu\left[\begin{array}[]{ccc}\psi_{r}^{2}-\frac{1}{2}\left(\psi_{r}^{2}+\psi_{z}^{2}\right)&\psi_{r}\psi_{z}&0\\ \psi_{r}\psi_{z}&\psi_{z}^{2}-\frac{1}{2}\left(\psi_{r}^{2}+\psi_{z}^{2}\right)&0\\ 0&0&-\frac{1}{2}\left(\psi_{r}^{2}+\psi_{z}^{2}\right)\end{array}\right].$$

$$\nabla \cdot {\mathbf{\Pi }}_M=[\frac{1}{r}\frac{\partial }{\partial r}[r({\Pi }_M{)}_{rr}]$$

$$\nabla \cdot (\mu \nabla \psi {)}_{i,j}$$

$$(\frac{\partial \psi }{\partial r}{)}_{i+1/2,j}=\frac{{\psi }_{i+1,j}-{\psi }_{i,j}}{\Delta }.$$

$$\frac{1}{r_{i+1/2,j}}\frac{r_{i+1,j}(({\Pi }_M{)}_{rr}{)}_{i+1,j}-r_{i,j}(({\Pi }_M{)}_{rr}{)}_{i,j}}{\Delta }+$$

$$\frac{(({\Pi }_M{)}_{rz}{)}_{i+1/2,j+1/2}-(({\Pi }_M{)}_{rz}{)}_{i+1/2,j-1/2}}{\Delta },\mathrm{and}$$

$$\frac{1}{r_{i,j+1/2}}\frac{r_{i+1/2,j+1/2}(({\Pi }_M{)}_{rz}{)}_{i+1/2,j+1/2}-r_{i-1/2,j+1/2}(({\Pi }_M{)}_{rz}{)}_{i-1/2,j+1/2}}{\Delta }+$$

$$\frac{(({\Pi }_M{)}_{zz}{)}_{i,j+1}-(({\Pi }_M{)}_{zz}{)}_{i,j}}{\Delta },$$

$$(({\Pi }_M{)}_{rr}{)}_{i,j}$$

$$(({\Pi }_M{)}_{rz}{)}_{i+1/2,j+1/2}={\mu }_{i+1/2,j+1/2}\left(\frac{{\psi }_{i,j}-{\psi }_{i-1,j}+{\psi }_{i,j-1}-{\psi }_{i-1,j-1}}{2\Delta }\right)\times$$

Full text

††journal: Journal of Engineering Mathematics

1 Introduction

This review presents three aspects of the mathematical analysis of the motion of superparamagnetic particles and drops in a surrounding fluid, under the influence of an external magnetic field: drug targeting, ferrofluid drop deformation, and instabilities and dewetting of a thin ferrofluid film. These problems are important for a number of biomedical applications. The motion of a magnetic particle is determined by a force balance of the magnetic force, drag, and stochastic forces. The motion of a ferrofluid drop is affected by its deformation. Simulations using the Volume-of-Fluidk method are presented. Once a drop reaches the target and coats it, a thin film approximation yields insight into instabilities and breakup. We refer the reader to recent publications, which in turn contain more comprehensive references and historical reviews.

An advantage of magnetic drug targeting is the delivery of drugs to specific tissues via the blood vessels as pathways under externally applied magnetic fields. Progress towards robust biomedical applications relies on the synthesis of appropriately coated nanoparticles, clusters and ferrofluids Liu2007 ; ATRRWPR ; Bala2014 ; Dung2017 . Pros and cons of directing the magnetic nanoparticles in blood vessels and tissues are described in Voltairas02 ; Neuberger2005 ; Buzea2007 ; Berry2009 ; Mishra2010 ; Nacev2011 . Recent reviews highlight the use of magnetic nanoparticles in many biomedical applications; for instance, drug targeting, magnetic fluid hyperthermia, design of devices to measure targeting efficiency, and tissue engineering Ito2005 ; Puri2014 ; Maguire2014 ; aljamal2016 ; Mohammed2017 ; Radon2017 .

In Suh2011 , the mathematical modeling of the transport of paramagnetic particles in viscous flow is described in terms of dipole points in an external magnetic field, with particle interactions for larger particles and higher concentrationsBanerjee2012 . Numerical approaches using Brownian dynamics for particle interactions in aggregation and disaggregation are being developed vanReenen2014 . Technological applications include the sorting of not only cells but proteins and other biological components by biochemically funcionalized paramagnetic beads that are manufactured to be attracted to specific targets. The mixture flows in a fluid, and as they move through a channel, magnetic fields deflect the desired cells to migrate toward predetermined exits. In Tsai2011 , paramagnetic beads floating in a liquid are subjected to a magnetic field normal to the flow. The deflection by the magnetic field steers the beads towards a target whose location depends on the field and on the properties of the bead. An application to removing pathogens from the blood stream is shown in small , reproduced in figure 1. Here the magnetic particles bind to pathogens such as bacteria. In addition to the motion of the particles, a model for the kinetics of this binding mechanism is formulated.

The design of the superparamagnetic nanoparticles involves uniformity in the small size. Of particular interest is the formation of a ferrofluid from chemically bonding many of them. Unlike the colloidal suspension of magnetite particles, this ferrofluid avoids the migration of the particles in the suspending fluid at high fields. This advantage is seen in the works of Mefford ; Mefford2008a ; Mefford2008b for the potential treatment of retinal detachment by magnetically guiding a ferrofluid drop to seal the site. In particular, figure 2(a) shows a sample transmission electron microscopy (TEM) image of the magnetite nanoparticles at $300000$ times magnification, and (b) shows the narrow size distribution. These particles are relevant to section 2. The magnetic nanoparticles are coated with biocompatible polydimethylsiloxane (PDMS) oligomers, and bonded to manufacture a ferrofluid. Figures 2(c-d) show experimentally observed states for the ferrofluid drop suspended in glycerol in uniform magnetic fields. The modeling in section 3.2 concerns the numerical simulation of stable equilibrium configurations under applied uniform magnetic fields. A series of experiments conducted under low and high magnetic field strengths are presented in ATRRWPR to understand the behavior of PDMS ferrofluid drops in uniform magnetic fields. Figure 3 shows the experiments in ATRRWPR of the drop at equilibrium for the magnetic field strengths varying from $6.032$ to $162.33$ kA m*-1*. These results are revisited in section 3.2. For an updated review of literature on ferrofluid drop deformation in uniform magnetic fields, see Rowghanian .

The numerical investigation of the motion of a ferrofluid drop may be complicated by the highly distorted material interface. An advantage of the boundary integral method for three-dimensional (3D) drop deformation in Stokes flow is that the three-dimensionality is turned into a two-dimensional surface computation; this together with a high order surface discretization has been able to track viscous drops with high curvature, up to the first pinch-off CristiniBLSG ; Janssen2008 . However, this method does not extend easily to non-Newtonian models such as a power-law model that is suggested for blood flow. Novel nonsingular boundary integral methods are improving the accuracy and numerical stability of the 3D drop deformation in viscous flow at low inertia, including the influence of a magnetic field for constant magnetic susceptibility Bazhlekov06 . The effect of rotating magnetic field is investigated with the boundary integral formulation in Erdmanis2017 . The study of rotating magnetic fields on ferrofluid drops is of importance in many fields including astrophysics, and we refer the reader to recent reviews Lebedev2003 ; Fengchen2016 . A limitation of the boundary integral formulation is that it is no longer applicable if the magnetic behavior is nonlinear.

Numerical algorithms with a diffuse-interface method have been developed for a droplet or layer in a uniform and rotating magnetic fields Fengchen2016 . Intricate equilibrium shapes arise in the Rosensweig instability which originally refers to a ferrofluid layer under magnetic force, surface tension force, gravity and hydrodynamics Cowley ; Rosensweig ; Lange2007 ; Kadau2016 . In Lavrova2012 , the Rosensweig instability is investigated numerically for the diffusion of magnetic particles. The magnetic particle concentration and the free surface shape are part of the solution. A finite element method is used for Maxwell’s equations and a finite difference method is used for a parametric representation of the free surface. A finite element formulation is developed and used in Kang2013 to simulate the interaction of two or more particles in uniform and rotating magnetic fields. The interaction is significant in dense suspensions or larger sizes. In section 3.2, we focus on a Volume-of-Fluid formulation of ARRRP which is capable of utilizing a non-constant susceptibility, time-dependent evolution, and interface breakup and reconnection.

Finally, thin films driven by a magnetic field have been studied recently Seric2014 ; Conroy2015 . These models are derived for the flow of a thin ferrofluid film on a substrate in the same spirit as the previous work by Craster and Matar Craster2005 for the case of a leaky dielectric model. The derivation for a nonconducting ferrofluid simplifies since there is no interfacial charge involved. The work of Seric et al. Seric2014 includes the van der Waals force to impose the contact angle using the disjoining/conjoining pressure model and studies the film break up and the formation of satellite (secondary) droplets on the substrate under an external magnetic field. In section 3.5, we focus on a thin film approximation of ferrofluid on a substrate subjected to an applied uniform magnetic field and present the results of the magnetic field induced dewetting of thin ferrofluid films and the consequent appearance of satellite droplets.

2 Mathematical modeling of magnetic drug targeting

In this section, we summarize the set-up and results of YLAR2012 . One of the simplest model problems for predicting the accuracy of drug targeting concerns superparamagnetic nanoparticles of at most several nanometers diameter, and flowing through a venule of order 10*-4* m radius; see figures 4 (a - b). The velocity is $(u(y,z,t),0,0)$ where $u$ denotes the axial component. The no-slip conditions at the walls $y^{2}+z^{2}=R^{2}$ and the Stokes equation are applied. For such a small radius, the effect of cardiac pumping on pressure fluctuations is negligible; therefore, $u=u_{max}(1-\frac{r^{2}}{R^{2}})$, where the flow rate $\frac{\pi}{2}u_{max}R^{2}$ is typically $10^{-10}$m3s*-1* House86 . This gives $u_{max}\approx 0.006$m s*-1*. This flow generates a Stokes drag force on a solid sphere: ${\bf F}_{v}=-(6\pi\eta a)(\frac{d{\bf x}}{dt}-(u,0,0))$ where $\eta$ is the viscosity, $a$ the radius, ${\bf x}$ is the location of the particle, and $D$ is the friction coefficient. For a larger blood vessel, the oscillatory component of the flow can be estimated from the systolic and diastolic blood pressure data.

The external magnetic field ${\bf H}_{e}$ without a particle is known from experimental measurements and is used to calculate the field in the presence of the particle. Under the assumptions of constant susceptibility $\chi$ (magnetic field not too high), the permeability is $\mu=\mu_{0}(1+\chi)$, where $\mu_{0}$ denotes the permeability of vacuum. Maxwell’s equations (curl ${\bf H}_{e}={\bf 0}$, div ${\bf H}_{e}=0$) give ${\bf H}_{e}=\nabla\phi$, where $\phi$ is a harmonic function. The external magnet is modeled by a magnetic dipole. We find the actual field ${\bf H}$ with the particle by using the known formula $\phi=\frac{-1}{4\pi}\frac{{\bf m}\cdot{\bf r}}{r^{3}}$, where the vector ${\bf r}$ points from the dipole to the particle. Provided that the external magnet is small, it is approximately a point dipole. This magnetic field exerts a force on the particle: ${\bf F}_{m}=\int{\mu_{0}({\bf M}\cdot\nabla){\bf H}_{e}dV}$ where ${\bf M}=\chi{\bf H}$ denotes the magnetization corresponding to the magnetic field ${\bf H}$ with a particle of volume $V=\frac{4}{3}\pi a^{3}$, namely the Clausius-Mossotti formula Rosensweig ; Cohen2007  ${\bf M}=\frac{3\chi}{3+\chi}{\bf H}_{e}$. This is used to calculate ${\bf F}_{m}$. The governing equation must satisfy ${\bf F}_{m}+{\bf F}_{v}={\bf 0}$. This is a system of nonlinear ordinary differential equations.

The superparamagnetic particles and clusters addressed in YLAR2012 are small enough that the Néel relaxation time and the magnetization are instantaneous compared with the motion to be solvedBala2014 . However, we need to keep in mind the that Brownian motion should be included if thermal energy $kT$ becomes comparable to magnetic energy $\mu_{0}MHV$, where $k$ denotes Boltzmann’s constant, $T$ is in degrees Kelvin, $H$ is the field and $V$ is the volume of the particle. This estimate gives the diameter to be smaller than $(6kT/\pi\mu_{0}MH)^{1/3}$ ((2.2) of Rosensweig ). In the examples below, this is diameter of order $10^{-8}$m. We incorporate Langevin’s model for Brownian motion through a stochastic forcing vector. This leads to the system

[TABLE]

where $\rho$ is the density of the particle/cluster, and ${\bf u_{b}}$ is the base flow. The left hand side term denotes acceleration which is estimated to be relatively small. $N_{1}(0,1)$, $N_{2}(0,1)$ and $N_{3}(0,1)$ denote independently generated, normally distributed, random variables with zero mean and unit variance . The random variables $N_{i}(0,1)$, $i=1,2,3,$ are constants over very short time intervals $dt$, and change randomly with a Gaussian distribution. Finally, the stochastic ordinary differential equation for the state vector $X(t)={\vec{x}}^{T}=(x(t),y(t),z(t))^{T}$ is,

[TABLE]

$\mathbf{I}$ is the $3\times 3$ identity matrix, $f(X(t))$ is the $3\times 1$ drift-rate vector, $g$ is a $3\times 3$ instantaneous diffusion-rate matrix, and $dW(t)$ denotes the vector $\sqrt{{\bf d}t}{\bf N}(0,1)$.

2.1 Numerical algorithm

The system (4) is integrated using the Euler-Maruyama method. At the $n$th time step, where $\Delta t$ denotes a fixed step size, $X_{n+1}=X_{n}+\Delta t\ f(X_{n})+g(X_{n})\Delta W_{n}$. As detailed in Higham2001 , the advantage of this method is the strong convergence under broad assumptions. We use the implementation in the SDE Toolbox for Matlab Picchini . The capture rate for a particle is defined to be the probability of hitting the tumor, which is situated at $|x|<R_{tumor}$ at the wall in figure 4. Figure 5(a) shows sample trajectories, initially at $x=-R_{tumor}$, where $R_{tumor}$ denotes the radius of the tumor at the wall. A trajectory can escape without hitting the wall ($-.-$), or escape after hitting the wall ($-$), or be captured at the target. Figures 5(b-d) use data from ARRRP ; House86 . The magnet is $d=0.03$ m outside the vessel wall; vessel radius $R=10^{-4}$ m, $u_{max}=\frac{2}{\pi}0.01$ ms*-1*, particle radius $a=10^{-7}$ m, $\eta=0.004$Pa s, $R_{tumor}=100R$, constant susceptibility $\chi=0.2$. (b) shows the capture rate versus dipole moment $m$. For larger $m$, a non-constant Langevin fit approximates the magnetization data ARRRP . The capture rate scales with the square of the dipole moment even when the velocity profile is slightly blunted at the centerline which may relate to blood flow. To optimize the capture rate, it is found that the distance $d$ should match the tumor size. (c) shows trajectories for $a=10^{-7}$ m, and compared with this, (d) shows that Brownian motion makes a significant difference for $a=O(10^{-8})$ m.

The capture rate is optimal for a magnet that produces a strong field and a focused effect, and therefore high field gradients. This is one of the technological challenges in magnetic drug targeting. Recent studies show that arrays of magnets, such as the Halbach array, are an improvement over the single magnet. The balance of magnetic force and hydrodynamic drag for Stokes flow is calculated for high field gradient arrays and compared with experimental data in Barnsley2017 . This gives guidance on how the particle trajectories depend on magnet shapes and arrays.

Non-Newtonian effects in blood flow have been studied with empirical formulas that replace the Stokes drag coefficient Sprenger2015 . For example, in the case of the larger arterial speeds, an empirical viscosity term which incorporates shear-thinning at the local shear rate is applied to write approximate equations for the motion of micron sized magnetic particles in Poiseuille flow Cherry2014 ; Mohren2017 . The particles in Cherry2014 are of micron size, much larger than the nanometers considered in YLAR2012 , and also the diameter of the artery and flow speeds are much larger, and therefore the Brownian effect for the particles is small. The use of the local shear rate of the background flow in the Stokes drag formula is an ad-hoc substitution. For a particle moving in non-Newtonian flow, such as in a shear-thinning fluid, the calculation of drag is tractable for quiescent flow, not for non-quiescent flow such as shear flow. Unlike Stokes flow, the equations are nonlinear so that solutions can not be superposed. The real velocity of the particle is spatially and temporally varying in response to the magnetic field and flow. A separate effect is the influence of red blood cells on the capture rate. Mohren2017 applies a formula for shear-enhanced diffusion in colloidal suspensions (equation (1) of Griffiths2012 ) to estimate that the red blood cells contribute an additive term proportional to the shear rate inside the square root in equation 4. Particle trajectories for the Cherry viscosity model and including the diffusion due to interaction with red blood cells are investigated in Mohren2017 . They find that the additional diffusion deteriorates the capture rate; figure 6 shows this.

3 Mathematical modeling of ferrofluid drops

To investigate the response of a ferrofluid droplet to an applied magnetic field, extra effort is needed to solve the details of the magnetic field, which in turn requires the knowledge of the shape of the interface. The non-local moving boundary problem which couples the spatial variation of the magnetic field with the interfacial shape of the drop necessitates a direct numerical computation. Numerical algorithms that explicitly track the interface and use finite element methods are developed to investigate equilibrium shapes of ferrofluid drops and interfacial instabilities in Lavrova ; Lavrova04 ; Bashtovoi ; Matthies ; Knieling . For the case of constant magnetic susceptibility, a boundary-integral formulation is used to study equilibrium shapes for high-frequency rotating magnetic fields Erdmanis2017 . Numerical algorithms that implicitly track the interface include the Volume-of-Fluid (VOF) methods, Level-Set methods and diffuse-interface (DI) methods Scar99 ; Sethian2003 ; AMW . Each of these methods can be applied to subcases of the results we show below. We specifically describe the VOF method because of its simplicity, efficiency and robustness for tracking topologically complex evolving interfaces. However, this as well as other approaches has its limitations. An advantage is that it conserves mass and with the recent improvements in the discretization of the surface tension force, it remains a competitive method for modeling interfacial flows. We present these results in section 3.2.

A recent study has shown that a ferrofluid droplet on a hydrophobic substrate breaks up into two daughter droplets and a smaller satellite droplet under a sufficiently strong magnetic field Timonen2013 . These experiments also show that if the magnetic field is increased even further, the droplets break up again and rearrange to form an assembly of droplets on the substrate. In the same spirit, we present a model derived in Seric2014 for the flow of a thin ferrofluid film on a substrate. The results in Seric2014 show the film break up and the deformation of sessile droplets on the substrate under an external magnetic field. We present this study in section 3.5

3.1 Governing equations

The equations governing the motion of an incompressible ferrofluid drop suspended in another fluid are the conservation of mass and the balance of momentum

[TABLE]

where ${\mathbf{u}}({\mathbf{x}},t)$ is the velocity field, $p({\mathbf{x}},t)$ the pressure, $\rho({\mathbf{x}},t)$ and $\eta({\mathbf{x}},t)$ are the phase dependent density and viscosity, respectively, {\mathbf{D}}=\frac{1}{2}\big{(}\nabla{\mathbf{u}}+(\nabla{\mathbf{u}})^{T}\big{)} is the rate of deformation tensor (where $T$ denotes the transpose), $\rho\mathbf{g}=-\rho g{\hat{z}}$ the body force due to gravity and ${\hat{z}}$ the unit vector in the $z$-direction (with gravitational constant $g$), and ${\mathbf{\Pi}}_{M}({\mathbf{x}},t)$ is the magnetic stress tensor

[TABLE]

where ${\mathbf{I}}$ is the identity tensor. The Maxwell equations for the magnetic induction ${\mathbf{B}}$, magnetic field ${\mathbf{H}}$, and magnetization ${\mathbf{M}}$ are the magnetostatic Maxwell equations for a non-conducting ferrofluid

[TABLE]

where $\mu_{d}$ denotes the permeability of the drop, and $\mu_{m}$ is the permeability of the surrounding fluid. Here, we assume that the surrounding fluid has a permeability very close to that for a vacuum, $\mu_{0}$. Therefore, we shall set $\mu_{m}=\mu_{0}$ throughout. A magnetic scalar potential $\psi$ is defined by ${\mathbf{H}}=\nabla\psi$, and satisfies

[TABLE]

For linear magnetic materials, the magnetization is a linear function of the magnetic field given by ${\mathbf{M}}=\chi{\mathbf{H}}$, where $\chi=(\mu/\mu_{0}-1)$ is the magnetic susceptibility. The magnetic induction ${\mathbf{B}}$ is therefore ${\mathbf{B}}=\mu_{0}({\mathbf{H}}+{\mathbf{M}})=\mu_{0}(1+\chi){\mathbf{H}}$. To describe the paramagnetic susceptibility quantitatively, the Langevin function $L(\alpha)$=coth $\alpha-\alpha^{-1}$ is used to describe the magnetization versus the magnetic field as

[TABLE]

where the saturation magnetization, $M_{s}$, and the magnetic moment of the particle, $m$, enter as parameters, $T$ denotes the temperature, and $k_{B}$ is the Boltzmann’s constant.

Let $[\![\cdot]\!]$ denote the difference, $\cdot_{\rm surrounding}-\cdot_{\rm drop}$, at the interface between the ferrofluid drop and the external liquid. Let ${\mathbf{n}}$ denote the unit normal outwards from the interface, and ${\mathbf{t}_{1}}$ and ${\mathbf{t}_{2}}$ denote the orthonormal tangent vectors. We require

The continuity of velocity, the normal component of magnetic induction, the tangential component of the magnetic field, the tangential component of surface stress,

[TABLE] 2. 2.

The jump in the normal component of stress balanced by capillary effects,

[TABLE]

where {\mathbf{\Sigma}}=p-\rho gz+\bigl{(}2\eta{\mathbf{D}}\bigr{)}+{\mathbf{\Pi}}_{M}, $k=-\nabla\cdot{\mathbf{n}}$ is the interface curvatures, and $\sigma$ is the interfacial tension.

3.2 One-fluid formulation and Volume-of-Fluid method

We take advantage of the axial symmetry present in the motion and deformation of a ferrofluid droplet placed in a viscous medium under an externally applied magnetic field and resort to the formulation in axisymmetric cylindrical coordinates $(r,z)$. The VOF method represents each liquid with a color function as

[TABLE]

which is advected by the flow. The position of the interface is reconstructed from $C(r,z,t)$. The one-fluid formulation of the governing equations (6) and (5) then reads

[TABLE]

where ${\mathbf{F}}_{s}$ denotes the continuum body force due to interfacial tension Brackbill92

[TABLE]

${\bf n}=\nabla C/|\nabla C|$, $\delta_{s}=|\nabla C|$ is the delta-function at the interface, $\rho=\rho(C)$, $\eta=\eta(C)$, and $\mu=\mu(C)$. The magnetic potential $\psi$ in axisymmetric cylindrical coordinates satisfies equation (12),

[TABLE]

The dimensionless variables, denoted by tildes, are defined as follows

[TABLE]

where $L_{c}$, $\tau$, $H_{o}$ and $g$ are characteristic scales of length, time, magnetic field and gravitational acceleration, respectively. The ferrofluid density and viscosity are ${\rho_{d}}$ and ${\eta_{d}}$, respectively. For the choice of viscous time scale,

[TABLE]

corresponding to the viscous length scale

[TABLE]

equation (18) becomes, dropping the tilde notation,

[TABLE]

where the Ohnesorge number $\mbox{Oh}=\eta_{d}/\sqrt{\rho_{d}\sigma L_{c}}$, the magnetic Bond number $\mbox{Bo}_{m}=\mu_{0}H^{2}_{o}L_{c}/\sigma$, and gravitational Bond number $\mbox{Bo}=\rho_{d}gL_{c}^{2}/\sigma$.

3.3 Spatial discretization

In axisymmetric coordinates $(r,z,\theta)$, the magnetic field is given by

[TABLE]

and the corresponding Maxwell stress is

[TABLE]

We have

[TABLE]

where $\bf{e_{r}}$ and $\bf{e_{z}}$ are unit vectors in $r$ and $z$ directions. The magnetic potential field is discretized using second-order central differences. In axisymmetric coordinates, the discretization of (20) at cell $(i,j)$ yields, on a regular Cartesian grid of size $\Delta$,

[TABLE]

where, for instance for the cell face $(i+1/2,j)$,

[TABLE]

The spatial discretization of the velocity field is based on the MAC grid in figure 7.

Therefore, the evaluation of the components of the magnetic stress tensor requires the evaluation of gradients at faces. In axisymmetric coordinates, the divergence of the magnetic stress tensor is discretized in the ${\bf e_{r}}$ and ${\bf e_{z}}$ directions as

[TABLE]

respectively, where, for example, components such as $(({\Pi}_{M})_{rr})_{i,j}$ and $(({\Pi}_{M})_{rz})_{i+1/2,j+1/2}$ are discretized as follows

[TABLE]

[TABLE]

respectively, where $\mu_{i+1/2,j+1/2}$ is computed using a simple averaging from cell center values. $({\Pi}_{M})_{rr}$ and $({\Pi}_{M})_{rz}$ at other grid locations are discretized similarly. Analogous relationships can be written for the other components of the magnetic stress tensor.

3.4 Drop deformation under a uniform magnetic field

When a ferrofluid drop is suspended in a non-magnetizable medium in an externally applied uniform magnetic field, it elongates in the direction of the applied magnetic field and assumes a stable equilibrium configuration achieved via the competition between the capillary and magnetic forces. Next, we show the numerical results of the VOF method for a ferrofluid drop suspended in a non-magnetizable viscous medium by Afkhami et al. ATRRWPR . The initial configuration for the computational study of a ferrofluid drop suspended in a non-magnetizable viscous medium is shown in figure 8. A uniform magnetic field ${\bf H}=(0,0,H_{o})$, where $H_{o}$ is the magnetic field intensity at infinity, is imposed at the top and bottom boundaries of the computational domain. In order to solve Laplace’s equation (20) in the presence of an interface, the following boundary conditions are employed: ${\frac{\partial}{\partial z}\psi}=H_{o}$ at $z=0,L_{z}$, and ${\frac{\partial}{\partial r}\psi}=0$ at the side boundary $r=L_{r}$. Note that a symmetry condition, ${\frac{\partial}{\partial z}\psi}=0$, at $z=0$, and ${\frac{\partial}{\partial r}\psi}=0$, at $r=0$ can be applied.

The drop radius is $R_{o}=1$ mm, interfacial tension is $\sigma=1$ mN m*-1*, and $\mu_{d}$ is chosen to be constant.

Figure 9(a) shows the results in ATRRWPR along with theoretical predictions. The details of the theoretical analysis is given in ATRRWPR . From the computational results alone, it is observed that at sufficiently large $\chi$, the drop equilibrium aspect ratio jumps to a higher value when the magnetic Bond number reaches a critical value, while for small values of $\chi$, the drop equilibrium aspect ratio increases continuously as a function of the magnetic Bond number, in agreement with the theories. Specifically, at $\chi=20$, the jump in drop shape occurs when the magnetic Bond number changes from $\mbox{Bo}_{m}=0.18$ to $0.19$. they also find that for a fixed Figure 9(b) shows the result for $b/a\approx 7$, $\mbox{Bo}_{m}=0.2$, and $\chi=20$ (see the line -.- in part a), where the shape exhibits conical ends. To provide further insight, the magnetic field lines inside the highly deformed drop as well as in the non-magnetizable surrounding medium, showing the effect of the magnetic field responsible for the appearance of conical ends.

3.4.1 Drop deformation under non-uniform magnetic fields

Here we present the motion of a hydrophobic ferrofluid droplet placed in a viscous medium and driven by an externally applied non-uniform magnetic field is investigated numerically in an axisymmetric geometry. This numerical investigation is motivated by recent developments in the synthesis and characterization of ferrofluids for possible use in the treatment of retinal detachment Mefford . Figure 10 shows a cartoon of the application of a small ferrofluid drop injected into the vitreous cavity of the eye and guided by a permanent magnet inserted outside the scleral wall of the eye. The drop travels toward the side of the eye, until it can seal a retinal hole.

To better understand the motion of the ferrofluid droplet moving in the eye, Afkhami et al. ARRRP present a numerical investigation of a more straightforward scenario where an initially spherical drop is placed at a height $L$ above the bottom of the cylindrical domain depicted in figure 8. In addition, the boundary condition in the figure is changed to reflect the presence of a magnet at the bottom, which instantly magnetizes the drop. The boundary condition on the magnetic field is reconstructed from the experimental measurements in Mefford , where the magnitude of $H(z)$ is measured as a function of distance from the magnet, $z$, in the absence of the drop. We then fit the data to a 5th degree polynomial, as shown in figure 11. The scalar potential therefore is a 6th degree polynomial $\phi(0,z)=P_{6}(z)$ along the axis of the cylindrical domain. In the absence of the drop, $\psi$ satisfies Laplace’s equation $\frac{1}{r}\frac{\partial}{\partial r}(r\frac{\partial\psi}{\partial r})+\frac{\partial^{2}\phi}{\partial z^{2}}=0$. If there is a solution, it is analytic and has $r^{2}$-symmetry. The ansatz $\phi(r,z)=P_{6}(z)+r^{2}P_{4}(z)+r^{4}P_{2}(z)+r^{6}P_{0}(z)$ yields

[TABLE]

This yields the boundary condition, and also approximates an initial condition when the drop is relatively small.

The magnetic potential $\psi$ is calculated from equation (20). The boundary conditions for $\psi$ on the domain boundaries $\partial\Omega$ are defined as

[TABLE]

where $\partial/\partial n_{b}={\bf n}_{b}\cdot\nabla$, and ${\bf n}_{b}$ denotes the normal to the boundary $\partial\Omega$. In order to impose the boundary condition in the numerical model, a transformation of variables is performed to $\zeta$: $\phi=\psi+\zeta$, where $\phi$ is the potential field without the magnetic medium. One can then rewrite equation (20) such that

[TABLE]

where $\nabla\cdot(\mu\nabla\phi)$ vanishes everywhere except on the surface between the drop and the surrounding fluid $\partial\Omega_{f}$ and

[TABLE]

Figure 12 shows direct numerical simulation results for one of the series of experiments conducted in Mefford : a droplet of 2 mm diameter placed 11 mm away from the bottom of the domain where the permanent magnet is placed. The figure shows the migration of the drop at times t = 120, 160, and 170 s along with the velocity fields. The travel time compares well with the experimentally measured one in Mefford . Additionally, the results show that at an early stage, the flow occurs approximately downwards and only in the region close to the droplet. As the droplet approaches the magnet, it elongates in the vertical direction and vortices induced in the viscous medium become stronger. When the droplet reaches the bottom of the domain, the flow inside the droplet is pumped outward from the center of the droplet, resulting in the flattening of the droplet and consequently a decrease in the droplet height. These observations are consistent with experiments in Mefford .

3.5 Magnetowetting of thin ferrofluid films

Figure 13 shows the phenomenon of magnetowetting for an experiment with magnetic droplets on a superhydrophobic surface, below which is a permanent magnet Timonen2013 . By gradually increasing the strength of the magnetic field and the vertical field gradient, the figure shows the transition of the the droplet from a spherical shape into a spiked cone and eventually splitting into two smaller droplets at a critical field strength.

In Seric2014 , a problem closely related to magnetowetting is investigated: the application of a uniform magnetic field to induce dewetting of a thin ferrofluid film. Here we present the study in Seric2014 , where thin film equations are derived using the long wave approximation of the coupled static Maxwell and Stokes equations and the contact angle is imposed via a disjoining/conjoining pressure model.

Figure 14 shows the schematic of the system of two thin fluid films in the region $0<y<\beta h_{c}$, with the ferrofluid film occupying the region $0<y<h_{c}$, and the nonmagnetic fluid occupying the rest of the domain. We denote the permeability and viscosity of the fluids by $\mu_{i}$ and $\eta_{i}$, where $i=(1,2)$ denotes fluids $1$ and $2$, respectively. The interface between two fluids is denoted by $y=h(x,t)$.

The unit normal and tangential vectors to the interface, n and t, respectively, are given by

[TABLE]

Given the small thickness of fluids, we ignore inertial effects and gravity. Hence, the equations governing the motion of the fluids are Stokes equations for continuity and momentum balance, $\nabla\cdot{\bf{u}}=0,\ \nabla\cdot\mathbf{T}=0,$ respectively, where $\mathbf{T}=-p\mathbf{I}+2\eta{\mathbf{D}}+\mathbf{\Pi}_{M}+\Pi_{W}\,\mathbf{I}$ is the total stress tensor and the disjoining pressure is specified by

[TABLE]

where $\bar{\kappa}={\sigma\tan^{2}\theta}/{\left(2M{h}_{*}\right)}$ and $M={\left(n-m\right)}/{\left[\left(m-1\right)\left(n-1\right)\right]}.$ In this form, $\Pi_{W}$ includes the disjoining/conjoining intermolecular forces due to van der Waals interactions. The prefactor $\bar{\kappa}$, that can be related to the Hamaker constant, measures the strength of van der Waals forces. Here, $h_{*}$ is the short length scale introduced by the van der Waals potential. We use $n=3$ and $m=2$. The contact angle, $\theta$, is the angle at which the fluid/fluid interface meets the substrate. The addition of the van der Waals forces is crucial in the present context, since it allows us to study dewetting of thin films under a magnetic field.

We nondimensionalizing the governing equations and boundary conditions using the following scales (dimensionless variables are denoted by tilde)

[TABLE]

[TABLE]

where $\delta\ll 1$. The initial thickness of the ferrofluid film is denoted $h_{c}$, and $u_{c}$ and $x_{c}$ are the characteristic velocity and horizontal length scale, respectively, given by

[TABLE]

Dropping the tilde notation for simplicity, the evolution equation for $h(x,t)$ then becomes

[TABLE]

where $\kappa={h_{c}\sigma\tan^{2}\theta}/{\left(2\mu_{0}\psi_{0}^{2}M{h}_{*}\right)}$ is a nondimensional parameter representing the ratio of the van der Waals to the magnetic force, $\mu_{r}$ is the ratio of the ferrofluid film permeability to the vacuum permeability, $\mu_{2}/\mu_{0}$, and $\alpha=0,1$ for Cartesian and cylindrical coordinates, respectively. The ratio $\beta\mu_{r}/(\mu_{r}-1)$ is inversely proportional to the magnetic force (note that $\beta$ is the nondimensional distance between the plates with constant potential, so the gradient of the potential is inversely proportional to $\beta$).

Here we show the thin film simulations of the steady state profiles obtained for a range of parameter values, in particular for nondimensional parameters $\beta$ and $\kappa$. We fix $\mu_{r}=44.6$ and $h_{*}=0.01$. The initial condition is set to a flat film perturbed around a constant thickness $h_{0}$, i.e. $h\left(x,0\right)=h_{0}+\epsilon\cos{\left(k_{m}x\right)}\,,$ with $\epsilon=0.1$, and $h_{0}=1$. Here, $k_{m}$ is the fastest growing mode computed from the linear stability analysis in Seric2014 . To put this in perspective, $\kappa=13.5$, $\beta=8$, and $\mu_{r}=44.6$ correspond to values of $\sigma$ and $\mu_{2}$ for an oil-based ferrofluid ($\sigma=0.034\mbox{ N/m}$, $\mu_{2}=\mu_{0}(1+\chi_{m})$, where $\chi_{m}=3.47\times 4\pi$), and $\psi_{0}=1.2\mbox{ A}$, and $\beta$ is chosen to produce a sufficiently strong magnetic field. Making use of the symmetry of the problem, the computational domain is chosen to be equal to one half of the wavelength of the perturbation, i.e. $L_{x}=\pi/k_{m}$. No flux boundary conditions are imposed at the left and the right end boundaries as $h^{\prime}=h^{\prime\prime\prime}=0$. Figure 15(a) shows the comparison of the steady state profiles for varying parameter $\kappa$ when keeping $\beta$ fixed. We note that as $\kappa$ decreases, i.e. the magnetic force increases compared to the van der Waals force, the satellite droplets start to appear. Similar behavior is observed in figure 15(b) for decreasing values of $\beta$; we observe the formation of satellite droplets for sufficiently small $\beta$, i.e. magnetic field dominates over the van der Waals interaction. It should also be noted that the satellite droplets are not present in the simulations where the effect of the magnetic field is ignored (i.e. when $\mu_{r}=1$). We also note that for $\beta=7.0$ shown in figure 15(b), a static drop cannot be obtained, unlike when $\beta>7.0$: here the height of the fluid approaches the top boundary and the assumptions of the model are not satisfied for the times later than the one at which this profile is shown.

The results of simulations suggest that the satellite droplets form when (i) the magnetic force is sufficiently strong, and (ii) the van der Waals force is sufficiently weak, relative to the magnetic force.

4 Conclusion

Biomedical technology is expanding the use of coated superparamagnetic nanoparticles and their suspensions in novel directions; for instance, the delivery of drugs to targeted cells, magnetic resonance imaging, early diagnosis of cancer, treatment such as hyperthermia, and cell separation. The forefront of research lies at the symbiotic development of nanomedicine, nanotechnology, engineering of strong and highly focused magnets, together with numerical algorithms. In this article, three simplified mathematical models have been presented with the aim of producing numerically generated solutions: stochastic differential equations linked to magnetic drug targeting, a volume-of-fluid computational scheme for the motion of a ferrofluid drop through a viscous medium under a magnetic field, and the evolution of a thin film of ferrofluid with numerical investigation of dewetting. An important future direction is the integration of open-source codes for the numerical simulations used with actual biomedical applications. A drawback of commercial software is that the location of numerical inaccuracy in any algorithm is difficult to pinpoint in a black box. The drive toward small scales and complex biomedical domains and materials necessitates a computational approach.