Understanding the dynamics of superparamagnetic particles under the influence of high field gradient arrays

We consider a spherical, superparamagnetic particle suspended in a flowing liquid under the influence of an approximately constant magnetic force generated by a linear Halbach array (figure 1). The magnetic force along the channel is simplified to a step-function bounded by a range, x_M from a starting position, x₀ and acting in the -z-direction

$$\begin{eqnarray}{{\mathbf{F}}_{\text{M}}}(x)=\left\{\begin{array}{*{35}{l}} -|{{F}_{\text{M}}}|\widehat{\boldsymbol{k}} & \quad \quad {{x}_{0}}\leqslant x\leqslant {{x}_{0}}+{{x}_{\text{M}}} \\ 0 & \quad \quad {{x}_{0}}>x>{{x}_{0}}+{{x}_{\text{M}}} \end{array}\right. \end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

In a steady flow, the forces acting on the particle are given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{m}_{\text{eff}}}\frac{\text{d}\mathbf{v}}{\text{d}t}={{\mathbf{F}}_{\text{M}}}+{{\mathbf{F}}_{\text{g}}}+{{\mathbf{F}}_{\text{buoy}}}+{{\mathbf{F}}_{\text{D}}}, \end{array}\end{eqnarray} \tag{ 2 }$$

where ${{m}_{\text{eff}}}$ is the mass of the particle including the added mass, $\mathbf{v}$ is the particle velocity and the four force terms on the right-hand side of the equation describe the magnetic force, gravity, buoyancy, and the fluid drag force respectively. We assume that ${{\mathbf{F}}_{\text{M}}}\gg {{\mathbf{F}}_{\text{g}}}+{{\mathbf{F}}_{\text{buoy}}}$ and ${{m}_{\text{eff}}}$ is sufficiently small that the particle reaches terminal velocity extremely rapidly (so ${{m}_{\text{eff}}}\text{d}v/\text{d}t\approx 0$). In a cylindrical channel carrying a laminar flow and positioned so that the central axis coincides with the x-axis, the flow velocity of the fluid is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathbf{v}}_{\text{F}}}=2{{v}_{\text{A}}}\left(1-{{\left(\frac{y}{R}\right)}^{2}}-{{\left(\frac{z}{R}\right)}^{2}}\right)\widehat{\boldsymbol{i}}, \end{array}\end{eqnarray} \tag{ 3 }$$

where v_A is the mean fluid flow velocity and R is the radius of the channel. Assuming a Newtonian fluid, the drag force on the spherical particle is given by Stoke's law:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathbf{F}}_{\text{D}}}=6\pi \mu r\left({{\mathbf{v}}_{\text{F}}}-\mathbf{v}\right), \end{array}\end{eqnarray} \tag{ 4 }$$

where μ is the dynamic viscosity of the fluid, r is the particle radius. Combining (2) and (4) gives

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathbf{v}={{\mathbf{v}}_{\text{F}}}+\frac{{{\mathbf{F}}_{\text{M}}}}{6\pi \mu r}={{v}_{\text{F}}}\widehat{\boldsymbol{i}}+{{v}_{z}}\widehat{\boldsymbol{k}}. \end{array}\end{eqnarray} \tag{ 5 }$$

The time-dependence of the position components, $x(t)$ and $z(t)$ is then determined by integrating the velocity components in (5) and setting a starting position vector described by $\mathbf{r}(t=0)=\left({{x}_{0}},y,{{z}_{0}}\right)$, resulting in

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & z(t)={{z}_{0}}+{{v}_{z}}t \\ {} & x(t)={{x}_{0}}+2{{v}_{\text{A}}}t\left(1-\frac{1}{{{R}^{2}}}\left(\frac{1}{3}{{\left({{v}_{z}}t\right)}^{2}}+{{z}_{0}}{{v}_{z}}t+{{y}^{2}}+z_{0}^{2}\right)\right). \end{array}\end{eqnarray} \tag{ 6 }$$

