A two-phase, two-way coupled model of targeted magnetic drug delivery for small Reynolds numbers

Abstract Targeted magnetic drug delivery (TMDD) is a promising approach relevant to multimodal cancer therapy. A majority of the TMDD models represent a mixture of blood and nanoparticles as a one-phase solution. However, in many cases it is an oversimplistic assumption. The existing two-phase models are usually one-way coupled, i.e. the blood flow has an impact on the MN flow. However, the inverse impact of the MNs on the the blood is not included. To eliminate these drawbacks the model is governed by two–way coupled momentum and temperature equations for the blood flow and the MN. The numerical procedure invokes the stream function–vorticity formulation and an efficient numerical procedure. The model, validated by experimental results, has been applied to analyze the formation of vortices generated by a combination of the magnetic force (MF) and the drag force (DF). The model also simulates the zones of TMDD and the corresponding changes in the vorticity. Finally, the model includes the impact of the size and concentration of the MNs on the temperature of the blood. These important scenarios cannot be analyzed by the earlier models.


Introduction
During conventional delivery, a drug travels through the entire circulatory system of the patient. However, when a high concentration of the drug is required to treat a particular area of the patient's body, e.g. chemotherapy, the procedure may have undesirable side effects. The most sensitive areas are the bone marrow, hair, skin, gastrointestinal organs, and the immune system. Nevertheless, recent studies have shown that TMDD is a powerful alternative to conventional drug delivery (Haverkort et al., 2009). It solves the main problem of chemotherapy, which is the inability to deliver the required drugs to the affected areas, reducing the toxic effects. A recent survey by Mamun et al. (2020) observes that in TMDD 'drug concentrations at the tumor are significantly higher, compared to drugs delivered by systemic delivery methods.' Among the blood vessels, TMDD through arteries is particularly useful for cancer therapies. The hepatic artery is used for the treatment of liver tumors (Jeon et al., 2016), the femoral artery is used for treating hind limb tumors (Alexiou et al., 2011), and the carotid artery is used for brain tumor treatment (Chertok et al., 2010).
The drug is placed on the MNs, injected near the tumor, and delivered using an MF created by an external CONTACT Stanislav S. Makhanov makhanov@siit.tu.ac.th or even implanted magnet. Although the approach reduces the side effects of anti-cancer treatment, it is a challenging and difficult process (Fernández et al., 2018;Lee et al., 2017;Luong et al., 2017;Rosiere et al., 2018;Sun et al., 2018a;Sun et al., 2018b;Wei et al., 2017). A weak MF may not achieve the required targeting, whereas a strong MF may cause leakage from the vessels through the internal organs. In addition, the MF can generate an undesirable vorticity which slows down the blood flow down to dangerous levels (Goodwin et al., 2001;Widder et al., 1981). The particle size and the intensity of the MF are two important delivery parameters. A strong MF combined with large particles may clog the blood vessels while a weak MF may not be able to deliver the drug. Widder et al. (1981) conducted an experimental study of Yoshida sarcoma tumors in rats. Albumin microspheres (placebo microspheres and microspheres with 0.5 mg/kg doxorubicin) coated by magnetite Fe 3 O 4 (100-200 nm) and doxorubicin (0.5-5 mg/kg) were injected. The external magnet (0.55 T) was held around the tumor for 30 min. The tumors decreased on average by 83%, and no deaths occurred for treated tumors. For placebo-treated animals, there were 80% deaths. The death rate for the third group treated by the same MNs without the external magnet was 100%. Lübbe et al. (1996a) and Lübbe et al. (1996b) conducted the first phase clinical trial of TMDD to deliver the anti-cancer drug 4'-epidoxorubicin to advanced solid tumors for 14 patients. The MNs were 50-150 nm in diameter and coated with anhydroglucose polymer on the surface where the drug was adsorbed. External magnets arranged on the skin close to the tumor provided an MF of 0.5-0.8 T. The results of the treatment were considered positive for 50% of the patients. To improve the results, a stronger MF and larger MNs (up to 1 μm in diameter) were suggested. The discussion of these trials is found in (2001) and (1999). Thereafter, Alexiou et al. (2000) and Alexiou et al. (2001) performed another extensive study on squamous cell carcinoma tumors in rabbits. They used MNs (50-100 nm) bound with anticancer agents. The MF was increased to 1.0-1.7 T. Accumulation of the chemotherapeutic agents with a permanent remission of the tumors was detected by visual inspection with MRI combined with histology. Goodwin et al. (1999) showed experimentally a high concentration of iron/carbon particles of 0.5-5 μm diameter in the targeted areas in animal experiments. However, in many cases, large particles block the vessels, causing a hemorrhage. Jain et al. (2003) recommended MNs of about 400-600 nm.
To date, TMDD has been limited either to using external permanent magnets to target shallow tumors (Krukemeyer et al., 2012;Liu et al., 2010;Pouponneau et al., 2011;Raut et al., 2010;Wilson et al., 2004) or to implanting magnetic wires, seeds or stents to reach a deep tumor (Cregg et al., 2012;Fernández-Pacheco et al., 2007;Forbes et al., 2008). Magnetic implants are promising for treating bone cancer. They can also be embedded in fatty tissue to treat obesity (Saatchi et al., 2017) and in the inner ear to treat deafness (Le et al., 2017). However, they are not suitable for every patient and every clinical condition. Potentially harmful procedures are often required for placing such magnets (Dobson, 2006). Therefore, application of these procedures remains limited (Cao et al., 2011;Hayden & Häfeli, 2006;Shapiro, 2009). Lübbe et al. (1996a), Nacev et al. (2011a, 2011b and Nacev et al. (2015) use strong permanent magnets (20×40 cm) to reach deep tumors. A magnet is placed outside the patient's body up to 15 cm from the tumor, to create an appropriate MF. Several authors have proposed that a permanent magnet shaped like a solid cylinder or rectangle can produce stronger MFs. Superparamagnetic particles constitute another solution (Kayal et al., 2011;Nacev et al., 2015;Rukshin et al., 2017;Sharifi et al., 2019). Nacev et al. (2011a), Lübbe and Bergemann (2005), Alexiou et al. (2000), and Alexiou et al. (2001) study TMDD using an advection-diffusion model verified by an experimental data. They predict TMDD with an MF of about 2 T using large blood vessels is possible when the distance to the tumor is less than 20 cm. The human clinical trials by Lübbe et al. (1996a), Lübbe et al. (2001), and Goodwin et al. (2001) detect the accumulation of MNs by visual inspection with MRI (magnetic resonance imaging) combined with histology in vivo. Unfortunately, the quantitative evaluation of the concentration of MNs was not performed since the MRI does not have the required resolution. The histology was performed, but the velocity of the blood flow was not measured. Testing the required MF per nanoparticle relative to the Stokes force was performed in vitro and in vivo experimentally by Widder et al. (1981), Ganguly et al. (2005), Xu et al. (2005), and Alexiou et al. (2000). In vitro studies (Ganguly et al., 2005;Xu et al., 2005) confirm that the MNs are captured even when the centerline DF exceeds the MF. In-vivo studies (Alexiou et al., 2000;Widder et al., 1981) confirm that when the DF exceeds the MF, TMDD is still possible. Nacev et al. (2011a) argue that a simple comparison of the DF and MF is not sufficient to derive a definite conclusion about the possibility of TMDD. Their example is an MF of 0.5 T acting on the MNs (250 nm) and a DF 7 times higher than the MF. Similar non-trivial experimental results by Alexiou et al. (2000) and Alexiou et al. (2001) imply that the TMDD requires numerical models capable of simulating the impact of the DF, the volume fraction, the size of the MN, and the temperature. Moreover, the conventional diffusion-convection approach (Grief & Richardson, 2005;Nacev et al., 2011a) seems to be oversimplistic.
An early work of Sud and Sekhon (1989) uses an analytical model for numerical analysis of the interaction between the MF and the blood flow in a multi-branching arterial system. Their results demonstrate that the rate of blood flow through the system was reduced by the MF. Kinouchi et al. (1996) include the Lorentz force into the Navier-Stokes equations. Their finite-element solution shows a reduction in the velocity of the blood flow by 5-10% under a strong MF.
Another early BFD (Biomagnetic Fluid Dynamics) model was introduced by Haik et al. (1996) and Haik et al. (1999). The model, including Ferro Hydro Dynamics (FHD) and Magnetic Hydro Dynamics (MHD), is solved numerically and verified by experiments. The fluid exhibits magnetization and vortices generated by the external MF. The blood flow has been significantly slowed down by an MF. Loukopoulos and Tzirtzilakis (2004) present the viscous, steady-state two-dimensional, incompressible, laminar BMF of a flow between two parallel plates. The FHD and magnetization are considered as functions of the temperature. The energy equation includes the magnetocaloric effect. In particular, the model simulates a vortex formed near the lower wall close to the operating magnetic device. Grief and Richardson (2005) propose a convection-diffusion model which includes interactions and collisions between the red blood cells. The model includes Brownian diffusion, shear diffusion, and particle convection. Bali and Awasthi (2007) present graphically the velocity resistance to the blood flow for different MFs. Tzirtzilakis and Loukopoulos (2005) analyze a one-phase FHD/MHD steady-state finite-difference numerical model with an external MF. The model generates two vortices rotating in the opposite directions. The vortices are generated even for relatively small magnetic numbers such as M nF ≈ 100. The MF is 8 T. The authors observe that a similar effect is achieved by an MF of about 1 T. Following this estimate, forthcoming numerical experiments consider an MF in the range of 0.5-1 T. Tzirtzilakis (2005) reports a 3D model consistent with the FHD and MHD, i.e. including magnetization and electrical conductivity of the blood. The application is a laminar, incompressible, three-dimensional, viscous flow of a Newtonian BMF. The blood flows in a rectangular duct under a spatially varying MF. Note that a sharp MF generated by a magnetic wire allows neglecting the Lorentz force. For the BMF under sharp gradients of the MF in large blood vessels, the dominant force is magnetization. Hence, the FDH model applies to this scenario. In particular, this assumption is valid for the blood flow in the vessels sized 10-20 cm and the MF exceeding 0.5 T (Haik et al., 1999). (These assumptions have been discussed in detail by Tzirtzilakis (2015)). Kenjereš (2008) and Kenjereš and Righolt (2012) present a numerical simulation of TMDD in a complex vascular system of the brain using a non-Newtonian model of blood. Their papers report simulations of TMDD in a complex brain vascular system. The MN flow is simulated by the Lagrangian tracking of the spherical double-layer particles. The particles are characterized by the diameter of the outer and inner (magnetic) cores. The momentum equation of the MNs is coupled with the momentum equation of the blood flow, whereas the momentum equation for the blood flow is independent. Kenjereš (2014) develops another version of the model to simulate the local deposition of a low-density lipoprotein. The model shows promising results, i.e. the TMDD is ten-fold, relative to conventional delivery.
Numerical simulation of TMDD in aerosols in the human upper and central respiratory systems under a non-uniform MF is presented by Kenjereš and Tjin (2017) and Rukshin et al. (2017). The particles are tracked in the Lagrangian frame. The model estimates the distribution of the MNs under a strong MF. Yue et al. (2012) develop a 3D model of a superparamagnetic cluster in a Poiseuille flow. Their model includes the Stokes DF, MF, and gravity. The numerical model of Habibi and Ghasemi (2011) shows a circulation in the region covered by the MF. Under certain conditions, the absorption of the 200-nm MNs decreases significantly compared to 2000-nm particles since the MF is too weak to overcome the DF. However, TMDD also depends on the concentration of the MNs and the strength of the magnet.
One-way coupled MN delivery has been used for stenosis aortic and vessel bifurcations by Larimi et al. (2014). The Lagrangian particle tracking is performed. The numerical results show that the MF increases the volume fraction of MNs in the target region. However if the blood flow is characterized by high Reynolds numbers, the efficiency of TMDD is low. Kandelousi and Ellahi (2015) apply the lattice-Boltzmann method. The results show that the MF affects the flow considerably. In particular, a backflow occurs near the MF region.
The validity of the above single-phase model is questioned by some authors. A recent survey (Said et al., 2021) suggests that the 'single-phase method suffers from a lack of precise simulations of thermophysical features related to nanofluids.' Therefore, Bose and Banerjee (2015) propose a two-phase one-way 2D model. The trajectories of the particles are obtained in the Lagrangian frame. The force balance includes the DF, the magnetic force, the buoyancy force, the Brownian force, and the thermophoretic force. The external MF affects the blood flow considerably, creating a strong recirculation zone near the location of the insert. Tzirtzilakis (2015) presents a BFD-FHD model of the blood flow in an aneurysm. Blood is considered to be an electrically non-conducting, homogeneous, non-isothermal Newtonian magnetic fluid. The paper presents an analysis of the effects of the MF on the blood that may or may not include MNs. The dominant force is the magnetization. The solution includes the stream function-vorticity formulation. The finite-difference solution is obtained on a curvilinear grid which is a variant of the algebraic grid generation (Gordon & Thiel, 1982), adapted to the shape of the aneurysm. Alshare and Tashtoush (2016) simulate magnetohemodynamics in stenosed arteries. Doubling the MF from 4 to 8 T increases the pressure drop by nearly 15%, but has a negligible impact on the wall shear stress. TMDD of pharmaceutical aerosols in the human upper and central respiratory is analyzed in (Kenjereš & Tjin, 2017). The airflow dynamics are based on the Eulerian approach, whereas the dynamics of the particulate phase are represented in the Lagrangian frame of reference. The deposition of the drug has been significantly improved by the MF. The most effective TMDD has been obtained for 5 × 10 −4 ≤ St ≤ 10 −1 with the MNs ranging from 0.3-5 μm.
An interesting single-phase model of the nanofluid flow designed to simulate the irrigation of a root canal in the human tooth is proposed by Ghalandari et al. (2019). An analytical single-phase model of the blood flow mixed with copper nanoparticles with a magnetic field is proposed in (Umadevi et al., 2021).
A recent two-phase model of TMDD for solid-liquid coupled BFD has been reported by Boutopoulos et al. (2020). The focus of the paper is an injection of MNs in different locations of the artery. The model is based on the Navier-Stokes equations for the blood flow and the advection-diffusion equation for the MNs. The experiments include different injecting scenarios. The model is one-way coupled, i.e. the blood flow affects the MNs via the corresponding convection-diffusion terms. However, the MNs do not have an impact on the blood flow. Therefore, the model is not capable of analyzing the possible impact of the size of the MNs on the characteristics of the blood flow. The model does not include the energy equations. Therefore, the analysis of the impact of the concentration and the size of the MNs on the temperature of the blood flow is excluded.
Our paper introduces a numerical, solid-liquid, two-phase Euler-Euler BFD model. The governing equations include continuity, momentum, and temperature equations for the blood flow and the MN flow in a 2D rectangular channel. The momentum and temperature equations of the blood flow and the MN flow are two-way coupled. The numerical procedure invokes the stream function-vorticity formulation and an efficient numerical method on a finite-difference grid.
The model has been analyzed with reference to single-phase models (Bianco et al., 2009;Tzirtzilakis, 2008) and a two-phase model (Boutopoulos et al., 2020). Note that the two-phase model is based on balancing the hydrodynamic and magnetic forces, where v is the velocity of the fluid and v mag is evaluated using the Stokes' drag law. In other words, the corresponding continuity, momentum, and energy equations of the solid phase are not included. In contrast, the proposed model includes the MNs having a two-way momentum exchange. This includes the reverse impact of the MNs on the blood flow. In particular, it is possible to analyze the impact of the size and the concentration of the MNs on the blood flow. It has been shown that under certain conditions, the geometric structure of the MN flow is different from that of the blood flow. Simulation of these effects is possible only with the proposed two-way coupled, two-phase model and its possible extensions. The formation, velocity, and size of the vortices are evaluated and discussed. We demonstrate that the size of the MNs has an impact on the temperature of the blood. Under certain conditions, the MF acting on the MNs increases the blood temperature up to 2.1°C, i.e. from 37 to 39.1°C. The size of this 'overheated' region depends on the size of the MNs. Note that this effect can not be simulated using the preceding models. Finally, a variety of engineering methodologies can be adapted to simulate the TMDD. Recent surveys (Afshari et al., 2021;Alsabery et al., 2020;Baghban et al., 2019;Granados-Ortiz et al., 2021) have shown a variety of algorithms designed to solve specific engineering problems. However, the different temporal and spatial scales of the medical and engineering applications must be considered.

Mathematical Formulation
Consider a viscous, unsteady, two-phase, two-way coupled flow consisting of the MNs and the blood in a 2D rectangular channel with lengthL and heighth ( Figure 1). The flow is assumed to be fully developed at the entrance of the channel. The upper and lower walls have a constant temperature. The magnetic source is a magnetic wire placed perpendicular to the (x,y)-plane at the point (a,b) near the lower wall. The governing equations of the two-phase flow are given with regard to the velocity V = (ū,v) and the temperatureTof the blood and the particles. The subscript p indicates the MN flow.
where the drag force F is given by: The solid phase, i.e. the MN flow is simulated by the momentum and energy equations, assuming that the impact of the dynamic viscosity is negligible. The MNs are assumed to be spheres with a uniform size and density. The fluid and the solid phases are electrically non-conductive. Therefore, the MF term is included in the momentum equations (Kenjereš & Tjin, 2017;Nacev et al., 2011). The temperature of the fluid and the MNs are affected by the MF and vise versa. Hence, the solid phase equations are given by: wherep is the pressure,ρ the density,μ the dynamic viscosity,c the heat capacity,k the thermal conductivity, C v is the concentration of MNs in the blood, and μ 0 = 4π × 10 −7 Hm -1 is the magnetic permeability of vacuum (Rosen, 2004). The additional termμ 0M ∇H in the momentum equation is the MF per unit volume, whereasμ 0T ∂M ∂T V · ∇H is the thermal power per unit volume from the magnetocaloric effect.
Following Tzirtzilakis (2008) and Tzirtzilakis (2015) we consider a current-carrying wire conductor positioned perpendicular to the (x, y) plane.H denotes the intensity of the MF given by: M is magnetization, approximated by (Matsuki et al., 1977)M K is an experimental constant, andT c is the Curie temperature.
The wire is at position(a, b) below the lower wall ( Figure 1). The strength of the MF is defined by the magnetic induction given byB(x,ȳ) = μ 0H ≡ μ 0 (H x ,H y ) (Tzirtzilakis, 2005), wherē and where γ is the intensity of the MF at the source.
Following (Bianco et al., 2009), the Stokes' resistance law for the small Re p and the sub micrometer particles is given byF where d is the diameter of the particle, f is the drag coefficient, and τ p is the particle response time. The Cunningham correction coefficient is given by where λ a , is the mean free path of the particle. Further, whereũ 0 (tilde) denotes the average, and index 0 refers to the inlet boundary. The drag coefficient is defined by: The energy equation for the spherical MNs is given (Ranz & Marshall, 1952).

Dimensionless Equations
The following dimensionless variables are introduced.
T 0 is the temperature of the MNs at the channel inlet,q 0 is the heat flux. All other notations are self-explanatory.

Grid Generation
The accuracy of the numerical algorithm in the regions with a strong MF requires an adaptive numerical grid. An adaptive grid can be considered as a discrete version of mapping ξ = ξ(x, y), η = η(x, y) from the physical region in the coordinates (ξ , η) to a computational region in the new coordinates(x, y), 0 ≤ x, y ≤ 1. We consider a particular case of a rectangular grid with stretching, i.e. ξ ≡ ξ(x), η ≡ η(y). The grid is defined by: where In the above expressions, x max =Lh is the dimensionless length of the channel, x 0 is the point where the grid clustering occurs, and τ is a parameter controlling the rate of clustering toward the ξ direction. It varies from zero (no stretching) to 1, to produce a dense grid near ξ = x 0 . Further, δ 1 and δ 2 are parameters to control stretching in the η direction. If δ 1 = 0, the mesh is adapted to η = 1, whereas, if δ 1 = 0.5, the mesh is adapted to η = 0 (Tzirtzilakis, 2008). Thus, the grid is adapted to the areas of large gradients of the MF. Figure 2 illustrates the techniques. Solving Equations (29) with regard to x(ξ ), y(η) yields: The partial derivatives in (23)-(28) are replaced as follows: where w x = dx dξ , w y = dy dη ., delete this line Some calculus yields The additional non-dimensional parameters are as follows

Numerical Method
The equations for the stream functions of the blood and the MN flow (34)-(35) are approximated by the standard finite-difference equations and solved by the block Gauss-Seidel method (M. Cervera et al., 1996). This is followed by a correction similar to the successive over-relaxation i.e. ψ k+1 GS is the numerical solution obtained by the block Gauss-Seidel method, ψ k GS is the solution by the GS at the previous iteration step k, and θ is the relaxation parameter. Note that Equations (34)-(35) depend implicitly on time because of J,J p at the right-hand sides. Hence, the equations must be solved at every time step. The functions J,J p are taken from the n-th time step. The numerical solution of Equations (36b)-(39b) is based on a similar procedure. Writing (36a)-(39a) in the standard form yields: where the coefficients A i , B i , C i , D i and F i depend on the other unknowns and their partial derivatives. We apply the standard second-order finite-difference approximation to the second derivatives and the upwind first-order approximation to the first derivatives. For instance, (36b) is approximated as follows where J n i,j = J i,j (x i , y j , nτ ),τ is the time step, Δx = 1 N , Δy For i = 0, N and j = 0, M (36c) is replaced by the boundary conditions. The resulting tri-diagonal matrix is inverted by the Thomas algorithm. The first-order upwind differences in (36c) result in the first-order approximation; however, the corresponding matrix is diagonally dominant. Consequently, the Thomas algorithm is stable (Bortoli et al., 2015). The system (36b)-(39b) is solved by the following iterative algorithm.
where U = J, J p , ψ, ψ p , T, T p , and ε is the required accuracy. 5. If ≥ ε return to step 1.

Results and Discussion
Testing hydrodynamic models of the TMDD is hampered by a lack of available data. Up to now, the measurements of the MN flow inside the vessels during the TMDD are not available. However, several experimental studies performed to measure the temperature of the nanofluid have been published. Therefore, the model is validated by 1) comparing it with the preceding models and 2) comparing it with the published results of the measurements of the temperature of the nanofluid.

