Numerical Solution of Blood Flow and Mass Transport in an Elastic Tube with Multiple Stenoses

The simultaneous effect of flexible wall and multiple stenoses on the flow and mass transfer of blood is investigated through numerical computation and simulations. The solution is obtained using the Marker and Cell technique on an axisymmetric model of Newtonian blood flow. The results compare favorably with physical observations where the pulsatile boundary condition and double stenoses result in a higher pressure drop across the stenoses. The streamlines, the iso-concentration lines, the Sherwood number, and the mass concentration variations along the entire wall segment provide a comprehensive analysis of the mass transport characteristics. The double stenoses and pulsatile inlet conditions increase the number of recirculation regions and effect a higher mass transfer rate at the throat, whereby more mass is expected to accumulate and cause further stenosis.


Introduction
Caro et al. [1] postulated that atherosclerosis, which is a narrowing of the artery as a result of plaque build-up may occur due to shear-dependent mass transfer mechanism between blood cholesterol and the arterial wall. Cholesterol exists in blood in the form of low density lipoproteins (LDLs) whose deposition along the walls of the artery is a key step in atherogenesis, which would lead to stenosis. Stenosis can affect the velocity of blood flowing through the artery, affecting blood pressure, collapsing the heart, which could in turn lead to disastrous consequences. us, an understanding of the behavior of local mass transport in arterial stenosis is important in the study of the formation and development of atherosclerotic lesions for appropriate assessment on the possible correlation between the site of atherosclerotic lesions and the pattern of mass transport.
Ethier [2] carried out computational modelling of mass transfer and studied its links to atherosclerosis. Other studies on mass transport and fluid flow in stenotic arteries of axisymmetric and asymmetric models have been carried out by [3][4][5][6]. In these studies, the arterial wall was considered as rigid and the artery is assumed to have single mild stenosis, in which the geometry of the stenosis is represented by the usual cosine curve along with a restriction that the ratio of the severity of stenosis and the radius of the artery is very small. In reality, this is not the case where in many medical situations, the patient is found to have multiple stenoses in the same arterial segment.
Investigations on the effect of multiple stenoses on blood flow have been carried out amongst others by [7][8][9][10]. ese studies showed that from both experimental results and theoretical calculations, the total effect of a series of noncritical stenoses is approximately equal to the sum of their individual effects where they can be critical and produce symptoms of arterial insufficiency. e flow energy loss due to the presence of the stenoses, which is directly related to the pressure drop across them, increases with the number of stenoses but is not strongly dependent on the spacing between them. e authors of [11][12][13][14][15][16][17][18][19] have also investigated blood flow through multiple stenoses; however, these studies have not considered the mass transfer.
Another aspect to be considered in arterial blood flow is the cyclic nature of the heart pump which creates pulsatile conditions in the arteries, giving rise to unsteady flow. It is observed that most CFD models of arterial hemodynamics make the simplifying assumptions of rigid walls and fully developed inlet velocities (cf. [13][14][15][16][17][18][19]). But the arteries are not rigid tubes. ey adapt to varying flow conditions by enlarging or shrinking. All of these physiological conditions make the modelling and consequently the solution to be almost impossible to be obtained analytically and challenging computationally. Nandakumar and Anand [20] studied steady and pulsatile flow of blood through a channel with single as well as double stenoses on the assumption that the pulsations of flow are damped in the small vessels; thus the flow is effectively steady in the capillaries and the veins while Liu and Tang [21] investigated the influence of distal stenosis on blood flow through curved arteries with two stenoses. But again, these studies on pulsatile flow have also not considered the mass transfer. In another study, Layek et al. [22] investigated the effect of multiple stenoses on the flow of Newtonian fluid in a rigid tube and opined that the disturbance created by the constrictions is mainly concentrated at the downstream of the last constriction. Considering the flow of Newtonian fluid in a two-dimensional channel having a single constriction, Layek and Midya [23] concluded that the maximum stress and the length of the recirculation region associated with two shear layers of the constriction do increase with the increasing area reduction of the constriction. ey further concluded that the flow-field separates after the symmetry breaking bifurcation, and the symmetry of the flow depends on Reynolds' number and the height of the constriction. e flow of a fluid having hematocrit-dependent viscosity past a tube with partially overlapped constriction has been investigated by Layek et al. [24]. ey observed that the peak value of the wall shear stress decreases with increasing haematocrit parameter while a reverse trend is observed for the flow separation region. ey also opined that the deformability of the wall does reduce the wall shear stress as compared to the rigid wall case. All these studies [22][23][24] ignored the flow pulsatility and/or consideration of multiple constrictions, and the mass transfer as well which plays a pivotal role in the genesis and evolution of atherosclerosis. Based on the gap established above, with regard to studies involving mass transfer, the following work seeks to analyze the flow and mass transfer characteristics of pulsatile blood flow through an artery with double stenoses. e fluid considered is Newtonian in an axisymmetric setting, while the pair of stenoses vary in severities, lengths, and distances between them. e equation for stenoses is given in an algebraic form which could represent both moderate and severe stenoses instead of the usual cosine function which could only describe mild stenosis.
e objective of the present study lies in the consideration of the transport of mass as well as momentum together through a tube with a flexible wall, resembling the flexibility of the artery in the presence of double stenoses. e flow pulsatility cannot be ruled out from the present investigation.

