Numerical Solution of Biomagnetic Power-Law Fluid Flow and Heat Transfer in a Channel

: The effect of non-Newtonian biomagnetic power-law ﬂuid in a channel undergoing external localised magnetic ﬁelds is investigated. The governing equations are derived by considering both effects of Ferrohydrodynamics (FHD) and Magnetohydrodynamics (MHD). These governing equations are difﬁcult to solve due to the inclusion of source term from magnetic equation and the nonlinearity of the power-law model. Numerical scheme of Constrained Interpolation Proﬁle (CIP) is developed to solve the governing equations numerically. Extensive results carried out show that this method is efﬁcient on studying the biomagnetic and non-Newtonian power-law ﬂow. New results show that the inclusion of power-law model affects the vortex formation, skin friction and heat transfer parameter signiﬁcantly. Regardless of the power-law index, the vortex formation length increases when Magnetic number increases. The effect of this vortex however decreases with the inclusion of power-law where in the shear thinning case, the arising vortex is more pronounced than in the shear thickening case. Furthermore, increasing of power-law index from shear thinning to shear thickening, decreases the wall shear stress and heat transfer parameters. However for high Magnetic number, the wall shear stress and heat transfer parameters increase especially near the location of the magnetic source. The results can be used as a guide on assessing the potential effects of radiofrequency ﬁelds (RF) from electromagnetic ﬁelds (EMF) exposure on blood vessel.


