1 Introduction

The cancer disease is one of the leading causes of death around the world,and the tumor-targeted drug delivery is one of the major areas in its treatment. The guiding magnetic iron oxide nano-particle with the help of an external magnetic field to its target is the principle behind the development of super paramagnetic iron oxide nano-particles (SPIONs) as novel drug delivery vehicles. The accurate prediction of local dynamical behavior of discrete particles released in the fluid flow is an important key to better understanding and optimization of many processes in the biomedical application,e.g.,deposition of hazardous particles in the human respiratory or cardiovascular system. There are many researches related to the magnetic drug targeting (MDT) technique for drug delivery in the literature. The localized medical drug delivery enables a significant local increase of the medical drug in regions affected by the disease and leads to a significant reduction of the always present negative side effects of aggressive medical treatments. Experimental studies on animals and preclinical studies on human patients demonstrated the potentials of this approach^[1]. Figueroa et al.^[2] developed a new method to simulate the blood flow in three-dimensional deformable models of arteries. They presented here the mathematical formulation of the method and discussed issues related to the fluid-solid coupling. The method couples the equations of the deformation of the vessel wall at the variational level as a boundary condition (BC) for the fluid domain. Torii et al.^[3] investigated the fluid-structure interaction of blood flow and cerebral aneurysm. They used pulsatile BCs based on a physiological flow velocity waveform and investigated the relationship between the hemodynamic forces and vascular morphology for different arteries and aneurysms. Bin et al.^[4] adopted the discrete trajectory model to simulate particle tracks with the Eulerian method for solving the continuous fluid flow. The results showed that the particle deposition and the concentration are mainly affected by the ventilation conditions. For the same particle properties,a displacement ventilated room had a lower particle deposition rate and a larger escaped particle mass than the mixing one,while the average particle concentration of the displacement case was higher than the mixing case. Haverkort et al.^[5] studied the magnetic particle motion in a Poiseuille flow. Exact analytical solutions were derived for the particle motion under the influence of a constant magnetization force and a force decaying as a function of the source distance,e.g.,due to a current carrying wire or a magnetized cylinder. They observed that the orientation of the vessels with respect to the magnetic force crucially affected deposition rates of particles leading to a heterogeneous particle distribution. Hournkumnuard and Chantrapornchai^[6] studied the concentration dynamics of weakly magnetic nano-particles dispersed in a fluid medium during the process of high gradient magnetic separation by using a computational approach. Alexiou et al.^[7] studied ferromagnetic stents for implant assisted-magnetic drug targeting in vitro. Their results showed that the fluid velocity,particle concentration,magnetic field strength,and stent material were important for capturing the magnetic drug carrier particle (MDCP) surrogates. Nacev et al.^[8] studied the behavior of ferromagnetic nano-particles in blood vessels under applied magnetic fields. They carried out a detailed analysis to better understand and quantify the behavior of magnetizable particles in vivo. They found that there are three types of behaviors,i.e.,the velocity dominated,the magnetic dominated,and the boundary-layer formation uniquely identified by three essential non-dimensional numbers (the magnetic-Richardson,mass Peclet,and Renkin numbers). A comprehensive mathematical model for simulations of blood flow under the presence of strong non-uniform magnetic fields was used in Ref. ^[9]. The time-dependency of the wall-shear-stress for different stenosis growth rates and the effects of the imposed strong non-uniform magnetic fields on the blood flow pattern were investigated by Kenjeres^[9]. Haverkort et al.^[10] numerically simulated the magnetic particle motion in arterial flows. Firstly,the blood flow and the particle motion in a straight cylindrical channel and a 90$^\circ$ bended tube were analyzed in their study,and then they performed an unsteady simulation in the more complex geometries of a left coronary artery and a carotid artery. Cherry and Eaton^[11] studied the success of magnetic particle trapping in straight vessels using a simulation that models all nontrivial forces on the particles. The data showed that the blood flow can be highly disturbed by the injection of magnetic particles. The success of particle retention varies widely and depends on the concentration,the volume,and the aspect ratio of the injected particle cluster. Kenjeres and Righolt^[12] simulated the magnetic capturing of drug carriers in the brain vascular system. They demonstrated that magnetically targeted drug delivery significantly increased the particle capturing efficiency in the target regions,and the application of the magnetic drug delivery technique leads to a significant increase in the captured medical drug within targeted regions. Larimi et al.^[13] simulated guiding magnetic iron oxide nano-particles with the help of an external magnetic field to its target in a two-dimensional bifurcation vessel. They studied the behavior of the blood flow and magnetic particles (Fe$_{3}$O$_{4})$ influenced by the external magnetic field. Also,the effects of the high magnetization number on the blood,flow in two cases,oxygenated blood,and deoxygenated blood were investigated. They illustrated that,by use of an external magnetic field,particles can be delivered to a target region,and with increasing the magnetic gradient,the volume fraction of particles delivered to the target region,is increased.