Formulation of the Problem
We consider a fully developed two-dimensional axisymmetric flow of an incompressible Newtonian fluid of density ρ in a tube. e relevant equations of motions in vector forms are the continuity, momentum, and mass as follows: with D/Dt is the material derivative, V � (u, 0, w) where u and w are the radial and axial velocity components, respectively, p is the pressure, μ is the constant viscosity, C is the mass concentration, and D m is the coefficient of mass diffusion.
In the cylindrical coordinate system, the corresponding equations (1)-(3) are written in a conservative form as follows: e schematic diagram for the double stenoses is given in Figure 1, where r � R(z, t) is the radius of the artery in the stenotic region and R 0 is the radius of the artery in the nonstenotic regions. δ 1 , δ 2 are the critical heights of the first and second stenosis respectively; l 0 is the inlet segment, l 02 is the distance between stenoses, l 01 , l 03 are the lengths of stenoses, and L is the length of the arterial segment under consideration.
e equations describing the stenoses are given by the following: e time-variant parameter a 1 (t) is given by a 1 (t) � 1 + k cos(ωt) with k representing the amplitude parameter and ω the angular frequency is given by ω � 2πf p , f p being the pulse frequency and d � l 0 + l 01 + l 02 . To the best of our knowledge equation (8) is the first equation to address double stenoses without any control on the severity of stenoses which has not been considered before.

Boundary Conditions
where U is the cross-sectional average velocity of the fluid and C s is a constant.

Solution Procedure
e solution procedure involves the nondimensionalization, radial coordinate transformation, and the finite-difference Marker and Cell method (MAC) initially proposed by Harlow and Welch [25]. Sarifuddin et al. [26,27] and Mustapha et al. [15,16] have used the method to solve blood flow problems.

Nondimensionalization of the Equations
e nondimensional variables and parameters introduced are as follows: Using (13), equations (4)- (7) have their respective nondimensional forms as follows (omitting bar): BioMed Research International e boundary conditions (9)-(12) reduce to their respective dimensionless forms: at z � 0, where Re is the Reynolds number, Sc is the Schmidt number and α is the Womersley number defined as follows:

Finite-Difference
Method. e solution procedure consists of discretization of the governing equations, combining the discretized forms of the momentum and continuity equations to obtain the Poisson equation for pressure, the successive overrelaxation (SOR) method, and the pressure and velocity corrections. e schematic computational domain is given in Figure 2. e velocities and pressure are calculated at different locations of the control volume, as indicated in Figure 3. e difference equations are derived at three distinct cells, each corresponding to the continuity, axial and radial momentum equations. e discretization of the time derivative terms is based on the first-order accurate two-level forward time differencing formula, while the convective terms in the momentum equations are discretized with a hybrid formula consisting of central differencing and second-order upwinding scheme (cf. Courant et al. [28]). e diffusive terms are discretized using second-order accurate three- where n refers to time and Δt is the time increment. e length and width of the (i, j) n cell of the control volume are represented by Δz and Δx, respectively. e discretized version of the continuity equation (24) at the (i, j) cell is where . Here (z li , x lj ) and (z i , x j ) represent the respective coordinates of the center of the cell and the cell faces as shown in Figure 3, while w t and w b stand for w-velocities at the top and bottom middle positions of the control volume of the continuity equation. e momentum equations (24) and (25) are written in the following forms: where conw n i,j , conu n i,j , diffw n i,j and diffu n i,j are the finite-difference representation of convective and diffusive terms of the axial and radial momentum at the nth time level.
BioMed Research International e Poisson equation for pressure is derived from equations (31)-(33) which takes the final forms: where D n+1 i,j represents the discretized form of the divergence of the velocity field at the (i, j) cell and the expressions for A i,j , B i,j , . . . , H i,j , S i,j are the same as Mustapha et al. [15,16].
e Poisson equation (34) for pressure is then solved using the successive overrelaxation (SOR) method to obtain the intermediate pressure field. e increment Δx is chosen to be 0.025 along x, while Δz � 0.1 along z. Δt is chosen to be equal to or less than a prescribed stability criterion as depicted in Figure 4, where here c is taken to be 0.05. (cf. [25,29]). e number of iterations is limited up to 10. e pressure and velocities then go through a correction stage to achieve  better accuracy. e process is described in [15,16,26,27]. When the velocity field has been obtained, the mass concentration is calculated from the respective discretized versions of equation (26) with the relevant boundary conditions (equations (27)- (30)). e values chosen for k, α, Sc are 0.05, 2, and 3, respectively.