Introduction
Biological fluid such as blood includes protein, cell membrane and, haemoglobin molecules that contribute to its paramagnetic behaviour, while the content of hydrogen, oxygen, nitrogen and carbon atoms in vessel tissues contribute to its diamagnetic property [1,2]. Biofluid that has magnetic characteristics and able to interact with an applied magnetic field is called biomagnetic [3]. When an external magnetic field is imposed on the biomagnetic fluid, there are two body forces acting on the fluid, magnetisation force and Lorentz force. These two forces correspond to the principle of Ferrohydrodinamics (FHD) and Magnetohydrodynamics (MHD) respectively. The interplay between its hydrodynamics and magnetic properties when external magnetic field gradient is applied has been the topic of great interest [4]. Motivation in these studies are driven by the biomedic application such as in magnetic drug targeting in tumor treatment, hyperthermia, reduction of bleeding during surgeries, and magnetic separation of cells [5][6][7][8][9][10][11]. However, at the same time, there are some drawbacks of However, most of the works mentioned are concerned with biomagnetic of Newtonian fluid. When blood flows through larger arteries, it behaves according to Newtonian model, but in small blood vessels such as narrow arteries, it behaves according to non-Newtonian flow model at low shear rates [50][51][52][53][54]. Unlike Newtonian fluid, non-Newtonian are a broad class of fluid where the relation between the shear rate and shear stress is nonlinear. Hence, no single constitutive equation can be used to single handedly predict all the flow of non-Newtonian models in various condition. There are several models characterizing flow of blood such as Walburn-Schneck, Casson, Carreau-Yasuda, Quemada, and power-law, where each has its own importance in theoretical and experimental of flow phenomena [55][56][57][58][59]. The power-law model is the simplest form and commonly used for modelling shear-thinning or shear-thickening behaviour of non-Newtonian fluids on certain range.
On the study of purely hydrodynamics power-law flow, Bell and Surana [60] investigated the flow of non-Newtonian power-law in lid driven cavity using p-version least squares finite element formulation and provided the benchmark results for shear thinning and shear thickening flow in lid driven cavity. Their results on power-law lid driven cavity were validated by Neofytou [61], who also carried out simulations for Bingham plastics, Casson and Quemada fluid flow. There are many works that have been reported on magnetic field effect on power-law fluid flow. Kefayati [44] studied the effect of MHD magnetic field on power-law in cavity with two-square concentric duct. Kefayati [45] studied the effect of magnetic field on mixed convection on shear thinning fluids in a cavity with sinusoidal boundary conditions. Adıgüzel and Atalık [46] investigated the effect of MHD magnetic field on power-law ferrofluid flow past a circular cylinder. Jahanbakshi [47] reported the MHD magnetic field effects on natural convection using power-law fluid in an L-shaped enclosure. Ikbal et al. [48] and Sankar et al. [49] investigated the effect of MHD magnetic field on power-law in a stenosed artery. All of these works conclude that the vortex formation decreases under the action of MHD magnetic field. A number of investigations reported works on MHD effect on non-Newtonian flow other than power-law model. Misra et al. [62] studied the works of MHD viscoelastic fluid in a channel with stretching walls and reported that increasing the strength of the magnetic field will have the fluid velocity profile decreases but the temperature profile increases. Liu [63] performed an analysis of electrically conducting fluid of second grade in a porous medium over a stretching sheet subject to a transverse magnetic field and found that when the Eckert number is large, the heat flow transferred from the fluid to the wall rather than from the wall to the fluid when Eckert number is small.
The works that combine both biomagnetic and power-law fluid flow has not been investigated. The reason is that the mathematical model of Navier Stokes with FHD and MHD is difficult to solve due to the magnetisation term FHD introduced in the governing equation. This source term can be very large and will result in instability and divergence of the numerical method without the use of a proper technique. Additionally, in contrast to the Newtonian case, the non-Newtonian fluid flow is generally more complicated to be solved due to its complexity contributed from the constitutive equation to the Navier-Stokes equations. This is due to the fact that the effective viscosity diverges to infinite shear rates and zero shear rates for shear thinning and shear thickening fluid. When solving non-Newtonian flow numerically, it may result in significantly time step constraint when explicit schemes are used, and disrupt the solution technique usually used to solve implicitly the viscous terms in the momentum equation [64]. Several methods have been used to solve the problem analogous in this work. Stream function-vorticity finite difference method (FDM) with stretched grid has been used in [33,65]. FDM central difference has also been used in [48]. In [44,45], finite difference lattice Boltzmann method (FDLBM) is used, while in [35,36] finite element method (FEM) has been applied. In [42,43], dual reciprocity boundary element method (DRBEM) has been used. Lattice Boltzmann method (LBM) has been used in [66,67] while control volume finite element method (CVFEM) is applied to solve the flow problem in [68][69][70]. More detailed review of numerical method on solving FHD and MHD problems are given by Rashidi et al. [71].
An explicit CIP (Constrained Interpolation Profile) with FDM scheme was developed in this work. This scheme is a very stable finite difference technique that can solve generalised hyperbolic equations with third order accuracy in space. Unlike the conventional schemes such as Lax Wendroff and Upwind, this scheme approximated the values of the advection term with its spatial derivative by cubic polynomials interpolation and attains high-accuracy by being less diffusive [72]. The CIP method was introduced by T.Yabe et al. [73,74] and has since been used to solve many various problems [75][76][77]. This method can be combined with other methods [78] as well since CIP is efficient in solving the hyperbolic equation.
As far as we know, no previous research has investigated the works on biomagnetic and non-Newtonian fluid flow using CIP scheme. Therefore the aim of this work is to study the extent of using non-Newtonian power-law on biomagnetic flow in terms of heat transfer, skin friction and flow characteristics where the governing equations are solved numerically using CIP scheme. For non-Newtonian power-law fluid, the viscosityμ a is assumed to be a shear dependent fluid given as where m is the fluid consistency, I 2 is the second invariant of the rate of strain tensor and n is the power-law index. For two-dimensional flow, I 2 is defined as whereū is the velocity alongx-coordinate, andv is the velocity alongȳ-coordinate. For values of n < 1, the viscosity is shear thinning, while for n > 1, the viscosity is shear thickening and n = 1 refers to the Newtonian viscosity model whereμ a = m. Biological fluids have shear thinning properties as blood has a power-law index at n = 0.6 and synovial fluid has a power-law index at n = 0.2 [46]. In addition, this study may give valuable information on the thermal and hydrodynamics effects on blood associated with RF exposure.

