Numerical and Analytical Study of Two-Layered Unsteady Blood Flow through Catheterized Artery

The pulsatile flow of blood in a catheterized blood vessel is analyzed. The flow of blood in vessel is modeled as the flow of two immiscible fluids. The fluid in the core region is characterized as a non-Newtonian viscoelastic fluid satisfying the constitutive equation of an Oldroyd-B fluid. The fluid in the peripheral region is treated as a Newtonian fluid. The catheter inside the vessel is modeled as a rigid tube of very small radius. The resulting differential system for velocity in each region is computed numerically by finite-difference scheme and analytically by Laplace transform. A comparison of numerical solution with Laplace transform solution is carried out. Various physical quantities of interest through the computed velocity are also analyzed.


Introduction
Medical catheters are used in medical centers to insert fluids or gases to patients or to deplete bodily fluids such as urine. Apart from that they are utilized for the diagnosis and treatment of various arterial diseases. Moreover, the monitoring of velocity and pressure gradient during any treatment/diagnostic procedure is also achieved by using catheters. Despite their advantages, there are certain consequences of clinical significance attached with their usage. It is observed that the insertion of catheter in an artery affects the flow field, disturb the pressure distribution and enhance the resistance to flow. Therefore it is necessary to study such increase in flow variables due to catheterization.
Very few attempts have been made in the past to study the blood flow through a catheterized artery. Back [1] estimated the resistance to the flow by assuming that the blood as a Newtonian behavior. He concluded that even a very small size of angioplasty guide-wire leads to sizable increase in flow resistance. In another study Back and Denton [2] provided the estimates of wall shear stress in coronary angioplasty and discussed its clinical implications. Mac-Donald [3] investigated the characteristics of blood flow in catheterized arteries using finite difference technique. Some possible effects of catheter on arterial wall are studied by Karahalios [4]. He observed that such effects become more significant when the annular gap between catheter and artery becomes narrow. The influence of catheterization on blood flow in a curved artery was discussed by Jayaraman and Tiwari [5].
The simultaneous effects of non-Newtonian nature of blood and catheter on flow characteristics are also reported by some researchers. For instance Dash et al. [6] utilized constitutive equation of Casson model to estimate the increase in flow resistance due to catheter in a narrow artery for both steady and pulsatile flow cases. The analysis of Dash et al. [6] was extended by Sankar and Hemalatha [7,8] for Herschel-Bulkley model.
On the other hand the unsteady blood flow through stenosed arteries is also attended by several investigators. For instance, Mekheimer and El Kot [9] analyzed the effects of Hall current on blood flow through an artery under stenotic conditions. Mustapha et al. [10] examined the characteristics of arterial constrictions and body acceleration in unsteady blood flow. In another attempt, Mustapha et al. [11] explored the blood flow in irregular multi-stenosed arteries. Here emphasis is given to the effect of applied magnetic field. A mathematical model for generalized Newtonian blood flow with irregular arterial stenosis is also developed and simulated by Mustapha et al. [12]. Abdullah et al. [13] addressed the influence of magnetic field on blood flow with irregular stenosis. A mathematical model for axisymmetric artery via cosine shaped stenosis is computed by Shahed et al. [14]. Apart from above mentioned studies, some investigators have also explored the effects of catheter on Newtonain and Non-Newtonian flows through stenosed arteries. Mentioned may be made to the studies carried out by Sarkar and Jayaraman [15], Dash et al. [16], Reddy et al. [17], Srivastava and Rastogi [18] and Muthu et al. [19].
It is noted from the available literature that blood flowing through catheterized arteries is taken as a single-phase Newtonian or non-Newtonian fluid. However, due to accumulation of red blood cells in the center of the larger arteries, a cell-depleted layer exists near the vessel wall. To account for such a situation, many researchers have treated the blood as a twophase fluid where the core region is modeled as non-Newtonian fluids while the cell-depleted layer is assumed as a Newtonian liquid. Additionally there are experimental evidences that blood exhibits viscoelastic properties under certain conditions [20,21]. Thurston [22] was the first to identify the viscoelastic nature of blood. He developed an extended Maxwell model, which is applicable to one-dimensional flow [23]. Yeleswarapu et al. [24] and Yeleswarapu [25] proposed an Oldroyd-B fluid model with three parameters for studying the blood flow. It is generally accepted that blood is slightly viscoelastic, and its effect was ignored in most of the computational fluid dynamics studies. At low shear rates, aggregate of RBCs behave like solid bodies and has abilities to store elastic energy. However, at high shear rates due to fluid-like behavior of RBCs, viscoelastic effects are also less prominent. Therefore, viscoelastic models are more adequate for blood flow at low shear stress and in oscillatory flow conditions [26].
Motivated by above facts, we are interested to investigate the flow characteristics of blood in a narrow catheterized artery by considering blood as a two-phase fluid. The aggregate of red blood cells in the core region is modeled by Oldroyd-B fluid while the periphery region is assumed to behave as a Newtonian fluid. The compliant nature of the artery is not incorporated in the present study. A priori estimate about how the compliant nature of the artery wall may affect the flow rate, stability etc may be difficult to obtain. However, the interested reader are referred to reference [27] for further details. The layout of the paper is as follows: Fundamental laws and geometry of pulsatile flow problem are described in section 2. The governing equations of the problem and appropriate boundary conditions are derived in section 3. Solution obtained via Laplace transform is presented in section 4. Numerical solution through finite difference method is obtained in section 5. Section 6 comprises results and discussion for various values of parameters of interest. Finally, some conclusions are drawn in section 7.