Results and Discussions
e influence of the pulsatile inlet is reflected in Figure 5 where the pressure drop in this case is higher than the one generated with the parabolic inlet. It decreases with increasing Re with a strong linear correlation between them in both cases of parabolic and pulsatile conditions (cf. Sarifuddin et al. [26]). Further, the pressure drop is seen to increase with the number of stenoses, which agrees with the experimental study of Talukder et al. [8].
e behavior of the axial and the radial velocity at the narrowest points (z � 10 and z � 19) for different Re are shown in Figures 6(a), 6(b), 7(a), and 7(b). e axial velocity has positive values, and it is noted that the parabolic case results in higher velocity. It is also observed that the axial velocity near the wall increases with increasing Re; however, there is a cross over at x � 0.65 and x � 0.73 for z � 10 and z � 19, respectively. Figures 7(a) and 7(b) show that the radial velocity corresponding to the pulsatile inlet assume positive values at z � 10 except the value on the wall at (x � 1) while negative values are observed at the narrowest point (z � 19) and the flow velocity increases with increasing Re near the centerline while it is reduced near the wall with increasing Re. Both figures reveal that the velocities in the case of the parabolic inlet are negative and substantially less than that with the pulsatile inlet. Figure 8(a) exhibits the axial velocity profiles at different locations of the stenosed arterial segment at Re � 300 for both parabolic and pulsatile inlet conditions. At (z � 19), the velocity in the parabolic case is higher than in the pulsatile case. A backflow occurs in the pulsatile case at the downstream of the narrowest point (z � 19) near the wall. e curves decrease from their individual maximum at the axis as one moves away from it and finally they approach a minimum value (zero) on the wall surface. Note that the curves of the axial velocity at (z � 10) and at (z � 15) are coincident. e axial velocity at the critical height of the second stenosis (z � 19) is considerably higher than that of the first stenosis (z � 10). Figure 8(b) shows that the radial velocity has positive values everywhere at (z � 5), (z � 10), and (z � 25), excluding the position on the wall. At (z � 15) and the narrowest point (z � 19), the radial velocity is observed to have all negative values. e nonzero values of the radial velocity near the wall clearly reflect the influence of the radial motion of the arterial wall in the pulsatile case. Figure 9(a) exhibits the distribution of the wall shear stress (wss) along the arterial segment for different Re considering pulsatile as well as parabolic inlet conditions. e results show that wss for both the parabolic and pulsatile inlet conditions attain their peaks at the critical heights of the stenoses (z � 10, 19). It is observed that BioMed Research International separation occurs (negative values of wss) only at the downstream of the second stenosis for parabolic inlet condition. In the pulsatile case, the separation zone occurs between the two stenoses with multiple separation regions at the downstream of the second stenosis. en, wss starts to increase slowly towards the wall surface (reattachment point).
e effect of different severities on wss is depicted in Figure 9(b). In the pulsatile case, when the two stenoses have the same severities (δ 1 � δ 2 � 0.2) and, (δ 1 � δ 2 � 0.4), flow separation occurs at two specific places at the downstream of the first and the second stenoses with different peak values. In the case of stenoses with different severities (δ 1 � 0.2, δ 2 � 0.4) and (δ 1 � 0.4, δ 2 � 0.2) a larger separated region is formed at the downstream of the more severe stenoses (cf. Johnston and Kilpatrick [12]). A smaller separation region is produced in the case of pulsatile inlet condition and the peak wss is much higher than the parabolic inlet condition. Figure 9(c) determines the effects of the length of stenoses on wss. Peak wall shear stress decreases with increasing the length of stenosis but it increases with the gap between stenosis and at this position, there is a potential that plaque would rupture whereas, at the low shear stress position, atherosclerotic development may be induced. ese phenomena of separation and reattachment are due to the adverse pressure in these regions and are believed to be responsible for the malfunctioning of the cardiovascular system having atherosclerotic plaque. Figures 10(a)-10(d) show the instantaneous patterns of streamlines governing the flow of blood through the stenoses in case of (δ 1 � 0.4, δ 2 � 0.2) for both parabolic and pulsatile inlet conditions at Re � 300 and Re � 500. In the parabolic case, only one recirculation zone developed between the two stenoses where separation occurs at z � 11 (c.f Figure 9(b)). In the case of the pulsatile inlet, a multiple recirculation region is noted in between the two stenoses (separation point z � 11).
us, an increase in Re and a consideration of the pulsatile flow increase the number of the recirculation region.   increase as the flow gets accelerated towards the throat leading to the increase of solute concentration. It is also observed that the concentration at the throat (z � 19) is much higher for pulsatile flow than the parabolic one. 12(c). e Sherwood number defined by Sh D � 2R 0 c l /D m ΔC where c l is the local mass flux to the arterial wall and 2R 0 is the inlet diameter of the artery. It is observed that Sh D increases with increasing Re while Sh D distribution appreciably changes specifically at the throat, between the two stenoses and at the downstream position. e highest mass transfer is experienced at the upstream, while the minimum value occurs at the downstream of the stenoses. Note that pulsatile flow increases the Sherwood number and that it is much higher in the case of double stenosis. e iso-concentration lines considering pulsatile as well as parabolic inlet conditions based on (δ 1 � 0.2, δ 2 � 0.4) are displayed in Figures 13(a) and 13(b). e iso-concentration lines for parabolic and pulsatile inlets have different distributions with multiple recirculation regions nearby the downstream of the more severe stenosis in the pulsatile case. e general trend of the iso-concentration lines is that they move away from the inlet region towards the upstream of the stenosis and correspondingly impair the mass transport in this region, while they adhere to the outline of the stenosis at both the upstream and downstream ends. At this region of low wss (compare with Figure 9), cholesterol may tend to accumulate and cause more severe stenosis. is observation conforms with Schneiderman et al. [30].