Problem Formulation
The continuity, momentum and energy equations for biomagnetic fluid flow are given as [24,71] where v is the dimensional velocity vector,t is the dimensional time, ρ is the fluid density,p is the dimensional pressure, τ is the stress tensor, J is the current density,B is the magnetic field,μ 0 is magnetic permeability,H is the magnetic field intensity,T is the dimensional temperature, c p is the specific heat at constant pressure, κ is the thermal conductivity, andσ is the electrical conductivity. The fluid is considered to be magnetic non-Newtonian and the flow is two-dimensional, laminar, incompressible and under the influence of an applied magnetic field. The body force terms incorporated are from the principles of Magnetohydrodynamics and Ferrohydrodynamics. The flow is considered to be taking place in between two walls. It is assumed that the flow is dominated by the viscous force and the convection is due to temperature gradient, hence the effect of gravity is insignificant. The flow at the entrance is assumed to be fully developed and the walls are kept at a constant temperaturē T w , while the fluid is at temperatureT f such thatT w <T f . The configuration of the flow is shown in Figure 1. With these assumptions, the governing equations in Cartesian form are given as ∂ū ∂x ρ ∂v ∂t In the above equations, v = (ū,v) is the velocity,p is the pressure,T is the temperature,ū(ȳ) and T(ȳ) is the parabolic velocity profile and temperature profile respectively for fully developed flow,σ is the electrical conductivity,μ 0 is the magnetic permeability vacuum, c p is the specific heat at constant pressure, κ is the thermal conductivity,H = (H x ,H y ) is the magnetic field strength,M represents the magnetisation, andB is the magnetic induction whereB =μ 0H ⇒ (B x ,B y ) = µ 0 (H x ,H y ). The bar symbol on u, v, T, t, x and y denote the dimensional form of these quantities.
The components of magnetic field intensityH x andH y along thex andȳ coordinates H = (H x ,H y ) are given by [40] and the magnitudeH is given bȳ where γ is the magnetic field strength at the coordinate point of (a, b) where the magnetic source is located. In the magnetic force terms, magnetisation and magnetic field intensityH derived experimentally and is given asM whereK is a constant andT c is the Curie temperature.
The termsμ 0M ∂H ∂x andμ 0M ∂H ∂ȳ in Equations (7) and (8) respectively represent the component of magnetic force per unit volume and dependent on the existence of magnetic gradient. The term in Equation (9) refers to the thermal power unit per unit volume due to the magnetocaloric effect. These three terms arise due to the principles of FHD. The termsσB 2 yū + σB xByv andσB 2 xv +σB xByū in Equations (7) and (8) respectively represent the Lorentz force per unit volume and included due to the electrical conductivity of the fluid. The termσ(ūB y −vB x ) 2 in Equation (9) represents the Joule heating. These three terms arise due to the principles of MHD. Therefore the principle of FHD and MHD are combined in the mathematical model with the application of non-Newtonian fluid.
The non-Newtonian fluid considered here can be represented by the Oswald-de Waele power-law model. The relationship between velocity field and shear stress may be expressed in terms of the effective viscosity asμ where m is the fluid consistency and n is the power-law index. The boundary conditions of the problem are

Non-Dimensionalisation and Transformation to Stream Function-Vorticity
Equations (7)- (9), with the boundary conditions of Equation (14) are transformed into dimensionless form using the following dimensionless variables whereH 0 =H(a, 0) and u r is the maximum velocity at the inflow.
To solve the governing equation, stream function-vorticity approach is used and the velocity v = (u, v) is transformed into dimensionless stream function ψ(x, y) and vorticity w(x, y) by the relation Using these transformation, the pressure p is eliminated, and Equation (6) is automatically satisfied and hence only the nondimensional variables ψ, w and T need to be solved numerically. The transformation results to the following system of equations The effective viscosity is nondimensionalised to be The magnetic field intensity is also nondimensionalised by using Equation (15) and is given by The boundary configuration for the flow problem is also nondimensionalised using Equation (15) where the length of the horizontal wall is considered to be 10 times the length of the vertical inlet and outlet. The boundary conditions are non-dimensionalised and transformed according to Equations (16) and (17), and given as