Flow Equations
The continuity and momentum equations governing the flow of an incompressible fluid are where u is the fluid velocity, ρ is the density, T is the Cauchy stress, b is the body acceleration and d/dt is the material time derivative given by The Cauchy stress tensor of an Oldroyd-B fluid is [26] T in which p is the pressure, I is the identity tensor and the extra stress tensor S satisfies the following expression where λ 1 , λ 2 are the relaxation and retardation times respectively, μ is the co-efficient of viscosity, A 1 is the first Rivilin-Ericksen tensor and D/Dt is the contravariant convective derivative. The expression for D/Dt and A 1 are defined by Where

Mathematical Formulation
Consider an axially symmetric, unsteady, uni-directional, incompressible and two-phase flow of blood through an artery of radius R in which a catheter is inserted. The radius of catheter is kR (0 < k < 1). The blood in core and peripheral regions is modeled by Oldroyd-B fluid and Newtonian fluid, respectively. We use cylindrical coordinates (r,θ,z) to formulate the flow under consideration. Following [24], it is assumed that the flow is subject to periodic acceleration in the z-direction. The schematic diagram of the catheterized artery is presented in Fig 1. Since the flow in the rigid catheterized tube is assumed to be unsteady and uni-directional, therefore we write velocity field in core and peripheral region as follows: The momentum equation in this case simplifies to where G(t) is periodic acceleration in z-direction. We divide our domain kR r R as kR r R 1 (a core region) and R 1 r R (a periphery region). For the present problem, the fluid in the core is an Oldroyd-B fluid while in the periphery region it is a Newtonian fluid. Thus the constitutive equation for the shear stress in periphery region is where μ 2 is the viscosity of the Newtonian fluid. The general form of the constitutive equation for shear stress in the core region is [26] 1 þ l 1 @ @t The pressure gradient can be put in the form [28]: Here A 0 is the systolic and A 1 is the amplitude (diastolic) components of the pressure gradient, ω p = 2πf p is the circular frequency and f p is the pulse rate frequency. Following [29], we can write.
in which A g is the amplitude, f b is the frequency (ω b = 2πf b ) and ϕ is the lead angle of G(t) with respect to the heart action. Eliminating S rz between (10) and (12), we get in the core region and for the periphery region.
where subscript denotes differentiation with respect to the indicated variable. The boundary and initial conditions for the present flow situation are [26]: The volumetric Flow rate and shear stress at the wall of the tube are respectively given by The resistance to flow or impedance experienced by flowing blood at any cross-section is [28]:

Dimensionless formulation
We are interested in numerical solution of Eqs 15 and 16 subject to conditions 17, for this we first normalize these equations by defining [29] In terms of new variable the momentum equations after dropping the bars takes the following form The dimensionless boundary and initial conditions become Similarly the volume flow rate and the shear stress in dimensionless form are The various parameter appearing in Eqs 21-25 are defined as follows [28]:

Analytical Solution
Eqs 21 and 22 subject to the initial and boundary condition given in Eq 23 can be solved analytically by employing Laplace transform. To this end, let us denote U 1 (r,s) and U 2 (r,s), respectively, the Laplace transform of u 1 (r,t) and u 2 (r,t). Then Eqs 21 and 22 reduce to following ordinary differential equations in transformed domain where The boundary and initial conditions in the transform domain read The condition (35) is established by taking Laplace transform of both sides of interface condition in Eq 17 and using the u 1t = u 2t at t = 0. The differential Eqs 27 and 28 can be readily solved in terms of Bessel's functions to get the expressions of U 1 (r,s) and U 2 (r,s) as The values of arbitrary constants c 1-4 can be easily calculated by implementing the boundary conditions (33)- (35). Due to the complicated nature of the solution in transformed domain, analytical inversion is difficult. Therefore, we have relied on the numerical inversion using Mathematica Package NumericalInversion.m. The command "Stehfest" is used to get the numerical values of the inverted solution.

Numerical Solution
Eqs 21 and 22 subject to the boundary condition given in Eq 23 are solved numerically using the finite difference method [30,31]. The uniformly distributed discrete points in radial direction are defined as r i = (i−1)Δr, (i = 1,2,. . ..,N c + 1) such that r ðN c þ1Þ ¼ b and r i = (i−(N c +1))Δr, i = (N c + 1, N c + 2,. . .,N + 1) such that r N+1 = 1, where Δr is the increment in the radial direction. Similarly we define t j = (j−1)Δt, (j = 1,2,. . ..) as discrete time points with Δt indicating the small time increment. For problem under consideration we have chosen Δr = 0.025 and Δt = 0.0001. This choice of Δr and Δt yields results convergent up to order 10 −7 . Let us denote the discretized value of u k (k = 1,2) at (r i ,t j ) by u j k i . In this notation, finite difference formulas for first and second order derivatives read @u k @r ffi u j k iþ1 À u j k iÀ1 2 Dr ¼ u k r ; For the time derivative, we define the approximation: Utilizing (38) and (39), Eqs 21 and 22 may be transformed to the following difference equations: The boundary and initial conditions in the discretize form can be written as follows: The partial differential equations governing the flow are solved both numerically and analytically. The details of the stability of numerical scheme used by us are given in the book by Hoffman, in which it is mentioned that the stability of explicit finite scheme depends on time and spatial step sizes. The numerical solution obtained by us fulfill all these requirements on spatial and temporal step sizes. Following Hoffman [32], we have chosen Δx = 0.025, Δt = 0.0001. This choice of Δx and Δt yield results convergent up to order 10 −7 . However it is observed during the computations that the numerical results are much sensitive to the choice of involved parameters. Thus for the present problem the stability of the numerical solution cannot be guaranteed only by full filling the requirements of the temporal and spatial step sizes. In fact with the present choice of spatial and temporal step sizes a stable numerical solution is possible only if material parameters λ 1 and λ 2 are such that 0.1 < λ 1 < 0.5 and 0.1 < λ 2 < 0.5.

Results and Discussion
In this section, we are interested to analyze the graphical results of velocity, flow rate, wall shear stress and resistance to flow results for different values of non-dimensional variables namely the catheter radius ratio (k), amplitude of pulsatile pressure gradient (e), relaxation and retardation time (λ 1 and λ 2 ) and generalized Womersley frequency parameter (α). In the present analysis the height of interface is assumed independent of the involved parameters and therefore it is not determined as a part of the problem. Instead we have taken it as a constant. The primary motivation for this assumption is based the evidence that RBC is migrate away from the wall leaving a cell-depleted layer near the wall. In such situation one can model blood flow in two stages (i) a peripheral layer modeled as a Newtonian liquid and (ii) a core region modeled as a non-Newtonian fluid. The thickness of both layer is assumed constant. Based on this fact several authors for instance Maji and Nair [33], Massoudi and Phouc [29], Ikbal et al. [31], Bugliarello and sevilla [34], Sulkla et al. [35], Akay an Kaye [36], Pralhad and Schultz [37], Sajid et al. [38] and Srivastava and Sexena [39] have used two-layered fluid model with constant peripheral thickness to study blood flow through arteries. We have followed these studies in our analysis. However, we have also shown the results for different interface positions in Figures (4, 10 and 11). The physiologically relevant values of the above mentioned parameters is listed in Tables 1 and 2. The radius ratio k, characterizing the radius of the catheter is assumed to be in the range 0.1-0.5. The length of separation point between peripheral layer and the core region is represented by β. In this way, we get the thickness of the peripheral layer as 1−β. In the present analysis, β are taken as 0.9. Fig 2 illustrates the profiles of velocity at different instants of time within a single cardiac cycle calculated through both numerical and analytical solutions. A pleasing agreement between both solutions is observed. This agreement also demonstrates the validity of our numerical scheme. Fig 2 further indicates an increase in velocity when t increases from 0.1 to 0.2 in systolic phase. However in diastolic phase i.e. at t = 0.3 a decreasing trend in velocity is observed. The increasing behavior of velocity in systolic phase while opposite trend in diastolic phase is clearly due to the pulsatile pressure gradient produced by the pumping action of heart. Fig 3 shows the effects of the catheterization at t = 0.3 on the velocity distribution in the two-fluid model of blood flow. It is noticed that at a given instant of time t = 0.3 with β = 0.9 and the increasing values of the catheter radius ratio k, the velocity decreases considerably. It means that due to the insertion of catheter, the magnitude of the velocity reduces significantly.
The velocity distribution for different values of the interface position β is shown in Fig 4 for k = 0.1. It is found that for a given value of k, the velocity increases significantly with decreasing β, whereas the behavior is opposite when the width of the peripheral layer decreases (i.e., when the value of β decreases). We have shown the effects of relaxation and retardation times on the velocity profile in Fig  5. It is readily observed that the role of λ 1 here is to increase the magnitude of velocity while λ 2 counters the effects of λ 1 and it decreases the magnitude of velocity. Generally it is clear that the flow rate increases as t increases from 0 to 2.3 and then achieves its steady state condition. It means that flow rate fluctuates around its mean value with constant frequency and amplitude after achieving steady state condition. It is further observed from this figure that the amplitude of oscillations in flow rate increases (decreases) by increasing λ 1 (λ 2 ). Similarly as expected that amplitude of oscillations in flow rate increases by increasing the amplitude of pulsatile pressure gradient (Fig 9). Moreover it also evident from Similarly, the variation of dimensionless flow rate for different values of catheter radius ratio k and interface position β is shown in Fig 10. It is noticed that the flow rate decreases by increasing both catheter radius and peripheral layer thickness.
In pulsatile flow, the non-dimensional wall shear stress τ w can be calculated from Eq 25. The variation of wall shear stress in a cycle of oscillation for different values of amplitude e and peripheral thickness β with k = 0.5 is shown in Fig 11. The wall shear stress fluctuates around Numerical and Analytical Study of Two-Layered Unsteady Blood Flow through Catheterized Artery its mean value as time increases. It is also noticed that amplitude of the wall shear is proportional to the magnitude of e. Its value increases with the increase in the magnitude of e while it shows opposite trend with the increase of peripheral radius length β. Fig 12 shows the variation of wall shear stress in a cycle of oscillation for different values of λ 1 , λ 2 and k with β = 0.9. This figure indicates that the magnitude of wall shear stress decreases by increasing λ 1 . while it is suppressed by increasing λ 2 . Moreover it is also observed through Fig 12 that the magnitude of the wall shear increases with an increase in the radius of catheter.
The longitudinal impedance (Λ) of the artery is calculated using Eq 19 and its variation during a flow cycle for different values of amplitude k with e = 0.5 and β = 0.9 is illustrated in Fig  13. It is observed that these profiles follow an opposite trend as compare to the flow rate profiles which is expected from the formula of resistance to flow given by expression (19). It is concluded from these profiles that the magnitude of resistance to flow increases with the increase of the size of catheter.

Conclusion
Two-layer mathematical model is developed for pulsatile flow of blood with the effects of body acceleration through a catheterized artery using Oldroyd-B fluid model in the core region and Newtonian model in the peripheral region. The analysis is based on the solution of a linear partial differential equation corresponding to each region. Assuming the continuity of velocity and no-slip conditions, the governing equation are solved in each region by employing finite difference technique. The validity of numerical solution is checked by comparing it with the analytical solution obtained through Laplace transform. An excellent agreement between both solutions is observed. After validation, numerical solution is further utilized to obtain various flow quantities. It is found that the resistance to flow or impendence, velocity, flow rate and wall shear stress are greatly affected by blood rheology, catheter radius and amplitude of pulsatile pressure gradient. More precisely a small increment in the catheter radius significantly decreases the magnitude of blood velocity and as a result of that a decrement in the flow rate is observed. As a consequence of decrease in velocity the wall shear stress and resistance to flow also increase considerably. Another important observation is that such quantitative measurement varies in magnitude by changing the rheological parameters of blood. In fact for a given catheter radius the flow rate of blood increases by increasing dimensionless relaxation time and as a result the resistance to flow decreases. Such observation may have certain implications for instance, the rheology of blood can be controlled to maintain the same flow in an artery with catheter as that in an artery without catheter. The results obtained through present computations may be related to hypertension. In normal circumstances/conditions, the amplitude of pressure gradient maintain the sufficient flow of blood in arteries for normal function of the organs. However, in situation where the person/human is subject to periodic body acceleration the results are quite different. In such situation the superimposed periodic oscillation due to body acceleration increase the amplitude of flow rate produced by pulsatile pressure gradient. This increase in the amplitude of the flow rate may results in hypertension, nausea, loss of vision and headache. On the contrary, the magnitude of the flow rate decreases with increasing Numerical and Analytical Study of Two-Layered Unsteady Blood Flow through Catheterized Artery  the catheter radius. In such situation the impedance/ resistance to flow rate increases resulting in hypotension. Although the significance of the present study is limited but still this study throws some light on the flow behavior of blood. In fact this study presents the effects interface position, non-Newtonian nature and an inserted catheter on flow features. Analytical as well a numerical solution is also presented in this study. The analytical solution may be used as a benchmarking for more sophisticated 3D CFD analysis accounting pointing by the worthy reviewer.

Author Contributions
Conceptualization: AZ NA MS TH.