Conclusion
e hemodynamics of the pulsatile flow and the transport of mass in an arterial segment having a couple of stenoses have been studied in relation to the distensibility of the vessel wall. Predicted results show that the pulsatile inlet and double stenoses with varying severity affect the flow characteristics significantly, especially the development of the recirculation zone and the peak value of the wall shear stress. It is also predicted that the concentration at the throat (z � 19) is much higher for pulsatile flow than the parabolic inlet condition. Moreover, the pair of stenoses contributes much to the mass concentration than the case of single stenosis. e mass flux to the arterial wall (Sherwood number) does increase with the increasing values of Re and here too, mass flux increases with the flow pulsatility and the presence of double stenoses. At the downstream, cholesterol may tend to accumulate and causes more severe stenosis. For severe stenoses, the peak value of the wall shear stress is higher in the pulsatile flow case and the isoconcentration lines show more recirculation regions nearby the downstream end and their lengths are longer. In conclusion, the results presented agree well with physical observations and provide an insight into the link between atherosclerosis, stenosis, and the pattern of mass transport.
ough the detailed knowledge of the dynamical variables is possible and provides useful elements, the mechanism of influence of the haemodynamical factors in the arterial disease is not clear. e characteristics of the red cells must be taken into consideration by including a shear-dependent viscosity in the diffusion terms in time-dependent flows highlighting the scope of further work. All these mechanical and biochemical aspects related to the biofluid dynamics are of some importance and demand further investigation. A great deal of work is needed to establish the rheological parameters for the physiological values and to understand the connection of the issues with biological facts.
Data Availability e data on blood flow parameters used for analysis and validation purposes are from previously reported studies and datasets, which have been cited.

Disclosure
To the best of the authors' knowledge, the content of this work is correct . Each author has contributed in part to the problem formulated, numerical computations and simulations, analysis of results, and manuscript writing.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.