Constrained Interpolated Profile (CIP) Based Finite Difference Method
To solve Equations (18)-(23) numerically, the Constrained Interpolated Profile (CIP) method is implemented. In this method, a cubic polynomial is used to interpolate the spatial profile of the value f and its spatial derivative ∇ f . In this case, f is f = (w, T). Once the values of f and ∇ f are calculated, the same values for the next time step are simply obtained by shifting the cubic polynomial.
To implement the CIP algorithm, the governing equations of Equations (19) and (20) are split into two categories [79] which are non-advection phase and advection phase.
Advection Phase: Non-advection Phase: The computation of advection phase, Equations (24) and (25) is conducted by the CIP method while the non-advection phase, labeled with superscript (*) is solved numerically using finite difference method. In order to adopt this method, the Poisson equation of Equation (18) need to be solved first, then followed by the non-advection phase and finally the advection phase.
Equation (18) is solved by using Successive Overrelaxation (SOR). The iterative procedure is carried out until the residue between two consecutive iteration converges to desired tolerance. The tolerance for SOR is suggested to be smaller than the tolerance for the steady state solution in which |ψ new − ψ old | < 10 −7 , and the relaxation factor used is 1.7.
For the non-advection phase, the explicit scheme is used. The central finite difference discretisation is implemented on first order and second order spatial derivatives terms. The magnetic terms ∂H ∂x ,

∂H ∂y
, H x and H y are computed analytically from Equation (22). It should be noted that the method has also been tested to be solved using Crank-Nicolson scheme and no significant difference is observed. However using the implicit scheme will require solving the two matrices equation in every time step for each Equations (19) and (20), and inadvertently it is time consuming due to exhaustive computer power required.
Once the non-advection phase are computed, the spatial derivative are calculated for both w * and T * . Hence by denoting f = (w * , T * ), the spatial derivative The terms f k , ∂ x f k i,j , and ∂ y f k i,j refer to the values of previous time step. The spatial derivatives of u and v in Equations (28) and (29) are approximated by using central order finite difference.
Once f * = (w * , T * ) and their respective spatial derivative are obtained, they will be used in CIP algorithm. Hence f k+1 = (w k+1 , T k+1 ) are advanced to the new time level after the interpolation of CIP via shifting such that f k+1 = f (x i − u∆t, y j − v∆t). Similarly in non-advection phase, the spatial gradient for the new time iteration is also required. The CIP scheme to solve Equation (24) and (25) is given as ∂ f k+1 i,j ∂y = 3C 03 ∆y 2 + 2C 02 ∆y + ∂ f * ∂y + C 21 ∆x 2 + 2C 12 ∆x∆y + C 11 ∆x.
The term C is the coefficient for the interpolation function Equation (30), and for ease of reference, the subscript on C denotes the power of its multiplication with ∆x and ∆y respectively. Since , The above equations are only for u < 0 and v < 0. For the case of u ≥ 0 , the pointer i + 1 is to be changed to i − 1 and ∆x is changed to −∆x. Similarly, when v ≥ 0, j − 1 and −∆y are used instead of j + 1 and ∆y. The calculated f k+1 = (w k+1 , T k+1 ) and their respective spatial gradients are then updated to be f k and will be used for the new time level. The whole procedure is repeated until a steady state flow is achieved in which the tolerance is set to be where φ k is the estimation of φ k = (ψ k , w k , T k ) at iteration k.