Validation of the Model
We compare the proposed two-phase, two-way model (TPM) with the single-phase model (SPM 1 ) of Tzirtzilakis (2008), the single-phase model of Bianco et al. (2009) (SPM 2 ) and the double-phase, one-way coupled model (OW) of Boutopoulos et al. (2020) Consider Case 1: 25 ≤ Re ≤ 500, β 2 ≈ 4, St ≈ 3.56 × 10 −4 , C v = 0.001, and d = 250. Since C v is small, the single-phase and the double-phase models show similar results. Figure 3 displays the streamlines and the contour lines of the magnitude of the velocity V for the steady-state solution obtained by the TPM. Figure 4 shows the profile of |V| for x = 5 (position of the magnet). The difference between the solutions by the TPM and the reference models does not exceed 10 −4 .
The above tests are considered as a partial validation of the model. However, for large C v = 0.03 the proposed model shows different results. As an example, consider Case 2: β 2 ≈ 2, St = 2.8 × 10 −2 , 25 ≤ Re ≤ 500, d = 20000 and C v = 0.03. Figure 5a-c shows the steady-state solution for SPM 1 , OW, and TPM. Figure 6 shows the profile of |V| for x = 5. Clearly, in the case of a high concentration of MNs, the results obtained by the proposed model are different since the MNs have a reverse impact on the dynamics of the blood. For instance, the maximum difference between the solution by SPM 1 and the proposed model for x = 5 is about 0.9 cm/s. Wen and Ding (2004) experimentally evaluated the average Nusselt number using a straight copper tube of 970 mm length, 4.5 mm inner diameter, and 6.4 mm  outer diameter in the laminar regime. The setup consisted of four units: the flow loop, the measuring unit, the heating unit, and the cooling unit. The flow loop further included four sections: the pump with a built-in flow meter, the test section, the reservoir, and the collection tank. The test section consisted of a straight copper tube heated by a silicon rubber flexible heater linked to a DC power supply. The nanofluid flowed through the copper tube. Thermocouples were used to measure the temperature at the inlet and the outlet flow. The Nusselt number is defined by: where h C is the diameter of the tube, andh is the local heat transfer parameter given bȳ T m (ȳ) is the mean temperature of the fluid,T w (ȳ) is the temperature of the wall, and q is the convective heat    (3), whereρ = 996 kg/m 3 ,c = 4180 J/kgK. The average least square error for SM 1 is 4.8%, for OW is 3.3%, and for SM 2 is 3.4%. The TPM shows the least square error of about 5%. Therefore, the accuracy of the proposed TPM is comparable with the preceeding models when the MF is not applied.
Our second test apllies to the experimental set-up by (2020). The experiment was carried out on the flow of water and Fe 3 O 4 through a circular tube under the constant MF. The set-up is similar to Wen and Ding (2004). A straight circular copper tube (2700 mm) with 7.7 mm outer and 0.7 mm inner diameter had a thermal conductivity of 385 W/mK. Six wires were installed to produce an MF of about 0.13 T. C v ∈ [0.005, 0.01], d = 250. The average Nu obtained by the reference models and the proposed TPM are compared with the experimental results by (2020) in Table 3.
The average least square error for SM 1 is 12%, for OW is 8.3%, and for TMP is 5.7%. Therefore, the TPM shows a better accuracy in the presence of the MF.

Numerical Experiments and Discussion
Following Nacev et al. (2011) a magnet is located at 5, 10, and 15 cm below the wall. The MF, characterized byB = 0.5 T andB = 1.0 T , is applied to the blood flow mixed with Fe 3 O 4. The MNs are characterized byM ≈ 100 kA/m and diameters of 250, 800, and 20,000 nm. Note that the diameter of a human blood vessel ranges from 7 μm in the capillaries to 3 cm in the vena cava. The numerical simulations are performed forh = 2.0 cm andL = 30 cm.These are a typical diameter and length of the large blood vessel considered by the reference BFMs. The MNs and the blood enter the vessel at the left boundary. The left boundary condition for u is a parabola such that max(ū) =ū r1 at the centerline andū = 0at the walls. The density of the MNs and the blood areρ = 1050 kg/m 3 andρ p = 5200 kg/m 3 , respectively. The dynamic viscosity of the blood isμ = 3.2 × 10 −3 kg/(m s). The Reynolds numbers Re p ≈ 9.8 × 10 −4 and 3.1 × 10 −3 , for d = 250 and 800, respectively. The drag factor f = 1. When d = 20000, Re p ≈1.0 and f = 1.1. The thermal conductivityk = 2.2 × 10 −3 J/(m sK). The heat capacity of the MNs and the bloodc = 14.65 J/kgK andc p = 670 J/kgK, respectively. Hence, the Prandtl number Pr =¯cμ k ≈ 21. The parameters to control the transfer of the momentum and energy between the MNs and the blood are given by (Gireesha et al., 2017).

Impact of the Reynolds Number
The following examples analyze the flow patterns generated by the MF and the DF. The vorticity of the blood flow is an important factor that affects the patient during and after the TMDD. It has long been established that the vorticity of the blood flow is the cause of the cardiac dysfunction. Moreover, the blood vortices indicate physiological changes in the surrounding system, and can provide early indications of the long-term cardiac outcome (Pedrizzetti et al., 2014). The effects of the MF on microcirculation and microvasculature are not clear or widely explored. However, many studies have indicated that MFs could trigger either vasodilation or vasoconstriction (Contijoch et al., 2020;McKay et al., 2007). In the context of this work, we also refer to hemolysis which is the mechanical damage of the red blood cells from an excessively high stress induced by high gradients in the blood flow (Bletsos et al., 2021). Many studies have been performed to experimentally detect vascular inconsistencies. The medical research shows that the preferred vascular blood flow mode is laminar (Pedrizzetti et al., 2014). However, during the TMDD, the blood flow is affected by the MF and the exchange of the momentum and the energy between the MNs and the blood. Given that cancer chemotherapy requires from 3 to 6 months, there is a serious concern regarding the long-term effects of the TMDD on the vascular system. Hence, the numerical analysis of the coupled blood-MN flow may become indispensable, to study the patterns of the coupled blood-MN flow.
Consider three cases of the two-phase blood-MN flow characterized by C v ≈ 0.001,M ≈ 10000,λ a = 200 μm, and C c ≈ 4 (Kenjereš & Tjin, 2017). The numerical solution for Re ≈ 25, 50 and 100 is shown in Figures 7 and 8. The magnet is located at (a,b) = (10, 10). The magnetic inductionB = 0.5. Figure 7 shows the streamlines of the MN-blood flow and the contour lines of the vorticity function.
The streamlines are shown on the top of the density plot.
Case 1:ū r = 3.8 × 10 −2 . In this case Re ≈ 25, Re p ≈ 9.8 × 10 −4 , St ≈ 3.56 × 10 −4 , M nF ≈3281. This set-up complies with the data published by Nacev et al. (2011). The circulation starts at t C ≈ 13 minand reaches the steady state at t S ≈ 1.5 h. The MF generates a lower-wall vortex. The blood near the upper wall responds, creating a backflow. The velocity in the lower part of the vortex is minimal which makes the TMDD possible (Figure 7 (a and b)). The velocity at the upper part of the vortex reaches 30 cm/s. Case 2:ū r = 7.6 × 10 −2 corresponds to, Re ≈ 50, Re p ≈ 4.1 × 10 −3 , St ≈ 7.12 × 10 −4 , and M nF ≈ 820.In this case, the size of the region where the velocity is low or close to zero (the TMDD regions) decreases (Figure 7b). Case 3: Finallyū r = 15 × 10 −2 , where Re ≈ 100, Re p ≈ 1.0 × 10 −2 , St ≈ 1.42 × 10 −3 and M nF ≈ 205. There is a further decrease in the size of the TMDD regions ( Figure 7c).
The graphs in Figure 7 show that the DF works against the MF, decreasing the efficiency of TMDD. Clearly, to improve the dynamics of the blood flow, the magnet must be placed closer to the TMDD zone. Alternatively, the intensity of the MF must be increased. Furthermore, the contour plots of the vorticity function in Figure 7a1-c1 show that the vorticity increases as the Re increases. However, the vorticity plot for the smallest Re ≈ 25 (Figure 7a1) has the most complex geometry.
The streamlines of the flow and the countour lines of the vorticity function for the steady-state flow are shown in Figure 8.
Case 1: The MN flow creates four vortices with the antiparallel velocity vectors near the position of the magnet (x = 5). The TMDD zones, characterized by the low magnitude of the MN flow, are near x = 5 and x = 6.5. The maximum velocity of the MNs is about 28 cm/s. Case 2: The topological structure of the MN flow is approximately the same as in Case 1 (Figure 8b). However, the TMDD zone becomes smaller.
Case 3: Figure 8c shows generation of the new vortices. A double-vortex is on the left side of the vessel in Figure 8c2. The area of the TMDD zones is small. However, the area of the regions with the large horizonal velocity is also getting smaller.
Although the Stokes number for the numerical experiments is relatively small, i. e., the maximum St ≈ 1.42 × 10 −3 , the solid phase (MNs) does not follow the blood flow because of the impact of the MF.

Impact of the Magnetic Field
In order to evaluate the impact of the MF, consider Re p ≈ 9.8 × 10 −4 and Re = 250. The vertical position of the magnet is 15, 10, and 5 cm below the wall. The topology of the blood and MN flows is visualized in Figure 9-12. We consider the TMDD to be efficient if the blood velocity is low and the MN velocity toward the magnet is high.
When B = 0.5 , the magnet located at 15 cm does not have a significant impact on the MNs. As a result, the flow becomes laminar with the exeption of the region affected by the MF (see Figure 9a). Moving the magnet closer to the vessel expands the region of blood circulation. The TMDD zone where the velocity is close to zero increases. The topological pattern of the flow remains approximately the same for b = 10 ( Figure 9b). Moving the wire closer (b = 5) creates 3 vortices with the maximum velocity of about 15 cm/s. The TMDD zones expand significantly. Figure 10 displays the streamlines of the MN flow for B = 0.5 . Clearly, the MF attracts the MNs to the region of the required TMDD. However, the topology of the flow is complex. It does not follow the streamlines of the blood flow when the impact of the MF is strong (Figure 10).
When B = 1.0 the source located at 10 and 15 cm below the vessel (Figure 11 (a and b)) generates a flow topologically similar to that shown Figure 9 (a and b). However, the TMDD region becomes considerably larger. When B = 1.0 b = 5 (Figure 11c), the impact of the MF is significant. However, the topology of the flow is complex, and is characterized by an increased vorticity. Figure 12 for B = 1.0 shows how the MF works against the DF. The topology of the MN flow in Figure 12 is characterized by four large vortices. The magnetic numbers increase from M nF = 3281 and M nF p ≈ 6677 when B = 0.5to M nF ≈ 6562 and M nF p ≈ 13354 when B = 1.0. Figure 12c shows that an MF with a high intensity positioned close to the blood vessel has a strong impact on the vorticity of the blood flow. The evaluation of its impact on the patient and the negative effects of the vorticity versus the positive results of the TMDD is an open problem. This requires further experimental research outside the scope of this paper. However, the proposed model can be used as a second opinion to evaluate possible side effects.

Impact of the Size of the MNs on the Blood Flow
Let (a,b) = (10, 10), C v = 0.04, B = 0.5 and Re ≈ 250. Three cases are considered.
The magnetic numbers M nF = 3281, M nF p ≈ 6677. The streamlines of the blood flow for Cases 1,2, and 3 are shown in Figure 13. Clearly, the size of the MNs has a significant impact on the blood flow. The circulation is characterized by vortices and backflows near the magnet. The TMDD zone does not change significantly. However, the region having the high velocity of about 13 cm/s (yellow highlight) grows as the size of the MNs increases. Figure 14b (d = 800) shows three well-defined vortices, whereas Figure 14c (d = 20000) shows two vortices moving in the opposite (clockwise and anticlockwise) directions. The MN flow is toward the TMDD. The maximum velocity is about 12 cm/s nearby the magnet. The blood/MN flow in Figure 13c and Figure 14c is suitable for the TMDD.    Further, Re p ≈ 0.01, Re = 250, u r = u r1 = 3.8, Pr = 25, Ec = 4.2 × 10 −5 , M nF = 3281 , and M nF p = 6677. The thermal conductivity of the blood increases with the increase of the size of the MNs. Recall that the exchange of the energy between the blood and the MNs is proportional to parameter β 2 (see Equations (38a)-(39a)). This parameter is evaluated by (Gireesha et al., 2017) as follows β 2 = 6 ρ pcpk d 2 (2.0 + 0.6Re p 1/2 Pr 1/3 )h u r .
The impact of the size of the MNs on the temperature is illustrated in Figure 15 and 16. The temperature increases in the vicinity of a strong MF owing to the magnetocaloric effect and the energy exchange defined by Equation (44). The temperature increase inside the lower vortex is from 37to 39.1 • C. Note that the case β 2 = 0 in Figure 15a (magnetocaloric effect) complies with the results obtained by Tzirtzilakis (2008). For d = 250, the temperature starts increasing in the MF region at t ≈ 13 minand reaches approximately39.1 • Cat the steady state t S ≈ 2 h. When d = 800, the temperature reaches 39 • .1 for 1.5 h whereas d = 20000 requires only 0.5 h. The maximum temperature of the MNs reaches approximately 38.8 • Cirrespective of d (Figure 16). Such a temperature increase may have a negative impact on long-term treatments such as chemotherapy. As an example, we refer to the effects of the MRI with the MF of 4 T (Schenck et al., 1992;Yamamoto et al., 2004). The patients report sensations of nausea, vertigo, metallic taste, or sleepiness after treatments.

Conclusions
The proposed two-phase, two-way coupled model makes it possible to simulate dynamics of MNs in the blood flow under the impact of strong MFs. The most important feature of the model is the possibility to evaluate the inverse impact of the MNs on the blood flow.
The numerical experiments show the possibility to evaluate the impact of the Reynolds number on the topology of the blood flow and the MN flow. For Re > 100 and the MF characterized by B ∈ [0.5, 1] T, the DF works against the MF. The TMDD region (low blood velocity) becomes small and the magnitude of the MN flow in this region decreases. The model is capable of simulating the impact of the MF and the position of the magnet. Positioning the magnet close to the desired TMDD region and applying a strong MF (B = 1 T) generates large TMDD regions. However, the topological structure of the blood flow is characterized by an increased vorticity. In particular, the closest position of the magnet (b = 5 cm) generates a complicated flow with several vortices extending across the entire vessel. Our conjecture (supported by the referenced experimental research) is that the increased vorticity may harm the patient in the case of a long-term treatment. The model shows that the size of the MNs has an impact on the dynamics of the flow. Large nanoparticles create a large TMDD zone where the MNs move faster to the target. However, it is well known that the large MNs damage the blood vessels and cause blood leakage. This effect is the subject of further research. The proposed two-way coupled energy equations make it possible to analyze the impact of the size of the MNs on the temperature of the blood flow. The model shows that under certain conditions, the large MNs induce a temperature increase of about 2°C over a relatively large region. Such an increase combined with the increased vorticity may lead to a negative impact on the patient.
The model has a number of limitations. The magnet is represented as a pointwise wire. However, the MF depends on the shape and the orientation of the magnet. This effect has not been considered. Further, the Lorentz force is not included in the equations. The dependance of the heat transfer on the MF has not been considered either. Furthermore, the model does not include a number of important effects such as the Brownian motion of the MNs, the acoustic radiation force, the lift force (Saffman and Magnus effects), and the thermophoretic force. Finally, the shape of the MNs (an important factor) has not been included in the proposed two-phase equations. Therefore, future research requires further validation of the model using new experimental results and modifying the model to include the above mentioned effects.