In this paper,the magnetic particle behavior in an unsteady flow,considering the effect of magnetic field in a bifurcation vessel is studied. The performances of drug delivery and blood flow behavior for two cases of diabetic patients and healthy persons under treatment are compared. To the best knowledge of the authors,no other work has studied the effect about the drug targeting before. The Euler-Lagrangian method is used for tracking particles and investigating the particle volume fraction in the target position besides the effect of external magnetic field position. All the simulations are performed by the open source CFD software,i.e.,Open-FOAM. The blood flow in the vessel is considered to be non-Newtonian,incompressible,and laminar. The forces that have effects on the particles are Brownian,drag,buoyancy and magnetic forces. The partial magnetic field is used in vessel walls,and its effect is compared for different volume fractions of particles.

2 Governing equations 2.1 Fluid phase

For the non-Newtonian,incompressible,viscous fluid flow of the density $\rho $ (blood),the governing equations are as follows.

The momentum equation is considered as

(1)

where $\frac{\mathrm D}{\mathrm Dt}$ is the material derivative,$U_{\mathrm f}$ is the fluid velocity ($u,v,w)$,$P$ is the pressure,and $\rho _{\mathrm f}$ is the density of fluid. In order to model blood's complex rheological properties,the macroscopic non-Newtonian Carreau model^[17] is investigated,

(2)

where $\eta _{0}$ and $\eta_{\infty }$ are the zero and infinite shear viscosities,respectively,$\lambda $ is the characteristic relaxation time,and $n$ is the flow index. The four parameters occurring in the Carreau model may be determined by numerical fitting of experimental data for a healthy person (diabetes-controlled) from Kostovaa et al.^[15],

(3)

2.2 Particle phase

In the Lagrangian frame of reference,the equation of motion of a nano-particle is given by

(4)

where ${U}_\mathrm p $ is the particle velocity,$ {F}$ is the force,including contributions from the drag,magnetic,virtual mass and Brownian force,

(5)

2.2.1 Magnetic force

Electromagnetic fields are classically described by Maxwell's equations. For the special case of magneto-static equations that are appropriate for stationary,or slowly varying magnetic fields^[17],

(6)

(7)

(8)

where $\mu _{0}$ is the magnetic permeability of vacuum,$H$ is the magnetic intensity,$B$ is the magnetic field,$J$ is the current density,and $\chi $ is a dimensionless parameter describing the magnetic susceptibility. A single ferromagnetic particle in a magnetic field will experience a force that depends on the magnetic field and its gradient^{[5, 9]},

(9)

where $M$ is the magnetization of nano-particles (magnetic dipole moment per unit volume). The magnetization of such particles can often be assumed to be approximately proportional to the applied magnetic field. Above a certain magnetic field strength,the magnetization saturates to a constant value $M_\mathrm {sat}$,leading to the following model:

(10)