We assume a no-slip condition so that the final resting position of the particle coincides with where the particle makes contact with the bottom surface of the wall (as the external force is directed towards  −z), which occurs when $z\left({{t}_{\text{f}}}\right)={{z}_{\text{f}}}=-R\sqrt{1-{{\left(\frac{y}{R}\right)}^{2}}}$. Therefore, the final resting position is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{x}_{\text{f}}}-{{x}_{0}}=2\frac{{{v}_{\text{A}}}}{{{v}_{z}}}R\left({{\left(\frac{{{z}_{\text{f}}}}{R}\right)}^{2}}\left(\frac{{{z}_{\text{f}}}-{{z}_{0}}}{R}\right)-\frac{1}{3}\left({{\left(\frac{{{z}_{\text{f}}}}{R}\right)}^{3}}-{{\left(\frac{{{z}_{0}}}{R}\right)}^{3}}\right)\right), \end{array}\end{eqnarray} \tag{ 7 }$$

which has one real solution for z₀. The condition for which a particle can be considered to be magnetically retained is that the final resting position is within the bounds of the range of magnetic force, i.e. ${{x}_{\text{f}}}-{{x}_{0}}\leqslant {{x}_{\text{M}}}$.

We define the capture efficiency as the proportion of magnetized particles injected into the flow channel that are retained by the magnet. The capture efficiency was determined numerically by considering the final resting positions of an ensemble of particles initially distributed evenly about the cross section of the channel at the starting position, x₀. This particular case is analogous to a 'capture cross-section' (Cregg et al 2008), i.e. the proportion of streamlines that terminate under the influence of the magnet. Other initial distributions are considered in appendix C. A number of different combinations of parameters were tested and, in all cases, the capture efficiency was well described by the following power law

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{C}.\text{E}.\left({{x}_{\text{M}}}\right)=A{{\left(\frac{{{v}_{z}}}{{{v}_{\text{A}}}}\frac{{{x}_{\text{M}}}}{R}\right)}^{0.522}}\times 100 \%, \end{array}\end{eqnarray} \tag{ 8 }$$

where A  =  0.601 and assuming that $\text{C}.\text{E}.$ does not exceed $100 \%$. (This expression has a similar form to one derived for a two-dimensional case by Nandy et al (2008), though their equation has an additional component which is sensitive to the spatial variations of the field source they considered). The accumulation distribution, which we define as the proportion of particles that accumulate in a range of positions between $x\pm \frac{1}{2}\text{d}x$, is then

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\text{dC}.\text{E}.(x)}{\text{d}x}=0.31{{\left(\frac{{{v}_{z}}}{{{v}_{\text{A}}}R}\right)}^{0.522}}{{x}^{-0.478}}, \end{array}\end{eqnarray} \tag{ 9 }$$

and has units m⁻¹.

Typically, the magnetic force profile of a realistic magnetic system is not well described as a step-function. In order to predict the capture efficiency and accumulation distribution in MDT situations with a more physiologically relevant set of flow conditions and physically attainable magnetic forces, numerical simulations were performed to generate the trajectories of particles in a flow channel, under the influence of a position-dependent magnetic force. Simulations were performed using console applications written in the $C\sharp$ programming language (Microsoft Corporation, Redmond, WA, USA). Numerical simulations can be performed in 3D with a much higher number of particles than comparable models solved using finite element methods (which can often be computationally prohibitive or highly mesh-sensitive and inaccurate when solving for fluid flow fields in 3D geometries (Nacev et al 2011a, Barnsley et al 2016)); typical simulations involving 2000 particle trajectories and 14000 unique time iterations completed in 2 h 40 min on a computer with an Intel(R) Core(TM) i7 processor and 8 GB RAM. The trajectories were determined at each time interval, $\Delta t$ by updating the position, ${{\mathbf{r}}_{i+1}}={{\mathbf{r}}_{i}}+\mathbf{v} \Delta t$, with $\mathbf{v}$ calculated from an expression of a similar form to (5), but also taking into account gravity, buoyancy, and the spatial variation of ${{\mathbf{F}}_{\text{M}}}$:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathbf{v}={{\mathbf{v}}_{\text{F}}}\left({{\mathbf{r}}_{i}}\right)+\frac{{{\mathbf{F}}_{\text{M}}}\left({{\mathbf{r}}_{i}}\right)+{{\mathbf{F}}_{\text{g}}}+{{\mathbf{F}}_{\text{buoy}}}}{6\pi \mu r}. \end{array}\end{eqnarray} \tag{ 10 }$$

The magnetic force from an arbitrary arrangement of magnetized elements was calculated using a model described and experimentally verified previously (Barnsley et al 2015), whereby the field, $\mathbf{B}\left(\mathbf{r}\right)$ was determined by breaking the magnet into a 3-dimensional lattice of evenly-distributed point moments, and summing the contributed dipole field at $\mathbf{r}$ from each moment. In a curl-free magnetic field, the magnetic force on a particle with a magnetization, $\mathbf{M}$ and volume, V is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathbf{F}}_{\text{M}}}\left(\mathbf{r}\right)=\mathbf{M}V\centerdot \nabla \left(\mathbf{B}\left(\mathbf{r}\right)\right). \end{array}\end{eqnarray} \tag{ 11 }$$

The normalized magnetic force, used here for convenience because it has the same units as the field gradient (T m⁻¹), is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{F}_{\text{M}}}}{{{M}_{\text{s}}}V}=\frac{M}{{{M}_{\text{s}}}}|\nabla (B)|, \end{array}\end{eqnarray} \tag{ 12 }$$

where M_s is the saturation magnetization of the particle.

For simplicity, $\mathbf{M}$ and $\mathbf{B}$ were assumed to be parallel, and the particle was superparamagnetic so that the magnetization could be described using a Langevin function, $L\left(\,y\right)=\coth \left(\,y\right)-1/y$,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} M={{M}_{\text{s}}}L\left(\frac{{{M}_{\text{s}}}VB}{{{k}_{\text{B}}}T}\right), \end{array}\end{eqnarray} \tag{ 13 }$$

where ${{k}_{\text{B}}}T$ is the product of the Boltzmann constant and the temperature (Bean and Livingston 1959, Williams et al 1993, Majetich et al 1997).

For simulations using cylindrical channels, ${{\mathbf{v}}_{\text{F}}}$ was described analytically using (3). When a channel had a square cross-section, the flow profile was approximated by a ${{(z+iy)}^{4}}$ expression, following the formulation of Lekner (2007) and setting the wall topology parameter $\beta =1$,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathbf{v}}_{\text{F}}}=2{{v}_{\text{A}}}\left(\frac{4{{\left(\,{{y}^{2}}+{{z}^{2}}\right)}^{2}}}{{{w}^{4}}}\cos \left(4\phi \right)-\frac{4\beta \left(\,{{y}^{2}}+{{z}^{2}}\right)}{{{w}^{2}}}+1\right)\widehat{\boldsymbol{i}}, \end{array}\end{eqnarray} \tag{ 14 }$$

where w is the width of the channel and ϕ is the azimuthal angle. A 'sticky' condition was imposed for the wall, whereby particles that have made contact with the wall were assumed to be stable and no further updates were made to their position.

Numerical simulations were performed in which an ensemble of spherical particles modelling polystyrene microbeads loaded with a volumetric ratio (α) of superparamagnetic iron oxide clusters were evenly distributed at an injection point (Tehrani et al 2014, Bose and Banerjee 2015b, Lunnoo and Puangmali 2015) in a laminar flow contained within a straight channel aligned with the x-direction, and allowed to flow past a linear Halbach array consisting of nine rectangular, permanent magnet elements (figure 1). The magnet array was positioned a distance ${{z}_{\text{mag}}}$ away from the central axis of the channel and was designed so that, in the region of the channel above the magnet, the magnetic force was approximately, spatially uniform, provided the channel radius, R was sufficiently small. The main component of the magnetic force was aligned to the -z-direction, and gravity was directed perpendicular to both this direction and the direction of flow. A particle was considered captured if, when it came to rest at the channel wall, the magnitude of the magnetic force overwhelmed all others. The accumulation distribution was estimated as the relative proportion of particles lying within a range $x\pm \frac{1}{2}\text{d}x$. This was done by taking the derivative, with respect to x, of a cumulative distribution function of the proportion of particles that were stable before reaching x (i.e. ${{x}_{\text{f}}}<x$). Inter-particle interactions were neglected.

Table 1 gives the default values for simulation parameters. The magnet array tilt angle, θ was defined as the angle made between the bottom of the channel and the upper surface of the magnet. Of the 18 parameters in table 1, 14 are used to calculate the final term in (10). In most cases, the default values of these parameters were chosen to be relevant to the experiments described in section 2.3. The effects of varying v_A, F_M (which depends on ${{z}_{\text{mag}}}$), R and θ on capture efficiency and accumulation are investigated in section 3.1, while section 3.3 explores the performance of different magnet array designs.

Table 1. Parameters used in numerical particle tracing simulations.

Symbol	Description	Default value	Unit
μ	Fluid dynamic viscosity	$8.9\times {{10}^{-4}}$	Pa.s
vA	Fluid mean flow velocity	0.1	m s−1
${{\rho}_{\text{F}}}$	Fluid density	1000	kg m−3
R	Fluid channel radius	$0.25\times {{10}^{-3}}$	m
V	Particle volume	$\left(4\pi /3\right){{r}^{3}}$	m3
r	Particle radius	$0.5\times {{10}^{-6}}$	m
ρ	Particle density	$\alpha {{\rho}_{\text{F}{{\text{e}}_{3}}{{\text{O}}_{4}}}}+\left(1-\alpha \right){{\rho}_{\text{polysty}}}$	kg m−3
${{\rho}_{\text{F}{{\text{e}}_{3}}{{\text{O}}_{4}}}}$	Fe3O4 density	5240	kg m−3
${{\rho}_{\text{polysty}}}$	Polystyrene density	1050	kg m−3
α	Volumetric ratio of Fe3O4	0.1
${{M}_{s\text{F}{{\text{e}}_{3}}{{\text{O}}_{4}}}}$	Fe3O4 saturation magnetization	$4.7\times {{10}^{5}}$	A m−1
${{M}_{\text{NdFeB}}}$	NdFeB magnetization	$1.018\times {{10}^{6}}$	A m−1
${{d}_{\text{mag}}}$	Magnet element dimensions	$2.5\times 7\times 2.5$	mm3
${{z}_{\text{mag}}}$	Distance to magnet	3	mm
xM	Magnet array length	$9\times 2.5$	mm
θ	Magnet array tilt angle	0	°
${{F}_{\text{g}}}+{{F}_{\text{buoy}}}$	Gravity and buoyancy	$\left(\rho -{{\rho}_{\text{F}}}\right)Vg$	N
g	Gravity acceleration	10	m s−2

The ability of different magnet designs to target magnetic particles to a partially occluded vessel was also investigated by simulating trajectories within a symmetrically converging-diverging (or 'pinched') channel, shown in figure 2(b). The flow velocity inside a pinched channel was approximated following the derivation from Akbari et al (2010):

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathbf{v}}_{\text{F}}}=2{{v}_{\text{A}}}{{\left(\frac{{{R}_{0}}}{R(x)}\right)}^{2}}\left(1-{{\left(\frac{y}{R(x)}\right)}^{2}}-{{\left(\frac{z}{R(x)}\right)}^{2}}\right)\left[\widehat{\boldsymbol{i}}+\frac{\text{d}R(x)}{\text{d}x}\left(\frac{y}{R(x)}\widehat{\boldsymbol{j}}+\frac{z}{R(x)}\widehat{\boldsymbol{k}}\right)\right], \end{array}\end{eqnarray} \tag{ 15 }$$

where R₀ is the radius at the channel inlet and $R(x)$ is the radius at x. The channel has a constant slope on either side of the occlusion, $\text{d}R(x)/\text{d}x$.

Zoom In Zoom Out Reset image size

Numerical predictions of trajectories were verified experimentally by flowing polystyrene magnetic microbead particles (2.0–$2.9\times {{10}^{-6}}$ m, Spherotech, Inc., Lake Forest, IL, USA) in a glass, square capillary channel primed with water (figure 3). Microbeads were diluted to a concentration of $4\times {{10}^{5}}$ ml⁻¹ and injected into a steady-flow field using a syringe pump. A dilute concentration was chosen to diminish the formation of agglomerated beads in flow (superparamagnetic beads are observed to interact weakly when separated beyond approximately six times their particle diameter (Oduwole et al 2016)). Acoustic focusing of the microbeads into a single trajectory was achieved by generating an ultrasonic standing wave (USW) field inside the fluid cavity of the channel using a lead zirconate titanate (PZT) piezoelectric transducer driven at $2.55\to 2.65$ MHz, using a linear frequency sweeping regime. Additional information about the design of the ultrasonic manipulation setup is reported in appendix B. The trajectory of beads was magnetically deflected using a linear Halbach array constructed from nine rectangular N42-grade NdFeB permanent magnet elements (Magnet Expert Ltd., Tuxford, UK) arranged as per figure 1, with dimensions given in table 1. Images of the region of the channel adjacent to the magnet were acquired by an inverted microscope (Ti Eclipse, Nikon Instruments Inc., NY, USA) with a  ×10 objective (Plan Fluor Ph1 DLL, Nikon) and a digital camera (Digital Sight Qi1Mc-U2, Nikon), operating in bright-field mode.

Zoom In Zoom Out Reset image size

The trajectory of the microbeads actuated by the Halbach array were extracted by a purpose-built MATLAB (The MathWorks Inc., Natick, MA, USA) image processing routine. Imaging scans along the channel length were analyzed frame by frame. Each frame was thresholded to separate magnetic particles from background. Objects greater than 250 pixels in area were discarded to exclude agglomerated magnetic beads and objects with less than 12 pixels were discarded as noise or refuse. Separate threshold levels were applied to isolate the top and bottom channel walls. Linear regression was used to extrapolate the top wall coordinates from the first 8 frames of the scan (before the Halbach array) since thresholding with magnetic particles accumulating at the wall of the channel was less reliable. The origin was taken as the midpoint between the top and bottom channel walls at the beginning of the first frame of the scan and particle centroid coordinates were recorded relative to this point. Non-magnetic particles that satisfied the size and thresholding criteria were rarely detected and manifest as outliers in the plots in figure 9.

Magnetic measurements were used to characterize the magnetic behaviour of the microbeads in suspension. Hysteresis loops of the microbeads were measured using a MPMS superconducting quantum interference device (SQUID) magnetometer (Quantum Design, Inc., San Diego, CA, USA) at a temperature of 270 K (in order to freeze the fluid) between  ±1 T. The effective superparamagnetic particle size and volumetric loading of magnetite on microbeads were determined by fitting (13) to the hysteresis loop in figure A1, and these parameters were used to numerically predict the experimental trajectories reported in section 3.2.

Numerical simulations were performed to ascertain the trajectories of magnetic microbeads in a steady, laminar flow, under the influence of the magnet array shown in figure 1, to investigate how the capture efficiency and accumulation distributions varied with the flow velocity, the magnetic force, the channel radius and the alignment of the magnet. The field and force profiles emitted by the Halbach array along the direction parallel with the upper surface, displayed in figure 4, show that, for a short range of z values, the 'step function' approximation used for the derivation in section 2.1 was reasonable. The results from the numerical simulations were compared with predictions using the expressions for trajectories, C.E. and accumulation distribution derived analytically.

Zoom In Zoom Out Reset image size

Figure 5 compares the analytically derived trajectories described by (6) and the simulated trajectories for average flow velocities of 10 and 20 mm s⁻¹ and ${{z}_{\text{mag}}}$ set to 3 mm. The normalized magnetic force for the analytical trajectories was taken as the force directly above centre of the array at ${{z}_{\text{mag}}}$ (${{F}_{\text{M}}}=93.9$ T m⁻¹). While the analytical expression made reasonable estimates of the final resting position of particles in flow, a few discrepancies are apparent, owing to the assumptions made in the derivation. The analytical trajectories generally underestimated x_f on the near (negative) side of the magnet, and overestimated it on the positive side, which could be attributed to the variation of ${{F}_{\text{M}}}(x)$ produced above the sides of the Halbach arrray. More importantly, the trajectory of the particle with the highest starting elevation shown in figure 5(a) (z₀  =  2R/3) demonstrates how a particle could be both actuated and captured outside of the bounds of x_M when a more physically realistic profile of ${{F}_{\text{M}}}(x)$ was utilized. Once a particle was in a low flow regime, close to the channel wall, it slowed down significantly, and relatively little external force was required to pull it to the wall.

Zoom In Zoom Out Reset image size

In order to examine a wide range of possible trajectories, numerical simulations were run by initially distributing an ensemble of particles evenly around the cross-section of the flow channel upstream from the magnet, so that each particle had a unique starting ${{y}_{0}},{{z}_{0}}$. Figure 6(a) shows how the travel distance of simulated particles (relative to the centre of the magnet, which defines x  =  0) is related to their the initial position in the y-z plane at the injection point, with an average laminar flow velocity of 20 mm s⁻¹. The colours indicate the final resting positions of the particles, x_f, relative to the centre of the magnet. The particles that were initially furthest from the magnet (coloured white), and therefore had furthest to travel before being trapped at the channel wall, tended to pass over the magnet, reaching the outlet uncaptured. Alternatively, particles that had a very slow velocity near the bottom of the channel (coloured black) were removed from the flow by gravity before they reached the magnet. Analytical predictions for x_f as a function of y and z₀ using (7) for the same parameters are given in figure 6(b), and are qualitatively similar to the results of the numerical simulation. In both data sets, the proportion of particles trapped by the magnet in each 5 mm segment along the length of the channel (bounded by the counter lines) decreases with increasing x, indicating that more particles were trapped on the leading edge side (negative x_f) of the magnet, and fewer made it to the trailing edge (positive x_f).

Zoom In Zoom Out Reset image size

Capture efficiencies were determined for a set of numerical simulations following the method described in section 2.2 and were compared with the values predicted by (8). Figure 7 shows how the capture efficiency varied as one of the experimental parameters, v_A, F_M, R or θ was changed while the other parameters were fixed to their default values (table 1). The agreement attained between the simulation results and the analytical predictions indicates that the formulation derived in section 2.1 can provide a good first-order approximation for the capture efficiency in regions where the magnetic force is (relatively) spatially uniform.

Zoom In Zoom Out Reset image size

The discrepancy in figure 7(a) at low flow velocities was due to the influence of gravity in the numerical simulations causing a relatively large proportion of particles to be removed from flow by gravity before reaching the magnet. At flow velocities less than 10 mm s⁻¹, no particle made it to the outlet of the channel; particles were either forced out of flow by gravity before reaching the magnet, or captured by the magnet once they reached it. If simulations with these conditions did not include gravity, the capture efficiency would be 100%, as suggested by the analytical expression. Figure 7(c) shows that, as the channel size increased, so did the deviation between the two models. This was attributed to the fact that, as the channel size increased, the assumption of a spatially uniform magnetic force was no longer valid, and the magnetic force at the bottom of the channel, closer to the magnet, became significantly stronger than the force at the top of the channel. The effect of varying the magnet array tilt angle, θ is displayed in figure 7(d). For $\theta >{{22.5}^{\circ}}$, the capture efficiency didn't vary significantly, but the accumulation distribution (figure 8(d)) did.

Zoom In Zoom Out Reset image size

Accumulation distributions, calculated for a subset of numerical simulations using the method described in section 2.2, are shown in figure 8, and are compared with the profiles predicted by (9). For all tested parameter sets, (9) goes to infinity at the leading edge of the magnet, x₀, highlighting a major discrepancy between the two models. Almost all distributions exhibited a peak in density in a region close to the leading edge of the magnet and, downstream from this point the agreement between the numerical simulations and the analytical derivation was reasonably good in most cases where the magnetic force was spatially uniform. Exceptions to this can be observed in figure 8(b), where the channel was far from the magnet, and the normalized force was small, and in figure 8(c) for large channel radii, where the forces at the top and bottom of the channel were significantly different. Notably, in many simulations, a significant proportion of particles were captured in a region prior to the bounds of the spatial extent of the magnet. Also of note, the peak in density usually occurred in a consistent position for almost all distributions ($x\sim -8.5$ mm) and, provided the simulations were performed in a space where the magnetic force could be approximated by (1), the shapes of the profiles overlapped very closely when normalized to the peak density. This suggests that, while the proportion of captured particles is sensitive to physiological flow conditions inside the body (encapsulated in the parameters v_A and R), the relative spatial distribution of particles that are successfully captured is determined by the profile of the magnetic force (${{F}_{\text{M}}}\left(\mathbf{r}\right)$). This means that particle distributions could be tuned and optimized for a given application by engineering the spatial profile of the magnetic force, so that the greatest density of captured particles would coincide with the desired region of interest. One attempt to do this is shown in figure 8(d), by varying θ. By tilting the magnet away from the channel, so that the force was localized in the region above the trailing edge of the magnet, the peak position of accumulation distribution was shifted closer to the position $x={{x}_{\text{M}}}$. However, fewer particles were captured overall.

The numerically predicted trajectories of particles under the influence of laminar flow and magnetic force were verified experimentally by flowing magnetic microbeads in a glass square capillary channel with an internal width of 0.3 mm and focusing them using acoustic radiation forces generated by a USW field, upstream from a linear Halbach array. The USW field was created by a PZT transducer coupled to the glass channel approximately 8 mm away from the Halbach array, to separate the regions of the channel exposed to acoustic and magnetic forces (appendix B). Thus, the contribution of acoustic radiation forces was considered negligible when modelling the trajectory of particles exiting the active area of the transducer. The fluid velocity inside the channel was described using (14), with the average velocity, ${{v}_{\text{A}}}$ varied between 2 and 15 mm s⁻¹ (a range of flow velocities seen in intratumoral blood flow (Shimamoto et al 1987)). Starting positions used to calculate numerical trajectories were determined by averaging the elevation and distance to the wall of particles upstream from the influence of the magnet. The magnet position relative to the channel was measured from scanned images recorded using a fluorescent microscope, with the magnet-to-channel distance, ${{z}_{\text{mag}}}$ set at 3.9 or 4.6 mm (corresponding to field gradients of 45 and 23 T m⁻¹ above the centre of the array).

Figure 9 shows the measured trajectories compared with the numerical predictions, based on a microbead with the magnetic properties measured in appendix A. For each condition, reasonable agreement was found between the measured and predicted trajectories. The shaded regions in figure 9 represent the region of the channel adjacent to the magnet, but for low flow rates and/or higher magnetic forces, it was both observed and predicted that the trajectory would begin to deflect towards the magnet before the particles reached this area. At faster flow velocities, particles were less well focused due to them spending less time in the acoustic field of the transducer, resulting in a greater spread of particle positions around the predicted trajectories. For example, at 2 mm s⁻¹ and 45 T m⁻¹ (figure 9(a)), focusing confined the particles within a trajectory about 0.032 mm wide (±$2\sigma$) before being actuated by the magnet, resulting in a 1.2 mm spread of positions for particles that were forced to the wall by the magnet. By comparison, a 15 mm s⁻¹ average velocity resulted in a 0.056 mm wide trajectory, and a 8.1 mm spread of final particle positions for the same magnetic force (figure 9(d)). Variations in particle size and magnetic loading may have also contributed to distributions in particle positions about the predicted trajectory, as they would impact on both acoustic and magnetic forces.

Zoom In Zoom Out Reset image size

Simulations were performed to investigate how accumulation distributions would differ with different permanent magnet designs, in order to determine which shapes and arrangements result in the greatest concentration of particles in a localized region of interest. Figure 10 shows the six magnet designs explored here. All designs were constrained to have the same magnet volume as the linear Halbach array and, in the case of the button and rod magnets, the same spatial extent along x. For the cone magnet, the diameter and height were set to be equal. The cube magnet (figure 10(d)) had a magnetization vector that pointed along a face diagonal, and was orientated on its edge so that $\mathbf{M}$ pointed along z. The arrays in figure 10(e) and (f) were designed using a force optimization routine for arbitrary shapes (Barnsley et al 2016) and a position of interest set ${{z}_{\text{mag}}}=3$ mm away from the upper surface of the magnet. One design consisted of a single, uniform magnetization, while the other was a Halbach array assembled using six segments with different, orthogonally aligned, magnetization directions.

Zoom In Zoom Out Reset image size

The field and force profiles generated at a distance ${{z}_{\text{mag}}}=3$ mm from the upper surface of each magnet are displayed if figure 11, which shows that the uniform button and cone magnets are not well suited to apply a high field gradient in a localized region (the field and force from the button magnet are strongest at the edges). The two optimized magnets generate the strongest fields in the centre of the channel, particularly the Halbach array, which has a full width at half maximum (FWHM) in force of 8.54 mm.

Zoom In Zoom Out Reset image size

Numerical simulations were performed using the different magnet designs with the default values for v_A, ${{z}_{\text{mag}}}$, R and θ (table 1). The channel was positioned so that the maximum flow axis was ${{z}_{\text{mag}}}=3$ mm from the upper face of the magnet array. Accumulation distributions in figure 12 indicate that the optimized arrays are better suited for localizing particles in greater density around the region directly above the centre of the magnet. Of the uniform magnets, the button magnet was notable because most accumulation took place above the edges of the cylinder, where the forces were strongest, and almost no particles were collected towards the middle. The rod behaved similarly to the step-force (linear Halbach array) above, capturing more particles at the leading edge of the magnet, but still capable of distributing particles in the region bound by x_M, suggesting that these designs could be practical for applications in which particles need to be delivered to a large, diffuse area. These include diagnostic or biosensing applications involving a combination of active and magnetic targeting strategies for precisely locating and imaging tumours (Gunasekera et al 2009, Schleich et al 2014, Xu et al 2016). The cone and cube magnets were able to concentrate particles to a confined region close to their upper edges, but the peak linear density was not competitive with that of the two optimized arrays, owing to the fact that the generated force profiles persisted only over a very short-range compared with the optimized arrays (Barnsley et al 2015).

Zoom In Zoom Out Reset image size

The ability of each magnet design to deliver particles to a partial occlusion is assessed in figure 12(b). In all cases, the accumulation distribution at the site of the occlusion is less than the 'no occlusion' case in figure 12(a), as the flow velocity in the narrow region of the channel is significantly greater, and the distance from the bottom of the channel to the magnet has increased, reducing the peak magnetic force inside the channel. Table 2 details the capture efficiencies from the simulations shown in figure 12. While the optimized magnets are still the best of the tested designs at delivering particles to the position of interest, they are also the designs most relatively compromised by the presence of the occlusion.

Notably, even when the accumulation distribution was tightly localized, it was still skewed slightly towards the leading edge of the applied field/force. Figure 13 shows how the accumulation distribution varies as a function of two physiologically-relevant parameters, v_A and R in a channel with no occlusion. Again, when v_A was varied (and all magnet related parameters were fixed), the shape of the accumulation distribution changed only very slightly. The accumulation distributions displayed in figure 13(a) overlap almost entirely when normalized to the peak density, with the position of the peak shifting by just  ∼1.4 mm over three orders of magnitude in v_A (the mean FWHM of all distributions was $8.9\pm 0.1$ mm). The shape of the distribution did change with R, particularly at large values toward 2.5 mm (figure 13(b)), owing to the fact that the bottom of the channel was only 0.5 mm from the magnet, and the force profiles at the top and bottom of the channel were very different.

Zoom In Zoom Out Reset image size