Boundary Conditions
For the numerical solutions for the system of Equations (18)- (20), boundary conditions for flow variables w, ψ and T are discretised according to Equation (23). The boundary conditions for w, ψ, and T on the entrance and the exit are easily implemented. On the walls, the boundary conditions for ψ and T are also implemented according to Equation (23). However, the boundary condition for w on the walls are not straightforward.
The boundary conditions for the velocities on the lower and upper walls are set as where J refers to the index of the lower and upper walls. By using boundary conditions of Equation (34) and central differences of Equation (17), we obtained Equation (18) is then discsretised using second order central difference, and substituting the relations from Equations (35) and (36), the boundary conditions for the vorticity on the lower and upper walls are then proved to be [33] where fictitious points are to be used accordingly. The discretisation of boundary conditions used for w, ψ, and T during the advection phase and non-advection phase are similar. The procedure for CIP-Finite Difference method is summarised as follows: • Set the initial and boundary conditions according to Equation (23 (26) and (27).

•
Calculate the spatial gradient of the non-advection phase from Equations (28) and (29). • Advance the time iteration by solving numerically the advection phase and their spatial gradient of Equations (24) and (25) by using the CIP algorithm from Equations (30)-(32).

•
Repeat the above steps until the steady state is achieved by using the convergence criterion in Equation (33).

Dimensionless Parameters
To implement the numerical method, it is necessary to assign the values of parameters in the governing equation. Other than the power-law index n, the nondimensional parameters appear in the governing equations are The special definition of parameters Pr and Re above can be converted to conventional Newtonian fluid by taking the power-law index n = 1. The parameter Mn F is the Magnetic number that arises due to FHD, whereas the parameter N is the Stuart number, defined as the ratio of magnetic to inertia effects.
As in the works of [23,33,65], the magnetic fluid is considered to be blood as blood can be paramagnetic, diamagnetic and non-Newtonian fluid. The parameters are calculated by considering that ρ = 1050 kg m −3 , m = µ a = 3.2 × 10 −3 kg m −1 s −1 , with maximum velocity u r = 1.524 × 10 −2 ms −1 . The distance of the plates are L = 2.0 × 10 −2 m. Hence, the Reynolds number is at Re = 100. The location of the magnetic wire for applying the localised magnetic field is set to be (a, b) = (2.5, −0.05), where the maximum intensity occurred on H 0 = H(2.5, 0). The temperatures of the plate is assumed to beT w = 37 • C while the temperature of the fluid isT f = 41 • C, in which making the temperature number to be 77.5. The Prandtl number is considered constant at Pr = 25 by taking C p = 3.9 × 10 3 J kg −1 and κ = 0.5 J m −1 s −1 . From these values, the Eckert number is set to be Ec = 1.49 × 10 −8 .
The important dimensionless parameter for this work is the Magnetic number Mn F and it can be written as whereB 0 andM 0 = 40 A m −1 are the magnetic field induction and the magnetisation at the point (a, 0) respectively. The Stuart number N is calculated subsequently using the parameters assigned as in Mn F , hence Equations (38) and (39) show the relationship of Mn F and N with Re, where increasing the magnitude of Magnetic number is the same as lowering the Reynolds number. Hence, the Reynolds number Re is fixed to be Re = 100 throughout the numerical experiment, as varying the value of Mn F and N are sufficient. The parameter for Mn F and N assigned are dependent on the reference magnetic field induction ofB 0 in Tesla (T) unit. The values ofB 0 are chosen to be 1 T, 4 T and 8 T and the corresponding values for Mn F and N are given in Table 1. This variation of Tesla number was implemented to study the effect of applied magnetic field on power-law flows, in the case of shear thinning (n < 1), Newtonian (n = 1) and shear thickening (n > 1). Another important parameter in this study is the power-law index, that determine whether the fluid is Newtonian or non-Newtonian. The power-law index n is varied from n = 0.6 for shear thinning fluid, to n = 1 and n = 1.4 for shear thickening fluid. This scenario is apt for the magnetic fluid such as blood since its power-law index is varied in the human biological system. This can be also interpreted as using different class of non-Newtonian magnetic fluid for the flow in a channel.
Flow and heat characteristics that are being compared in this study are the local skin friction and the local rate of heat transfer coefficients. These quantities are defined as where C f is the local skin friction coefficient and Nu is the Nusselt number that describes the local rate of heat transfer coefficient. These coefficients are calculated both on the lower wall (y = 0) and upper wall (y = 1) with different Magnetic number Mn F and power-law index n number.

Grid Independence and Method Validation
Numerical solution of the system of Equations (18)- (20) has been developed using Matlab R . An extensive grid independence test is performed to check for the grid independent solution. The test is based on the u-velocity profile along the various x location of the channel, where the focus is emphasised near the location on the magnetic source. Five different size of grids configuration have been performed on the case flow problem n = 1, Mn F = 1312 and N = 0.056. The magnetic source is located at (a, b) = (2.5, −0.05). As shown in Figure 2, when the grid size is refined from 400 × 40 to 500 × 100, there is a significant change on velocity profile, particularly near the magnetic source and in the interval of the vortex size. The numerical simulation is sensitive on the vertical grid, while no significant change occurred when refining the horizontal grid. For the temperature profile, the local heat transfer on the lower wall of the channel has been used as a reference to test the grid independence. Only the lower wall is shown since most of the significant changes occurred near the magnetic source. The same grid independence test result is observed as seen in Figure 3, where the grids size have reached grid independence for temperature profile. These grids sizes have also been tested with smaller time step and no significant changes is accounted. It is observed that the grid size configuration of 400 × 80 is sufficient to capture the velocity profile change near the magnetic location.   The performance and the validation of present numerical method, are done separately on the benchmark cases for non-Newtonian flow in lid driven cavity, and for the flow in a channel undergoing the localised magnetic field. The results validation regarding the flow in a channel under the action of localised magnetic field is performed as in [33]. From Figure 4, it can be seen that the present numerical method developed is in good agreement with the compared studies. The second validation of the code is based on the non-Newtonian comparison. The computation involving non-Newtonian of power-law model in a lid driven cavity is validated with the works of Bell and Surana [60] and Neofytou [61]. In accordance with these works, simulations have been performed for Re = 100 with different values of the power-law index. For Newtonian flow case, the comparison is made with the benchmark works of Ghia et al. [80]. The grid size used is 150 × 150, with the dimension of the wall is 1 × 1 where the top wall is moving with velocity u = 1. Figure 5a shows the results for u-velocity along the vertical centreline and the Figure 5b v-velocity along the horizontal centreline. The non-Newtonian cases for n = 0.5 and n = 1.5 show very good agreement when compared with the results by Bell and Surana [60] and Neofytou [61]. For n = 1, the result is in very good agreement with the result by Ghia et al. [80].

Results and Discussion
The time step is chosen to be ∆t = 10 −4 for Newtonian and shear thickening cases while for shear thinning cases, the value of ∆t = 10 −5 is required for the numerical method to converge. For shear thinning fluid the lowest number of power-law index in this work that converged is n = 0.6 while the shear thickening cases were tested until n = 1.5, although only the results for n = 0.6, n = 1 and n = 1.4 are shown.
The flow is purely FHD by setting N = 0, and will be purely MHD if Mn F = 0. When both Mn F and N are zero, there is no magnetic source and the flow is purely hydrodynamics hence the flow in the channel has horizontal streamline corresponds to the parabolic profile for u-velocity. When the magnetic source is applied, a vortex arises near the area of the magnetic source. This phenomenon of how the flow is affected by the applied magnetic source with the variation of the power-law index is the main interest of this work. When the flow is tested on pure MHD, the effect of the magnetic field is not significant as the value of Stuart number N that is required to induced the vortex in the area of magnetic source has to be in the order of hundreds and the change is negligible even when N is increased into the order of thousands. Figure 6 shows the streamlines contours for variations of power-law index n = 0.6, 1.0 and n = 1.4, with variation of reference magnetic fieldB 0 = 1 T,B 0 = 4 T, andB 0 = 8 T. The notable effect of applying magnetic field on the flow is the formation of arising vortex near the location of magnetic field. As discussed in previous literature, this vortex extension increases with the increment of magnetic field from 1 T to 8 T. Interesting results are observed when the power-law index is increased from n = 0.6 to 1.4, the disturbance in the flow field is also affected. For the case of shear thinning fluid n = 0.6 it can be seen that the extension of the vortex is large and extended towards the downstream, compared to the case of when n = 1.0 and n = 1.4. Additionally when n = 1.4, the vortex formation is smaller although its size increases when the Magnetic number is increased. The measured range of sizes of the vortex for each cases summarised in Table 2.  Figure 7 shows the corresponding temperature contours for variations of power-law with variation of reference magnetic corresponds to the case presented in Figure 6. It is shown that in the area of magnetic source, the maximum temperature is displaced towards the upper plate. The effect is more pronounced when the Magnetic number is increased, regardless the power-law index. The region of the lower vortex is shown to be kept cooler and this region is extended towards the downstream when the Magnetic number is increased for all the cases of different power-law index.
From the previous discussion, it can be seen that the influence of increasing power-law index number n has the opposite effect to the increasing Magnetic numbers on the flow field. Table 3 shows the minimum velocity for different Magnetic numbers at various power-law index. The negative symbol indicates the direction of the velocity and the reverse flow of the vortex region near the magnetic source. From Magnetic number 1 T to 8 T the minimum u-velocity is increased for every n. Overall, the shear thickening flow has the highest velocity followed by Newtonian and shear thinning flow.    The variation of wall shear stress and heat transfer parameters of various Magnetic numbers for n = 0.6, n = 1.0 and n = 1.4 are shown in Figures 8-10 respectively. The wall shear stress are shown to have more influence on the lower plate, particularly at the location of the magnetic source. The variation of this parameter is qualitatively the same as magnetic field are varied from 1 T to 8 T, where the higher magnetic field results in greater variation of the wall shear stress irrespective the value of n. The variation of wall shear for the upper wall (Figures 8a, 9a and 10a) for all cases changes smoothly from increasing to decreasing near the area of the magnetic source, and its sign does not change. For the upper wall, the lower values of this parameter indicates the higher rates of skin friction. For the case of the lower wall (Figures 8c, 9c and 10c), the rapid changes of this parameter is largest for the case of shear thinning flow, followed by Newtonian and shear thickening flow as shown by the difference of the extrema values. For the lower wall, the negative values of the wall shear stress refers to inversion of the flow, and when wall shear stress reach its original value at x = 0, it refers to fully developed flow. This indicator is used to verify that all the flow cases managed to revert back to fully developed flow when it reaches downstream. The heat transfer parameters for the upper wall (Figures 8b, 9b and 10b) and lower wall (Figures 8d,  9d and 10d) similarly shows the same qualitative variation when magnetic fields are varied from 1 T to 8 T, for each case flow. On the lower wall, it is shown that the rate of heat transfer increases when Magnetic number increases, regardless the power-law index. For the case of n = 0.6, the overall heat transfer to lower wall is greater than the case of n = 1.0 and n = 1.4. The heat transfer coefficient is relatively high near the location of magnetic source, and another second region reaching the downstream due to the disturbance of the temperature field is extended far downstream. These second location of high heat transfer coefficient are observed at region 5.5 < x < 8.5, 4.5 < x < 6 and 4 < x < 5 for n = 0.6, n = 1.0 and n = 1.4 respectively. On the upper wall, the lower values of this coefficient refers to higher rates of heat transfer to the wall and it exhibits a similarity with the case of the lower wall.  Another quantitative comparison can be made from Figure 11, for different power-law index n and magnetic field at 8 T, where it shows that the wall shear stress and heat transfer parameters are highest for the case of shear thinning flow, followed by Newtonian and shear thickening. (c) C f Lower plate

Conclusions
In this paper, the effect of magnetic field on non-Newtonian biomagnetic power-law flow in a channel has been studied numerically by Constrained Interpolated Profile Finite Difference Method (CIP-FDM). This study has been conducted for the parametersB 0 = 1 T,B 0 = 4 T andB 0 = 8 T for power-law index n = 0.6 (shear thinning), n = 1 (Newtonian) and n = 1.4 (shear thickening). The investigation is performed based on the study of the streamline, temperature, vorticity contours, wall shear stress and heat transfer parameters. A good agreement validated with previous numerical investigation demonstrates that the CIP-FDM is an appropriate method for the power-law flow in a channel under an applied magnetic field. The enhancement of Magnetic number increases the length of vortex formation regardless the power-law index. The consideration of power-law index, causes the size of vortex formation to expand in shear thinning flow, and shrink in shear thickening flow. The amplification of Magnetic number increases both the wall shear stress and the heat transfer parameters dominantly near the location of magnetic source. The opposite effect is observed when the power-law is included as the wall shear stress and heat transfer parameters are decreased. The work reported here shows the severe drawbacks associated with RF radiation to blood vessel locally on the displacement of the high temperature and the blockage due to back flow. For future research, it is recommended to perform an analysis on fluid structure interaction by calculating the elastic deformation of the vessel.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: