Two-Phase Biofluid Flow Model for Magnetic Drug Targeting

Abstract

Magnetic drug targeting (MDT) is a noninvasive method for the medical treatment of various diseases of the cardiovascular system. Biocompatible magnetic nanoparticles loaded with medicinal drugs are carried to a tissue target in the human body (in vivo) under the applied magnetic field. The present study examines the MDT technique in various microchannels geometries by adopting the principles of biofluid dynamics (BFD). The blood flow is considered as laminar, pulsatile and the blood as an incompressible and non-Newtonian fluid. A two-phase model is adopted to resolve the blood flow and the motion of magnetic nanoparticles (MNPs). The numerical results are obtained by utilizing a meshless point collocation method (MPCM) alongside with the moving least squares (MLS) approximation. The numerical results are verified by comparing with published numerical results. We investigate the effect of crucial parameters of MDT, including (1) the volume fraction of nanoparticles, (2) the location of the magnetic field, (3) the strength of the magnetic field and its gradient, (4) the way that MNPs approach the targeted area, and (5) the bifurcation angle of the vessel.

1. Introduction

Nanoparticle-mediated drug delivery systems (DDS) appear to be a novel approach for drug delivery to diseased parts of the vasculature system [1]. The most commonly used drug nanocarriers for the MNPs are protein molecules, peptides, and deoxyribonucleic acid (DNA). Targeting agents and permeation enhancers are mounted on the surface of MNPs to improve the therapeutic outcome [2]. For atherosclerotic cardiovascular disease, the therapeutic goal is to reduce the likelihood of myocardial infraction and stabilize the atherosclerotic plaque [1].

In biomedical applications, there is a variety of different nanoparticles used [3,4,5,6]. The most commonly used are two iron oxides, namely magnetite $\left({\mathrm{Fe}}_3{\mathrm{O}}_4\right)$ and maghemite $\left(\gamma -{\mathrm{Fe}}_2{\mathrm{O}}_3\right)$ [7]. Their size ranges from 10 to 1000 nm [8], and they exhibit superparamagnetic properties [9]. An externally applied magnetic field, enhances their magnetic flux (their magnetic moments align with the magnetic field). They are not toxic, they are biocompatible at normal concentration (iron inherently occurs in many human organs including the heart) [6], they exhibit high saturation magnetization, and they have high magnetic susceptibility [10]. One of the most interesting properties of the MNPs is that, if an external magnetic field is applied, they can be easily removed from the blood stream [9,11]. The main disadvantages of MNPs are oxidation and aggregation [11,12]. Additionally, although they exhibit superparamagnetic properties, they tend to agglomerate due to their high surface energy. To reduce agglomeration, MNPs are coated with specific materials [11,12,13].

Biomagnetic fluid dynamics (BFD) studies the flow of biological fluids under the influence of an external magnetic field. There is an increasing interest in BFD studies, mainly due to the numerous bioengineering applications [14,15,16]. BFD modelling is closely related to ferro-hydro-dynamics (FHD) [17,18,19]. In BFD flow studies, the flow is strongly affected by the fluid magnetization under the presence of an applied magnetic field. In BFD models, the biofluid is a poor conductor and the Lorentz force is omitted. The first BFD model was introduced by Haik et al. [14]. Loukopoulos and Tzirtzilakis [15] studied the BFD in a channel under the influence of an external magnetic field, applying a single-phase fluid model. Tzirtzilakis and Loukopoulos in [20] extended the BFD flow model to Magnetohydrodynamics (MHD), such that the magnetization and electrical conductivity of the blood were included in the governing equations. Tzirtzilakis [21] reported on the biomagnetic flow in a channel with stenosis. The flow symmetry breaks downstream far from the stenosis, while the vortex located close to the magnetic field source increases in length. The flow reattachment is shifted far from the stenosis. Khashan et al. [16,22] combined the hydrodynamics and magnetic forces using a two-phase fluid model. They reported that the flow separates into a region located close to the magnetic source with high particle concentration, and a region far from the magnetic source with low particle concentration. In the work of Habibi et al. [23], the flow was considered pulsatile, the blood as non-Newtonian, and the nanoparticles were injected into the bloodstream. They report that the size of the nanoparticles does not affect the absorption (capture) of the particles. Additionally, during the systole MNPs were washed-out from the blood stream, decreasing the concentration of the nanoparticles. Comparing with their previous study [24], they concluded that there are major differences between steady and pulsatile flow in the flow regime, the particle absorption and the shear stress. Bose et al. [25] studied the magnetic particle capture in stenosed aortic bifurcation considering particle-fluid coupling. They reported that the magnetic field affects the flow regime. Larimi et al. [26] studied the flow of magnetic nanoparticles in a bifurcation vessel. They reported that the MNPs volume fraction in the region of interest depends on the intensity of the magnetic field.

The different phases, blood (continuum), and nanoparticles (dispersed phase) are modeled using a Eulerian or a hybrid or Eulerian–Lagrangian method. In the Eulerian formulation, both phases are treated as a continuum, while in the Euler–Lagrange formulation, the blood is treated as a continuum and the nanoparticles as discrete particles, such that the equation of motion is solved individually for each particle. In most of the BFD and FHD studies [14,15,20] assume that the base fluid-nanoparticles mixture interacts with the magnetic force as single-phase fluid. The base-fluid and the particles compose a colloidal suspension with infidelity strong coupling. Under the presence of an external magnetic field, the particles are not subjected to drifting and they do not re-distribute in the base fluid. On the other hand, the two-phase models in [22,27,28] include magnetophoretic motion relative to the base fluid, which results in a re-distribution of the particles when the Navier–Stokes equations are coupled with the advection-diffusion concentration equations.

In the present study, we investigate the two-phase (Euler-Euler) flow [27,28] of a pulsatile, non-Newtonian biomagnetic fluid (blood and dispersed MNPs) under the influence of a non-uniform magnetic field. In our simulations we use iron oxide $\left({\mathrm{Fe}}_3{\mathrm{O}}_4\right)$ as MNPs. We take into consideration three different geometries which simulate cardiovascular vessels, a planar microchannel, a microchannel with stenosis, and a microchannel with bifurcation. We examine significant parameters that affect the MDT such as: (1) the volume fraction of nanoparticles (2) the location of the magnetic field (3) the strength of magnetic field and its gradient, (4) the way that MNPs approach the targeted area, and (5) various bifurcation angles of the vessel. In any case, the time average volume concentration (TAVC) and the time average wall shear stress (TAWSS) indices are presented for the period of one cardiac cycle. Finally, contours of MNPs’ dimensionless concentration and stream function are presented at peak systolic velocity time moment.

2. Mathematical Model

2.1. Governing Equations

We study the transient, incompressible, laminar, non-Newtonian blood flow in two-dimensional flow domain (microchannels). We use iron oxide (Magnetite, ${\mathrm{Fe}}_3{\mathrm{O}}_4$) as the dispersed MNPs, which are spherical in shape and uniform in size. We assume that both blood and ${\mathrm{Fe}}_3{\mathrm{O}}_4$-MNPs are electrically non-conducting materials in thermal equilibrium. The ${\mathrm{Fe}}_3{\mathrm{O}}_4$-MNPs display linear magnetization without any initial or remnant magnetization. Buoyancy is negligible.

The flow equations in their primitive variable (u-p) formulation are as follows:

$$\nabla \mathit{u}=0$$

(1)

$${\rho }_{eff}\left(\frac{\partial u}{\partial t}+\left(\mathit{u}·\nabla \right)\mathit{u}\right)=-\nabla p+\mu {\nabla }^2\mathit{u}+\mathit{F}$$

(2)

Equation (2) in stream function-vorticity ($\psi -\omega$) formulation is written as [29]:

$$\frac{\partial \omega }{\partial t}+\frac{\partial \psi }{\partial x}\frac{\partial \omega }{\partial y}-\frac{\partial \psi }{\partial y}\frac{\partial \omega }{\partial x}=\frac{\mu }{{\rho }_{eff}}\left(\frac{{\partial }^2\omega }{\partial x^2}+\frac{{\partial }^2\omega }{\partial y^2}\right)+\frac{1}{{\rho }_{eff}}{\left(\nabla \times \mathit{F}\right)}_z$$

(3)

with $\omega =\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$ and $\psi$ computed using $u=\frac{\partial \psi }{\partial y}, v=\frac{\partial \psi }{\partial x}$, where ${\mathit{v}}_f=\left(u, v\right)$ is the fluid velocity. The nanofluid density ${\rho }_{eff}$ is a function of the volume concentration $\phi$ of the MNPs [30]:

$${\rho }_{eff}=\phi {\rho }_p+\left(1-\phi \right){\rho }_f$$

(4)

with ${\rho }_p, {\rho }_f$ being the density of the MNPs and blood, respectively. The force term $\mathit{F}=\left(F_x, F_y\right)$, which represents the effect of the magnetic field on MNPs concentration, is written as:

$$F_x=n_p{F}_{mag}^x$$

(5)

$$F_y=n_p{F}_{mag}^y$$

(6)

where $n_p$ is the volumetric MNPs density defined as $n_p=C_v/V_p$ ($C_v$ is the volume concentration and $V_p$ is the volume of MNP). The force volume components $F_x$ and $F_y$, are computed by multiplying the number of particles per unit volume $n_p$ with the magnetic force $F_{mag}^x$ and $F_{mag}^y$ acting on a single particle. Under the presence of an external magnetic field, magnetic particles will be magnetized. According to [16], [31] force ${\mathit{F}}_{mag}$ can be expressed as:

$${\mathit{F}}_{mag}={\mu }_0\left(\mathit{m}·\nabla \right)\mathit{H}$$

(7)

The force ${\mathit{F}}_{mag}$, considering the volume $(V_p)$ of MNPs, is expressed in term of the volume magnetization M as:

$${\mathit{F}}_{mag}={\mu }_0{V}_p\left(\mathit{M}·\nabla \right)\mathit{H}$$

(8)

The magnetization property $\mathit{M}$ for isothermal cases is described by a linear equation $\mathit{M}=\chi \mathit{H}$, with $\mathit{H}$ being the external magnetic field [16,32]. Equation (8), using $\nabla ·\mathit{H}=0$ and $H^2=\left(\mathit{H}·\mathit{H}\right)$, is written as:

$${\mathit{F}}_{mag}=\frac{1}{2}{V}_p{\mu }_o\chi \nabla H^2$$

(9)

with ${\mu }_0$ being the magnetic permeability of vacuum $\left({\mu }_0=4\pi \times {10}^{-7}\right)$, $\chi$ the magnetic susceptibility of the MNPs $\left(\chi =0.04\right)$, and $H$ the magnitude of the external magnetic field. The equation for the non-dimensional volume concentration is written as:

$$\frac{\partial C}{\partial t}+\nabla ·\left(C{\mathit{v}}_p\right)=\nabla ·\left(D\nabla C\right)$$

(10)

where $D$ is the diffusion coefficient of the MNPs and it is determined using the Einstein relation [33]:

$$D=\frac{k_{B}T}{6\pi \mu r_p}$$

(11)

$k_{\beta }$ is the Boltzmann constant, $T$ is the temperature of the blood $\left(=310 \mathrm{K}\right)$, $\mu$ is the dynamic viscosity, and $r_p$ is the MNPs’ radius. The MNPs velocity $({\mathit{v}}_p)$ is calculated as [16,23]:

$${\mathit{v}}_p={\mathit{v}}_f+{\mathit{v}}_{mag}$$

(12)

considering the superposition of the carrier fluid velocity $v_f$ and the relative magnetic particle’s velocity ${\mathit{v}}_{mag}$. Assuming that, after a short time period, there is a static balance of magnetic and viscous forces, the magnetic velocity can be calculated form the equation ${\mathit{v}}_{mag}=\frac{{\mathit{F}}_{mag}}{6\pi \mu r_p}$, by taking into account the Stokes drag law. The magnetic field intensity is given [32] by:

$$H_x=\frac{\gamma }{2\pi }\frac{x-x_{mag}}{{\left(x-x_{mag}\right)}^2+{\left(y-y_{mag}\right)}^2}$$

(13)

$$H_y=\frac{\gamma }{2\pi }\frac{y-y_{mag}}{{\left(x-x_{mag}\right)}^2+{\left(y-y_{mag}\right)}^2}$$

(14)

$$H=\sqrt{\left(H_x^2+H_y^2\right)}=\frac{\gamma }{2\pi }\frac{1}{\sqrt{{\left(x-x_{mag}\right)}^2+{\left(y-y_{mag}\right)}^2}}$$

(15)

with $\gamma$ being the magnetic field strength at the point where the magnet source is located $\left(x_{mag},y_{mag}\right)$.

2.2. Blood Viscosity—Power Law Model

The blood is considered as a non-Newtonian fluid. In our simulations we use the Power Law model [23,34], where the relation between the shear stress τ and the shear rate $\dot{\gamma }$ is given by $\tau =\mu \dot{\left(\gamma \right)}$. The shear rate is as follow $\dot{\gamma }=\left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)$. The tensor of shear is $\mathit{\tau }=\mu \left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)$. The viscosity $\left(\mu \right)$ for the Power Law model is given by the relation $\mu =m{\left(\dot{\gamma }\right)}^{n-1}$. Considering all the above, we conclude:

$$\mu =m{\left\{2\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]+{\left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)}^2\right\}}^{\frac{n-1}{2}}$$

(16)

which involves two material parameters, the consistency coefficient $m$ and the power-law index $n$.

2.3. Dipole-Dipole Interaction of MNPs

As mentioned above, the agglomeration of iron oxides is one of the main disadvantages. Gregg et al. [35,36] developed a model that includes both magnetic dipole-dipole and hydrodynamic interaction between magnetic drug carriers. The magnetic effect on n MNP of the other n-1 MNPs can de expressed as:

$${\mathit{F}}_{m{m}_n}=\left(m_n\nabla \right){\mathit{B}}_{tota{l}_n}$$

(17)

$${\mathit{B}}_{tota{l}_n}=\mathit{B}+{\displaystyle \sum }_{\begin{array}{c}i=1\\ i\ne n\end{array}}^{N}d{\mathit{B}}_i$$

(18)

where $m_n$ is the total magnetic moment, $\mathit{B}$ is the magnetic flux density as a result of the applied external magnetic field and $d{\mathit{B}}_{\iota }$ is the alteration of the consequential magnetic flux density. Nevertheless, in this study the effect of dipole–dipole interaction is negligible, hence it is not considered.

3. Algorithm Verification

We numerically solve the flow equations using the meshless point collocation (MPC) method. We represent the flow/spatial domain with a set of nodes, uniformly or randomly distributed over the interior domain and on the boundaries. We approximate the unknown field functions and their spatial derivatives using the moving least squares (MLS) method [37]. MLS method has been successfully applied to 2D flow problems [38,39,40,41,42]. We use a 4th order Runge–Kutta (RK4) time integration scheme [29] to deal with the transient term.

We demonstrate the convergence of the numerical method through a grid independence analysis study for all the flow cases considered. In this section, we present the convergence analysis study for the flow in a planar microchannel, the test case in the Section 4.1. The flow is laminar and non -Newtonian (following the Power-Law model). The density of the blood and MNPs are set to ${\rho }_f=1010 {\mathrm{kg}\mathrm{m}}^{-3}$ and ${\rho }_p=5180 {\mathrm{kg}\mathrm{m}}^{-3}$, respectively. At the inlet of the microchannel, we set a pulsatile velocity profile [43]. Additionally, we apply a non-uniform magnetic field (generated by a wire carrying electric current) with magnetic field strength $\gamma =3\times {10}^3\mathrm{A}$. The magnetic source is located at $\left(x_{mag},y_{mag}\right)=\left(0.01\mathrm{m},0.003\mathrm{m}\right)$. The MNPs volume fraction is $\phi =0.03$ and their radius is $r_p=5 \mathrm{nm}$.

We use successively denser grids and we examine the convergence of both stream function and vorticity field. Additionally, we examine the convergence of derived quantities, such as time average wall shear stress (TAWSS), which involves computation of spatial derivatives of the unknown field functions. We use five grid configurations, starting with a coarse grid with spatial resolution $h={10}^{-4}\mathrm{m}$ ($200\times 10$ uniform grid with 2000 nodes in total), and a denser grid with resolution $h=1.25\times {10}^{-5}\mathrm{m}$ ($1600\times 80$ uniform grid and 128,000 nodes). Figure 1 shows the TAWSS for all grid configurations used. The difference for the two finest grids is less than 1%, indicating that the resolution higher than $h=2\times {10}^{-5}\mathrm{m}$ ensures a grid independent numerical solution.

We demonstrate the accuracy of the meshless method, by comparing the numerical results obtained with those reported in Loukopoulos and Tzirtzilakis [15], using the finite difference (FD) method. Figure 2 shows the velocity profile at $x=3.5$ computed using the MPC and the FD method [15].

We also highlight the accuracy of our meshless scheme by comparing our time accurate solutions with those computed using the commercial CFD software ANSYS FLUENT^® (ANSYS Inc., Canonsburg, PA, USA) [44]. The flow domain is a rectangular channel of length $L=0.02 \mathrm{m}$ and width $H=0.001 \mathrm{m}$. The fluid flow is incompressible, laminar, non-Newtonian (power-law model with $n=0.8,\mathrm{m}=0.012$ (Equation (16)). The density is set to $\rho 1010 {\mathrm{kg}\mathrm{m}}^{-3}$ and the viscosity to $0.00309 \mathrm{Pa}\mathrm{s}\le \mu \le 0.02 \mathrm{Pa}\mathrm{s}$. We apply a sinusoidal velocity $u=0.0281sin\left(2.82485t\right)+0.0281 {\mathrm{ms}}^{-1}$ at the inlet. We compare the centerline velocity of the two methods for different time instances and the difference (for all three velocity vector components) is less than $1\%$. Figure 3 shows a comparison of the two methods at time $t=0.22243 \mathrm{s}$.

Finally, we collate our results for the two-phase flow model (Figure 4) with the numerical results computed using the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) (Figure 9a in [23]). The numerical results are in good agreement.

4. Results and Discussion

We numerically study MDT focusing on the flow regime and the concentration patterns of MNPs. For cardiovascular atherosclerosis, the therapeutic goal is to reduce the myocardial infraction and stabilize the atherosclerotic plaque. Biocompatible MNPs coated with drugs are navigated with a magnetic field to the atherosclerotic plaque.

In all flow cases considered, the flow is laminar and the blood is modeled as a non-Newtonian fluid using the power law viscosity model (Equation (17)), with material parameters values $n=0.8$ and $m=0.012$ (Equation (16)). Based on the experimental data in [23], viscosity obtains minimum and maximum values of $0.02 \mathrm{Pas}$ and $0.00309 \mathrm{Pas}$, respectively. We use the pulsatile velocity profile reported in [43]. The density of the blood and the MNPs is set to ${\rho }_f=1010 {\mathrm{kg}\mathrm{m}}^{-3}$ and ${\rho }_p=5180 {\mathrm{kg}\mathrm{m}}^{-3}$ [7], respectively. The radius of the MNPs, $r_p$, is equal to $5 \mathrm{nm}$.

We present results for different values of the solid volume fraction ($0.03\le \phi \le 0.09$), the magnetic field location ($0.001\mathrm{m}\le b\le 0.003\mathrm{m}$)—b is the vertical distance from the upper wall-, and the magnetic field strength (${10}^3 \mathrm{A}\le \gamma \le 7\times {10}^{3 }\mathrm{A}$). We examine three flow geometries, namely a planar microchannel, a microchannel with stenosis and a microchannel with bifurcation.

We calculate the time average volume concentration ($\mathrm{TAVC}={\int }_0^{T}Cdt,$ where T is the time period of one cardiac cycle), and the time average wall shear stress ($\mathrm{TAWSS}=\frac{1}{T}{\int }_0^T\mathrm{WSS} dt$, where WSS is the instantaneous wall shear stress, that is, the force per unit area exerted by the fluid on the wall in a direction on the local tangent plane every time). The TAVC shows the distribution for non-dimensional concentration of the MNPs. We present results along the upper wall of microchannel, where the magnetic source is located. High values for TAVC demonstrate increased MNPs capture close to the magnetic source. A drawback of MDT is the increased likelihood for embolization, since nanoparticles accumulate and block the blood flow [45]. High values for TAWSS indicate a high velocity gradient and maximum shear stress [26]. In arterial blood flow the WSS, which is the frictional force acting on the endothelial cell surface due to blood flow, is considered to be a major factor for the early formation of an atherosclerotic plaque [45,46]. We present stream function contours to illustrate the effect of MNPs to the flow field, and to detect recirculation zones (vortices). Shear rate and flow recirculation zones are associated with distinct pathogenic and biological pathways relevant to atherogenesis and thrombosis [47].

4.1. MDT in a Planar Microchannel

In this section, we study the biomagnetic flow in a planar microchannel. The flow domain and the boundary conditions are shown in Figure 5. The MNPs (drug carriers), enter the channel through the inlet. A detailed description of the concentration the boundary conditions is presented in [16,22]. We represent the spatial domain using 50,000 nodes, which ensure a grid independent numerical solution. We set the time step to $dt={10}^{-5}$ s and the MNPs volume fraction at the inlet to φ = 0.03.

We examine the effects of varying volume fraction $\left(\phi \right)$, distance $\left(b\right)$ of the magnetic field source from the upper wall of the microchannel, and magnetic field strength $\left(\gamma \right)$, affect the flow field, the distribution of dimensionless concentration along the upper wall, and the wall shear stress on the upper wall.

4.1.1. Effect of the MNPs Volume Fraction

Figure 6 shows the stream function contours for different values of magnetic MNPs volume fraction $\phi \left(0.03, 0.06 \mathrm{and} 0.09\right)$ at the peak systolic velocity ($t=0.8 \mathrm{s}$). The external magnetic field, with strength $\gamma =3\times {10}^3 \mathrm{A}$, is located at $\left(x_{\mathrm{m}},y_{\mathrm{m}}\right)=\left(0.01\mathrm{m},0.003\mathrm{m}\right)$. We observe two vortices forming in the vicinity of the magnetic source. As $\phi$ increases both vortices increase in size and in strength.

Figure 7a shows the variation of TAVC along the upper wall of microchannel versus volume fraction $\phi \left(0.03, 0.06, 0.09\right)$. We observe that as the volume fraction $\phi$ increases, there is an increase in the deposition of MNPs on the upper wall. Figure 7a shows a steep decrease of TAVC close to the magnetic source.

Figure 7b shows the TAWSS along the upper wall of the microchannel versus volume fraction $\phi$. We observe elevated values of shear stress close to the magnetic source.

4.1.2. Effect of the Magnetic Field Location

We examined the effect of magnetic field location over the flow field and the distribution of MNPs non-dimensional concentration. The magnetic field source is located at three different positions above the midpoint of the upper wall $b=0.001,$ $0.002,$ and $0.003 \mathrm{m}$. The volume fraction (φ) of MNPs is set $0.03$ at the inlet.

Figure 8 depicts the stream functions contours at the peak systolic velocity for t = 0.8 s. Strong recirculation zones, upstream and downstream of the location of the magnetic source are created for the two closest locations of the magnetic source. As the b increases, we observe that the magnitude of the disturbances decreases. For $b=0.003 \mathrm{m}$ the effect of the magnetic field on the flow field is limited.

The profile of the TAVC along the upper wall of microchannel is shown in Figure 9a, for three specific locations of the magnetic field source. It is observed that when the magnetic field source is located closed to the upper wall, $b=0.001 \mathrm{m}$, there is an increase of TAVC upstream of the magnetic field source. Similar behavior of the TAVC of MNPs is followed when the distance b is increased, $b=0.002 \mathrm{m}$ and $b=0.003 \mathrm{m}$, and thus the influence of the magnetic field is decreased. Figure 9b shows the effect of the magnetic field location on TAWSS. We observe that when the magnetic source is far from the upper wall of the channel, $b=0.003 \mathrm{m}$, the effect of the magnetic field on the TAWSS is negligible. When the distance $b$ is decreased, the TAWSS increases closed to the magnetic source. The highest TAWSS values are observed when the distance b is minimum and equal to $0.001 \mathrm{m}$.

4.1.3. Effect of the Magnetic Field Strength

Next, we examine the effect of the magnetic field strength and magnetic field gradients over the flow field and the distribution of MNPs non-dimensional concentration. The external, non-uniform magnetic field source is located at $\left(x_{\mathrm{m}},y_{\mathrm{m}}\right)=\left(0.01\mathrm{m},0.003\mathrm{m}\right)$ while the strength of the magnetic field takes values in the region (${10}^3 \mathrm{A}\le \gamma \le 7\times {10}^{3 }\mathrm{A}$). The volume fraction of magnetic MNPs $\phi$ is set $\phi =0.03$ at the inlet.

In Figure 10 we present the stream function contours for four values of the magnetic field strength γ at the peak systolic velocity when t = 0.8 s. The increase of the strength causes the formation of two recirculation zones (vortices) upstream and downstream the magnetic source. Moreover, the increase of γ increases the size and the strength of the vortices and thus causes further disturbances of the flow close to the magnetic source.

Figure 11a shows the variation of TAVC along the upper wall versus the intensity of the magnetic field γ for one period of the cardiac cycle. We observe that, when the strength of the magnetic field increases, the deposition of MNPs on the upper wall increase upstream the location of the magnetic field source. Moreover, the diagram shows a steep decrease of TAVC close to the magnetic source. The maximum TAVC value appears close to the magnetic source and is observed for the magnetic field strength $\gamma =5\times {10}^3 \mathrm{A}$. Figure 11b presents the TAWSS along the upper wall of the channel versus magnetic field strength. The strength of the magnetic field affects the TAWSS on the upper wall of the channel. Close to the magnetic source we observe high negative and positive levels of shear stresses. Higher levels of TAWSS can be noticed when the strength of the magnetic field increases.

4.2. MDT When MNPs Are Injected from the Upper Wall of the Microchannel

In this section, we investigate the MDT when MNPs are injected into the blood stream from the upper wall of microchannel. Figure 12 shows the flow domain and the boundary conditions.

We apply no-slip boundary conditions on the walls, and fully developed flow at the outlet ($\partial u/\partial x=v=0$). At the inlet, we apply the flow waveform used in [43]. We apply Neumann boundary conditions at the walls and the outlet for the non-dimensional concentration. The MNPs are steadily injected with velocity $v=0.03 {\mathrm{m}\mathrm{s}}^{-1}$ in the blood stream, having a volume fraction of φ = 0.03 at the inlet. The cross-section area $A$ of the injection site is $A=2.34\times {10}^{-4}{\mathrm{m}}^2$. The distance between the injection site and the inlet is $1.8\times {10}^{-3} \mathrm{m}$. The magnetic source is located at $\left(x_{\mathrm{m}},y_{\mathrm{m}}\right)=\left(0.01\mathrm{m},0.003\mathrm{m}\right)$, with magnetic field strength set to $\gamma =3\times {10}^3 \mathrm{A}$.

Figure 13 shows the distribution of the dimensionless concentration of MNPs in the vessel for various times during a cardiac cycle. When injection begins, it takes about $0.3 \mathrm{s}$ before the MNPs approach the targeted tissue. The concentration close to the targeted tissue increases during the cardiac cycle when the velocity decreases. As a consequence, MNPs create an obstacle in front of the targeted tissue, and this obstacle is removed at the systole time when the hydrodynamic forces overcome the magnetic forces.

Figure 14a shows the stream function contours computed at the peak systolic velocity ($t=0.8\mathrm{s}$) when MNPs enter the microchannel from the inlet, while Figure 14b when the MNPs are injected from the upper wall. In the former, two recirculation zones (vortices) are formed upstream and downstream to the magnetic source, while in the latter only one recirculation zone (vortex) is formed upstream the magnetic source at the upper wall.

Figure 15a shows the variation of TAVC along the upper wall for a time period of one cardiac cycle for the two ways by which the MNPs approach the targeted tissue. The values of TAVC are higher when the second way is used. Figure 15b presents the TAWSS along the upper wall of the microchannel for one period of the cardiac cycle for the two ways by which the MNPs approach the targeted tissue. We observe minor differences for the values of TAWSS for the two approaches.

4.3. MDT in a Microchannel with Stenosis

In this section, we investigate the biomagnetic fluid flow in a microchannel with a stenosis located in the middle of the upper wall. We set the volume fraction $\phi$ of MNPs, entering at the inlet, to $\phi =0.03$, the strength of the external magnetic field $\gamma ={10}^3 \mathrm{A}$ and the magnetic source in a distance $y_{Mag}=0.002 \mathrm{m}$ above the upper wall. We change the horizontal coordinate for the magnetic source location to $x_{Mag,1}=0.009 \mathrm{m}$ (in front of the stenosis), $x_{Mag,2}=0.01 \mathrm{m}$ (in the middle of the stenosis) and $x_{Mag,3}=0.011 \mathrm{m}$ (after the stenosis). We use a grid of approximately 58,000 nodes, that ensures a grid independent numerical solution. The time step was set to $dt={10}^{-5}$ s. The geometry and the boundary conditions of the problem are presented in Figure 16.

The stream function contours are presented in Figure 17 for three values of horizontal position of the magnetic source at the peak systolic velocity when t = 0.8 s. We see that two recirculation zones (vortices) are formed upstream and downstream of the stenosis. As the magnetic source is shifted from the left side to the right side of the stenosis, the first vortex is minimized, and the second vortex is enlarged.

The profile of TAVC along the upper wall for one period of cardiac cycle is shown in Figure 18a. As the magnetic source is shifted towards the right, we observe variations of the TAVC which are located before and after the stenosis. Figure 18b presents the TAWSS along the upper wall of the microchannel for one period of cardiac cycle. The main differences for the three different positions of the magnetic field are found before the stenosis. In the middle of the upper wall of the channel where the stenosis is formed, we observe high levels of shear stresses for the three cases.

4.4. MDT in a Microchannel with Bifurcation

The branching angle and the diameter ratio in epicardial coronary artery bifurcations are two important factors of atherogenesis [46]. In this section, we examine the effect of different branching angle ($30°, 45°, 60°$) of the bifurcation vessel on the MNPs concentration and the flow field. We present numerical results for the MDT biomagnetic (blood and dispersed MNPs) fluid flow in a microchannel (vessel) with bifurcation. Details on the geometry and the boundary conditions are in Figure 19. We set the MNPs volume fraction to $\phi =0.03$, and the intensity of the external magnetic field to $\gamma ={10}^3 \mathrm{A}$. The magnetic field is located at $\left(x_{Mag},x_{Mag})=(0.00125 \mathrm{m},0.0002 \mathrm{m}\right)$. We simulate two cardiac cycles. We use a grid of 120,000 nodes, which ensures a grid independent numerical solution. The time step was set to $dt=5\times {10}^{-6}$ s.

Figure 20 shows the stream function contours for three bifurcation angles ($30°$, $45°$, $60°$) at t = 1.6 s. As the bifurcation angle decreases, a vortex is generated at the lower wall of the bifurcation.

Figure 21 shows the distribution of the dimensionless concentration of MNPs in the vessel at different time instances. At the upper branch wall of the bifurcation, the MNPs capture is increased, causing a delay for MNPs to wash out at the outlet of the upper branch. The opposite is observed at the lower branch of the bifurcation, where the intensity of the magnetic field is weaker.

Figure 22a shows the variation of the TAVC along the upper wall of the bifurcation for various bifurcation angles. We observe that as the bifurcation angle decreases, there is an increase of TAVC of the magnetic field source. Finally, Figure 22b shows the variation of the TAWSS along the upper wall of the three bifurcations. The bifurcation angle affects the TAWSS upstream of the magnetic source, while downstream of the magnetic source the TAWSS maintains their values to the same order of magnitude.

5. Conclusions

We presented a numerical model to investigate the MDT flow problem. We examined the concentration distribution of MNPs under the presence of an external non-uniform magnetic field. The flow model used in the MDT flow simulation was based on the principles of FHD and BFD. The bio-magnetic fluid (blood and dispersed ${\mathrm{Fe}}_3{\mathrm{O}}_4$-MNPs) was modeled as two-phase, non-Newtonian fluid and the flow as pulsatile and laminar. The following key points were derived:

• MNPs loaded with medicinal drugs can be carried to a tissue target in the human body under the applied magnetic field.
• Recirculation zones are formed upstream and downstream of the targeting area (which is near to the location of the magnetic field), as the volume fraction of MNPs increases.
• The capture of MNPs on the tissue target increases with the volume fraction of MNPs and high values of the TAWSS are presented.
• When the magnetic source is placed close to the tissue target, recirculation zones are formed upstream and downstream of the targeting area.
• The increase of the strength of magnetic field causes the increase of TAWSS close to the tissue target.
• The MDT flow case in a microchannel with bifurcation shows that the decrease of the bifurcation angle results to the formation of one recirculation zone at the lower branch of the bifurcation, while at the upper branch, MNPs capture is larger.