where $\widehat {H}=\frac{{H}}{|H|}$. Therefore,the magnetic force can be obtained by^{[5, 9]}

(11)

As a property of a vector field,we have

(12)

Then,for a curl-free magnetic field,it can be concluded

(13)

Therefore,

(14)

By substituting Eq. (14) into Eq. (11) ,we have

(15)

where $V$ and $a$ are the volume and the radius of the magnetic particle,respectively.

2.2.2 Drag force

When the magnetic force of Eq. (15) is applied to a particle,it will make the particle accelerate in the direction of the force until it reaches to an equilibrium velocity,$ U _\mathrm R $,relative to the surrounding fluid,

(16)

The opposing Stokes drag force on a spherical particle is given by

(17)

where $\rho _\mathrm f $,$\rho_\mathrm p $,$m_\mathrm p $,and $d_\mathrm p $ are the density of the fluid,the density,the mass,and the mean diameter of particles,respectively. $ U_\mathrm f $ is the fluid velocity,and $ U_\mathrm p $ is the particle velocity in each time step. $C_\mathrm D $ is the drag coefficient which can be obtained from the following equation^[18]:

(18)

where the Reynolds number of particles is defined as^[18](19)

Here,$\mu _\mathrm f $ is the fluid viscosity. With combination of Eqs. (17) and (18) with (19) ,the drag force acting on a particle will be equal to

(20)

2.2.3 Buoyancy force

The following equation is used to calculate the force on a particle corresponding to the buoyancy force^[18]:

(21)

where $g$ is the gravity acceleration.

2.2.4 Brownian motion

The components of the Brownian force are modeled as a Gaussian white noise process^[18] with the spectral intensity $S_{n,ij}$,

(22)

where $\delta _{ij} $ is the Kronecker delta operator,and

(23)

where $T$ is the absolute temperature of the fluid,$v$ is the kinematic viscosity,and $C_{\mathrm c}$ is known as the Cunningham correction factor. $K_{\mathrm B}T$ is known as the thermal energy,in which $K_\mathrm B $ is the Boltzmann constant. Amplitudes of the Brownian force components are of the form,

(24)

Here,$\zeta $ is the Gaussian random number with the zero mean and unit variance. With combination of Eqs. (22) and (23) with Eq. (24) ,the Brownian force acting on a particle will be equal to

(25)

2.3 Particle motion

By considering Eqs. (15) ,(19) ,(20) ,(21) ,and (25) and this point that the applied magnetic field is not sufficient to saturate the nano-particles,Eq. (5) can be written as

(26)

where

(27)

2.4 Particle-wall interactions

The interaction of a particle with a solid surface can lead to adherence and sedimentation (accumulation) and the severe effect of surface separation (erosion) and chemical accumulation (corrosion). In fact,the mechanical interaction is adherence. The adherence model considers the elastic characteristics of the particle and the surface under dry conditions until the particle sticks. The next step would be known whether the particle gets separated from the surface or remains on it which is based on the critical momentum theory. This step is called separation and is the dynamical interaction of the fluid. Dahneke^[19] experimentally studied the effect of the velocity of the particle collision on the reflected velocity for the spherical particles. When these velocities reduce,the importance of adhesive force increases,and the reflected velocities drop significantly. In other words,when the normal velocity of the particle reduces,the reflected velocity decreases and eventually approaches to a point where no return occurs,and all particles are precipitated. This velocity in which the particle sedimentation takes place is called the capture velocity. Brach and Dunn^[20] presented a relationship for calculating the capture velocity of a particle using a semi empirical model,

(28)

where $E$ is named as the compound Young module and is determined based on the Young module of a particle. The normal velocity of a particle $v_{\mathrm n}$ is compared to the capture velocity. If the normal velocity is lower than the capture one,the particle sticks to the surface,otherwise,it has got separated from it,i.e.,

When the particle is separated,it keeps its route until it leaves the area or recollides the surface. The compound Young module is defined as follows:

(29)

where $E_{\mathrm s}$ is the surface Young module which is about 0.45 MPa for a vessel wall of 3-4 millimeters thick,$v_{\mathrm s}$ is the surface Poisson's ratio considered as 0.45 (Ilic et al.^[21]),$E_{\mathrm p}$ is the Young module of the particles which is 250 GPa for Fe$_{3}$O$_{4}$ particle,and $v_{\mathrm p}$ is the particle Poisson's ratio which is 0.3 for Fe$_{3}$O$_{4}$ particle (Ouglova et al.^[22]). $\rho _{\mathrm p}$ and $d_{\mathrm p}$ are the particle density and diameter,respectively. The Young module has a significant effect on the particle capture velocity. When the Young module increases,the capture velocity decreases (Ilic et al.^[21]). Decomposing the particle velocity into two tangential and normal directions after hitting the wall would lead to the following equation:

(30)

where $ U _\mathrm p $ is the particle velocity after collision,$ U _\mathrm {pn} $ is the velocity normal to the hitting direction after collision，and $ U _\mathrm {pt} $ is the tangential one. The magnitudes of the normal and tangential velocities after hitting the wall could be obtained through the following correlation:

(31)

(32)

where $ u _\mathrm {pn} $ is the normal velocity before collision,and $ u _\mathrm {pt} $ is the tangential one before collision. Here,$\varsigma $ and $\psi $ are the wall restitution coefficient and the wall friction coefficient,respectively,and $\varsigma$,$ \psi \in [0,1] $. To calculate the particle collision velocity after colliding to the vessel wall,the restitution and friction coefficients have been considered to be 0.5 and 0.2,respectively^[22]. The values of some parameters used in this paper are introduced in Table 1.

3 Geometry and mesh independency

The geometry that is considered in this study is a bifurcation vessel with an inlet (3 mm diameter) and two outlets. BCs and all sizes are shown in Fig. 1.

In order to obtain a mesh in the most optimum computational cost,four different meshes are used with 400 000,600 000,800 000 and 1 000 000 grids,respectively,with the Reynolds number of 700 and the magnetic field strength of 4 T. In Fig. 2,the shear stress is plotted along the upper wall from the inlet,up to the outlet 1 (along the green line which is shown on the upper wall in Fig. 1) . According to Fig. 2,the shear stress magnitudes for the first and second meshes have large differences compared with other results. However,the third and the fourth ones are almost completely in agreement with each other and have no diversity. \vspace*{2mm}

4 Results and discussion

For all simulations,the particle diameter is considered to be 0.5 μm,the magnetic susceptibility of nano-magnetic particle is considered to be 3,and the particles density is equal to 4 600 kg/m$^{3}$ (see Ref. ^[5]). Figure 3 shows the nano-particles distribution under influence of the magnetic field. Particles are released in all time steps in order to obtain equal conditions in all cases and to observe the effect of the releasing time on the particles. Therefore,the numbers of particles in all conditions are the same,and results are only affected by BCs and the inlet conditions. Particles with 0.5 μm diameter are released in the inlet with a uniform distribution throughout the whole cross section. The Reynolds number also is considered to be 500. In Fig. 3(b),the volume fraction of the nano-particles is shown in four cross sections,respectively.

In Section (I00) ,which is close to the inlet and the particle injection region,the nano-particles distribution is uniform and almost similar throughout the whole cross section. According to the fact that the particles with the initial velocity of zero have been located on the inlet and in this section,the nano-particles have not yet approached the fluid velocity,the particles accumulation in this region is higher than other sections,as a result,the volume fraction in this region is higher than other sections. Entering the region of the field effect which could be observed in the figure in Section (I01) ,a fraction of nano-particles is moved towards the upper wall closer to the field and gets further from the lower wall.

In Section (I10) which is the place of the maximum strength of the field,the distance between nano-particles is larger than that in the lower section due to the increase in the magnetic field strength which verifies the non-uniform distribution of the volume fraction of the nano-particles in the cross section. This non-uniform distribution in the region of applying field approves the motion of nano-particles towards the upper branch and the increase in the deposition in this region. In the last section (I11) ,which is selected almost at the end of the field affected area,through getting out the field affected region,the nano-particles are returned to the lower wall,and the nano-particles distribution would change into a uniform state in the upper branch. In Fig. 4,the effect of the field strength on the way drug carrying nano-particles move through the branches is investigated. For this purpose,the upper branch cells are selected,as shown in Fig. 4. Cells have been selected the way that the precipitated particles in the location of maximum strength would not interfere in the calculations of volume fraction. The maximum strength of the field has been selected in four cases (1 T,2 T,3 T,and\linebreak 4 T). As shown in Fig. 4,the contour of magnetic field distribution in the upper part and the effect of magnetic field on the average volume fraction are illustrated. All conditions of particles injection,including the nano-particles diameters,the number of injected particles in each time range,the initial velocity,and all boundary and initial conditions related to the fluid solution,are considered to be the same for all the four cases. The calculations have been done after a temporal cycle,i.e.,when the first particle gets out the geometry outlet. The nano-particles volume fraction in the upper branch will increase by increasing the field strength,and in fields of 1 T,2 T,3 T,and 4 T,the average volume fraction of nano-particles carrying drug in the upper branch is increased by 30,83,105,and 140 percent,respectively,compared to the case without field.

Figure 5 shows the deposition of magnetic nano-particles in the wall affected by the magnetic field. Except for the magnetic field strength,all the conditions are the same for four geometries,and the results are investigated after a cycle,i.e.,the first particle gets out of the geometry outlet. The results show that,in the region affected by the magnetic field,deposition of particles increases by increasing the field strength. This characteristic could be used to increase the nano-particles deposition in the tissue or the targeted location in the medical treatments,such as the tumor location or regenerative medicine delivery to the region where the bone has faced fracture.

The effect of the flow Reynolds number on the targeted drug delivery efficiency is shown in Fig. 6. In order to investigate this effect,the constant magnetic field of 2 T is used. The investigations are assessed in four different Reynolds numbers,and in all four cases,the diameters of magnetic nano-particles are considered to be 0.5 μm,and the number of particles and their entering time are accounted similar. The location and strength of the magnetic field are considered to be equal and similar in all four cases as shown in Fig. 7(a). The results show that through a rise in the Reynolds number,the average volume fraction of nano-particles reduces in a branch where the magnetic field source is located. In lower Reynolds numbers,the dominant force located in the region with the field strength concentration is a magnetic force,which could increase the movement of nano-particles towards the upper branch and increase the nano-particles volume fraction in this region. By increasing the Reynolds number,the inertial and drag forces of the particle become dominant over the magnetic force,and the nano-particles fraction in the upper branch reduces to the point that the volume fraction of the magnetic nano-particles in the upper and lower branch becomes equal. Therefore,it could be realized that in the vessels with the higher Reynolds number,the targeted drug delivery method will have lower efficiency.

In addition to the changes in the nano-particles volume fraction in the upper branch,the Reynolds number can have influence on the particles deposition in the wall affected by the magnetic field. In order to investigate this effect,an area of the upper wall affected by the magnetic field,is selected (see Fig. 7(a)).

As shown in the diagram,an increase in the Reynolds number could lead to a reduction in the average nano-particle volume fraction deposition in the specific region. The distribution of the particles color in the specified regions could also determine the magnetic field strength applied on the particle entering to it. One of the most common diseases in human in the world is diabetes. Diabetes is a type of disease in which the body fails to regulate the amount of glucose necessary for the body. Diabetes does not allow the body to produce or properly use insulin. The insulin is a hormone that is needed to convert sugar,starches and other food into energy needed for daily life. Diabetes has widespread fallout,with more than 16 million people affected by it in the United States alone (by the national diabetes statistics report,2014) . Even more concerning is the fact that a wide number of people affected by the disease remain unidentified or are unaware that they have it. Diabetes is one of the major metabolic disorders that affect the human body. Therefore,in this part,the effect of magnetic field on particle deposition and comparison of MDT efficiency in a patient with diabetes disease and a healthy human is investigated. For this purpose,as mentioned in the previous section,the results of Kostovaet al.^[15],for two samples are considered,i.e.,a diabetics patient (diabetes-diabetes) and a healthy control subjects (diabetes-controlled). The relationship between strain rate and blood viscosity is estimated by the Carreau model^[14]. The parameters $n$,$\eta _{\infty }$,$\eta _{0}$,and $\lambda $ are dependent upon the fluid. $\lambda $ is the time constant,$n$ is a measure of the deviation of the Newtonian fluid,and $\eta _{\infty }$ and $\eta _{0}$ are,respectively,the zero- and infinite-shear viscosities. The calculated parameters of the Carreau model based on the data from Kostovaa et al.^[15] for both samples are presented in Table 2.

For comparison of the simulation results in the two samples,all particles properties,for example,the time,the number of particles injected in each time step,the size and the injection position of particles and the Reynolds number ($Re=500$),are considered to be similar in both cases to make the results independence of these parameters.

Figure 8 shows the streamlines and nano-particles deposition in the area affected by the magnetic field. In Figs. 8(a) and 8(b),the volume fraction of nano-particles in a cross section near the affected magnetic area (Section I10 of Fig. 3) ,is shown for both samples. It is observed that\linebreak

the distribution of nano-particles (volume fraction) in both samples is almost the same but concentration and compactness of particles for the diabetes-controlled sample in the central position are more in comparison with the diabetes-diabetes sample. This can be due to more deposition of nano-particles in the wall affected by the magnetic field for the diabetes-controlled sample. Also,in Figs. 8(d)-8(e),the deposition of nano-paricles in the same legend is shown. The deposition on the diabetes patient is more than the health sample. Figure 9 shows the deposition of nano-particles in a vessel wall affected by the magnetic field (see Fig. 4(a)). As it is obvious that the maximum of the nano-particles deposition of the diabetes-diabetes patient is more than that of diabetes-controlled (about 44 percent higher),it may be because of lesser fluid momentum in the diabetes-diabetes sample. Also,as shown in Fig. 10,the pressure magnitude of the diabetes-diabetes sample is more than that of the diabetes-controlled patient which can be another reason for more depositions of the diabetes-diabetes sample.

5 Conclusion

In this paper,the MDT efficiency,influenced by various variables,is investigated. In the case of considered bifurcated vessel,the maximum pressure occurs at the divider wall. Furthermore,the maximum wall shear stress occurs at the beginning of the bifurcation due to the higher velocity gradient in this region. The bifurcation angle reveals that the maximum pressure is higher for the larger bifurcation angles. In the range of magnetic fields and diameters that are considered for drug targeting goals in this paper,it is illustrated that the fluid behavior does not be affected by the magnetic field,due to lower magnetism of blood in comparison with nano-particles. Incorporating the magnetic field at the bifurcation site results in the higher volume fraction of nano-particles at the upper branch and the increase of nano-particles settling at the region of maximum magnetic field strength. Simulation results reveal that the guided drug delivery efficiency is directly proportional to the magnetic field strength. However,the efficiency decreases for higher Reynolds numbers. The higher magnetic field strength helps particles settle at the site of exposure,but it is important to remember that the applicable magnetic field strength is limited by medical applications,sensitivity of exposure site to the magnetic field,patient age and condition and location of exposure with respect to the heart. Also,the efficiency of MDT is compared for two samples,and results reveal that the MDT for a patient with the diabetes disease is better than a healthy person.