Numerical Simulation of Dispersed Particle-Blood Flow in the Stenosed Coronary Arteries

1Department of Mathematics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand 2Centre of Excellence in Mathematics, Commission of Higher Education (CHE), Bangkok 10400, Thailand 3Department of Mathematics, College of Information and Communication Technology, Rangsit University, PathumThani 12000, Thailand 4Department of Mathematics, Faculty of Science, Maejo University, Chiang Mai 50290, Thailand 5School of Electrical Engineering, Computing and Mathematical Sciences, Curtin University, Perth, WA 6845, Australia


Introduction
Atherosclerosis is a disease narrowing a coronary artery due to plaque buildup.Generally, there is no symptom until it severely narrows the artery causing serious problems including heart attack, stroke, or even death.Critical information of blood flow in the stenotic coronary arteries is a principle factor of the development and progression of atherosclerosis.Figure 1 presents an angiogram of a critical proximal left anterior descending artery (LAD) in a patient with Wellens' syndrome.Atherosclerosis is often associated with some forms of abnormal blood flow in the blocked coronary arteries.Dealing with the pathogenesis of coronary artery diseases (CAD), various practical treatment of CAD including drug delivery, stent replacement, and coronary artery bypass grafting (CABG) have been developed through a number of in vivo and in vitro experiments.Due to the high rate of stent and graft failures, development of vascular drug delivery, one of the key rubrics of targeted therapeutics and nanodevices, becomes more and more important [1].Recently, several drug delivery approaches are undergoing clinical testing and medical industry development.
Over decades, many researchers have carried out experimental models and computational simulations to explore the flow phenomena in the stenotic arteries in order to optimize medical methods of treatment.Due to the difficulty and limitation in determining the critical flow conditions for both in vivo and in vitro experiments, the exact mechanisms involving these treatments are not well understood.Thus, mathematical modelling and numerical simulation are chosen to be a better alternative to analyze the problem.Complex phenomena of blood flow in arteries subject to various physiological conditions has been extensively analyzed using various mathematical models [2][3][4][5][6][7][8][9][10][11][12].The flow phenomena includes asymmetric flow, unsteady laminar-to-turbulent flow [13][14][15][16].The governing equations include the Navier-Stoke equations and the continuity equations subjected to the selected inlet velocity, no-slip condition at the artery wall, and stress-free condition at the outflow surface.The unsteady flow also is characterized by a high pressure drop and high wall shear stress (WSS) in the stenotic artery [1,[17][18][19].
For the requirement of simulation, the constant density of blood is assumed to be equal to 1055 kg m −3 .The vessel wall may be considered as rigid or elastic and the Power law non-Newtonian model or Carreau-Yasuda non-Newtonian model are normally applied to describe the viscosity of blood.Considering the blood flow as turbulent, the two-equation turbulence model, so-called the standard k- model or the k- model, has been commonly employed for the analysis.Young (1973) studied blood flow through an occluded tube under a pulsatile pressure gradient.Mazumdar et al. (1996) investigated the unsteady Newtonian blood flow through a stenosed artery and observed that the pressure gradient attains the maximum at the throat of stenosis and decreases with an increase of hematocrit parameter.Sanyal and Maiti (1998) proposed a mathematical model of blood flow in the artery with mild stenosis.They reported that the pressure gradient increases with an increase in hematocrit value where there is higher value in systolic and lower value in diastolic pressure.Deplano and Siouffi (1999) performed experimental and numerical study of pulsatile flows through a 75% severity stenosis to determine the wall shear stress temporal evolution downstream from the stenosis.The result shows high wall shear stress values of about 120 Pa (or dyn cm −2 ) during the cardiac cycle at the throat, and low values downstream from the stenosis of about -2.5 Pa (or dyn cm −2 ). Lee and Xu (2002) studied blood flow through a rigid mild stenosed tube.Wiwatanapataphee and Wu (2012) investigated the unsteady non-Newtonian blood flow through the real right coronary artery bypass graft system under the real pulsatile condition.They reported that the existence and intensity of a stenosis in the artery have significant effect on the blood flow behaviour.Gupta and Agrawal (2015) simulated the blood flow passing through an irregular stenotic descending aorta using a Finite Volume method.The results demonstrate that the formation of wall shear stress in the stenotic region by the irregular stenosis model is much complex than by regular stenosis.High oscillation of wall shear stress appears behind the irregular stenosis in which the Reynold number (Re) is between 130 and 540.
As the mechanical interaction between blood and arterial wall has an important role in the propagation of pressure wave from the heart to the whole body, many researchers have investigated the fluid-structural response to pulsatile blood flow through a stenosed vessel subject to various physiological conditions.Chan (2006) investigated the fluidstructural response to the pulsatile non-Newtonian blood flow through an axisymmetric stenosed vessel using ANSYS.The solid model was set to have isotropic elastic properties.The Fluid-Structural Interaction (FSI) coupling was twoway and iterative.It is found that interaction between vessel wall and the blood gives reasonable results.Due to the wall expansion, the axial velocity decreases and the recirculation effect of the flow increases.Torii et al. ( 2009) studied blood flow in the deformable cerebral artery using an FSI model and reported that the maximum wall shear stress tends to decrease when the blood flow impinges strongly on the wall.
The consideration of a non-Newtonian behaviour of blood in small arteries gives more relevant results and prediction and understanding of pressure distribution and wall shear stress have significant importance in disease diagnosis and surgical planning.The models with full threedimensional FSI problem increase the computational work.Hence, development of the computational framework built upon image-based CFD and discrete particle dynamics modelling is a big challenge.Bernad et al. (2013) investigated the particle motion in coronary serial stenoses to analyze the hemodynamic significance of three serial stenoses, named ST1, ST2, and ST3, in the right coronary artery (RCA) constructed from multislice computerized tomography images.Blood was assumed to be an incompressible Newtonian fluid and the artery walls was rigid and impermeable.Results illustrate that pressure drop increases with an increase of percentage stenosis.During the systolic phase, the pressure drop is higher about 32.84 mmHg and 36.78 mmHg for the stenosis ST1 and ST3 during the systolic phase while it is lower about 4.62 mmHg and 4.81 mmHg during the diastolic phase.For stenosis ST2, the pressure drop is not significant during the systolic and diastolic phases.They reported that wall shear stress distribution has a close reflection of the outline of the stenosis and the formation of recirculation zone.Range of wall shear stress varies from 7 to 262 Pa.Three intense regions of wall shear stress appear downstream at each stenosis, and its value is low in the recirculation zone.Mukherjee and Shadden (2017) studied embolic particle dynamics and transport through swirling chaotic flow structures of various vasculature beds.The results show the complex interplay of particle inertia, fluid-particle density ratio, and wall collisions, with chaotic flow structures, which render the overall motion of the particles to be nontrivially dispersive in nature.These researches motivate the present study to deal with the particle motion in a pulsatile blood flow through a stenotic artery.However, a few attempts have been made to study the flow phenomena in the stenotic artery in which the fluid-particle interaction, the particle-particle collision, and the particle-wall collision are considered.
The aim of this study is to investigate the hemodynamic parameter, the pressure distribution, and wall shear stress in a stenosed artery using Reynolds-averaged Navier-Stokes equations for turbulence fluid modelling and Newtonian equations for the translation and rotational motion of the bioparticles.Using the pulsatile boundary conditions based on a physiological waveforms of flow velocity and blood pressure, the results are compared with those obtained from the model with no particle to highlight the role of particle motion in the turbulence model of blood flow in the left coronary artery (LCA) connecting to the stenosed left anterior descending artery (LAD) (Figure 16) and the normal left circumflex coronary artery (LCX) (Figure 15).The rest of this paper is organized as follows.In the following sections, a mathematical model describing the turbulence flow of fluid and movement of dispersed phase under pulsatile condition is taken into presentation.The governing equations of dispersed particle-blood flow are theoretically presented in Section 2. Section 3 concerns numerical investigations to analyze velocity field, pressure distribution, and wall shear stress distribution along two investigated axial lines with respect to different degrees of stenosis severity.At the end of this paper, some discussion and conclusion are given in Section 4.

Mathematical Model
The computational domain is modelled by using two coordinate systems, i.e., an Eulerian frame Ω(X 1 ;X 2 ;X 3 ) for the fluid flow and a Lagrangian frame Ω  (x 1 ;x 2 ;x 3 ) for particle movement.The carrier fluid is assumed to behave as a non-Newtonian incompressible fluid.Both phases of dispersed particles and blood exchange momentum are allowed to have the fluid-particle interaction.In this study, we focus on two-way coupling model in which fluid phase influences particulate phase via drag and turbulence, and particulate phase influences fluid phase via source terms of mass and momentum.The dispersion of particles due to turbulence in the fluid phase is predicted using the stochastic tracking (random walk) model including the effect of instantaneous turbulent velocity fluctuations on the particle trajectories through the use of stochastic methods [20].The fluid flow influences the particle trajectories and the dispersed particles with the particle-particle collision and the particle-wall collision have a significant effect on the flow turbulence.

Laminar Flow of Non-Newtonian Incompressible Fluid.
The flow of a non-Newtonian incompressible fluid is governed by the continuity equation and the Navier-Stokes equations as follows: where   and  ext  denote, respectively, velocity component of fluid and the external force,  is the fluid pressure,   is the fluid density, and  is the viscosity of the non-Newtonian fluid based on the Carreau model; i.e., where  0 and  ∞ represent zero shear viscosity and infinite shear viscosity,  is constant shape parameter,  is consistent index,  is time constant, and   = (  /  +   /  )/2.

Turbulence Flow of Non-Newtonian Incompressible Fluid.
The turbulence flow of a non-Newtonian incompressible fluid which is formulated on a moving reference frame is governed by the mean continuity equation, the Reynolds-averaged Navier-Stokes equations, and Menter's SST - model [21] as follows: with the operator where   denote mean velocity component of fluid and  is the fluid viscosity based on the Carreau model (3).In (5), where   is the Reynolds stress term.According to the concept of Reynolds decomposition, the individual instantaneous velocity component   ,  = 1, 2, 3 (the dependent variables) of the continuity equation are decomposed into Navier-Stokes equations.Consequently, the mean part   and the fluctuation part    are defined as where    , ( = 1, 2, 3) are constant functions of time obtained from stochastic model and their random value is kept constant over an time interval given by the characteristic lifetime of the eddies.
Two blending functions  1 in ( 7) and  2 in ( 12) are defined by with  +  = max{2 −1 ,2  −1 (/  )(/  ), 10 −10 } and  is the distance to the next surface.In (6), where  * ∞ = 0.09 and   = 8.Other four parameters , ,   , and   in the SST - model are determined by which are the relationships that transform the constants of the original - model into the constants of the SST - model.The external force  ext  in ( 5) is formulated by examining the change in momentum of the solid particle passing through each control volume computed as a function of the particle mass flow rate ṁ and the time step Δ; i.e., where   is the particle density and Re ] is the dimensionless relative particle Reynolds number defined by The velocity component of the th particle, V   , is the arithmetic average nodal values on the face based on Green-Gauss node-based method defined by where  face represents number of nodes on face and Ṽ  is nodal value on the face.The drag coefficient,   , is based on the spherical drag law [22,23]; i.e., for To be more specific, the remaining constant model parameters introduced in this subsection are declared as follows: where x  denotes the particle position vector, k  and Ω  are, respectively, Lagrangian velocity and angular velocity of the  particle,   is the particle mass, and   the moment of rotational inertia given by   = (/60)   5  for the sphere particle with diameter of   and density of   .The right side term of (26) is the resulting torque depending on the rotational drag coefficient   and the relative particle-fluid angular velocity Υ defined by F   are the sum of various forces including drag force F  , Saffman force F  , and the Magnus or rotational lift force F  acting on the particle by carrier fluid: where   =    2  (18) −1 ⋅ 24(   V ) −1 ,  = 2.594, and   is the deformation tensor along the path of particles defined by Regarding (30),    is the th projected particle surface area, V is the relative fluid-particle velocity, and   denotes the rotational lift coefficient.Taking into account the effect of the rotational Reynolds number Re  and the particle Reynold number Re ] , the coefficients   and   are assigned as = 6.45 √Re  + 32.1 Re  ; The term F  in (25) represents the contact forces due to the particle-particle collision and the particle-wall collision; i.e., F  = F  + F  .Based on the Spring-Dashpot Collision Law [24], F  is given by where K is the elastic collision coefficient,  is the overlap of any two particles,  is the damping coefficient, k  is the relative velocity, and e  is a unit vector.For the position vectors x  , x  and the radii   and   of the particles  and , we have The damping coefficient  in (34) depending on the mass loss in the collision process   and the collision time scale   are defined by where ,   , and   denote the damper restriction coefficient, the mass of particle , and the mass of particle , respectively.The contact force due to the particle-wall collision F  is calculated in the same way as F  .

Initial and Boundary
Conditions.Initialization of the particle properties in the domain is not required and there is no interface boundary condition for point-mass particles.To complete defining the boundary value problem, the following boundary conditions are required.In this study, the pulsatile mass flow inlet and the pulse pressure outlet are used.The wave forms of the mass flow rate () and the pressure (t) are calculated by the following Fourier series: where  denotes the angular frequency defined by  = 2/ with a cardiac period , and all values of the parameters are given by Wiwatanapataphee et al. [16].On the inflow boundary, we have the pulsatile velocity as for where () is the pulsatile flow rate and  is the inflow surface area.The turbulence kinetic energy and the specific dissipation rate at the inflow are with the percentage of turbulence intensity  = 0.16(  ()/) 1/8 and turbulence length scale  = 0.07 for  = 1.13 1/2 .
On the outflow boundaries Γ LCA including Γ LCA1 , Γ LCA2 , and Γ LCA3 , the boundary conditions are set to  ⋅ n = ()n and the normal pressure gradient field is corresponding to /n = 0.For an inert and point-mass particle, interface condition between fluid and solid particle is not required.Furthermore, it is also assumed that the wall has an infinite radius and zero velocity, and no-slip condition is applied on the arterial wall.

Numerical Investigation
This section presents numerical simulation of particle movement and turbulence non-Newtonian fluid flow.The simulation of discrete particles is carried out using the discrete element method.With this method, particle trajectories are calculated through the simulation domain in Lagrangian reference.To reduce time-consuming simulation, a limited number of representative trajectories is calculated.The simulation of the continuous phase is carried out using the Finite Volume method.By this method, the mean continuity equation, the Reynolds-averaged Navier-Stokes equations, and the Menter's SST - model are solved in Eulerian reference frame.

Validation Study.
Studying the turbulence fluid-solid (two-phase) flow in the coronary artery with stenosis requires a reliable model that can fully describe the complex phenomena occurring in the artery with nonlinear response.The first task is undertaken for evaluating the suitability of the mathematical model using ANSYS 18.2.The simulation of the turbulence two-phase flow in the normal curved tube and 75% stenosis-curved tube was setup using the discrete model coupled with the mean continuity equation, the Reynoldsaveraged Navier-Stokes equations, and Menter's SST - model for a turbulent viscous incompressible non-Newtonian fluid.The computational domains of both tubes are displayed in Figure 2. Five complete pulses of pressure and flow velocity were used in each simulation.The results as shown in Figures 3 and 4 indicate that pressure drop presents in the tube with restricted area.The model with particle movement gives variation of pressure in the area occupied by the particle.These results show that our proposed model can capture important phenomena in the flow channel without and with restricted area.

Fluid-Particle Flow through the Human Left Coronary
Artery with Stenosis.To begin the numerical simulation, a 3D computational domain of a human left coronary artery with its branches including the LAD and LCX is firstly constructed by replicating the multislice computerized tomography (CT) image as shown in Figure 5. Taking into account the flow direction, there are a single inlet at the beginning of LCA and three outlets at the tail ends of LAD and LCX.
Next, the effect of coronary stenosis on the flow of fluid (blood) and discrete particles (bioparticles) is taken into investigation.In this simulation, 414 particles are tracked to determine the behaviour of the dispersed phase.The distribution of the dispersed phase, bioparticles, in the midst of blood inside the arterial vessel is closely focused when there exists a stenosis.Three different degrees of stenosis severity including 25%, 50%, and 75% are assigned at the proximal part of LAD as can be seen in Figure 6(a).After grid independent test, we obtained suitable domain mesh with 10 boundary layers with the size of the first cell 5 × 10 −3 mm for four cases including the normal artery and the artery with 25%, 50%, and 75% stenosis details as shown in Table 1.
The flow pattern, the pressure distribution, and the wall shear stress (WSS) distribution are analyzed.To investigate the effect of particle motion on pressure distribution and wall shear stress distribution, 500 bioparticles are injected into the LCA inflow surface at the same speed of the pulsatile velocity in the first cardiac cycle, t = 0.2 s.Assigning values of model parameters as in Table 2 and running the simulation corresponding to various numerical options as in Table 3, the numerical results along the left coronary artery connecting to the critical LAD with various-degree stenosis and the normal LCX simulated using ANSYS 18.2 for the turbulent dispersed particle-fluid flow are obtained in Tables 2 and 3.
Focusing on the proximal part of LCA as shown in Figure 7(a), the results obtained from the turbulent flow      Green -Gauss node -based the particle-particle and the particle-wall collisions on the blood flow pattern.
In the fluid flow model, small turbulent flow appears in the transition area connecting the LCA with its branches as shown in Figures 7(b) and 7(c).In the dispersed particle-fluid flow model, high turbulent flow appears in the LCA region and the particles continue to proportionally flow from LCA to LAD and LCX as displayed in Figure 7(d).Considering the blood velocity inside LCA, LAD, and LCX, the vector plots in three different planes are exhibited in Figure 8.The higher degree stenosis reduces number of particles flowing through the downstream LAD and increases number of particles flowing into the downstream LCX.The model with 0%, 25%, 50%, and 75% degree of stenosis severity on the proximal LAD allows 12.0%, 6.4%, 3.4%, and 1% of all particles flowing through the downstream LAD, respectively.Particles are almost completely blocked in the LAD with critical 75% stenosis.11 and 12, respectively.The results indicate that a pressure drop across the 75% stenosis is significant at the peak systole t = 0.65 s.In the fluid flow model at the peak systole in a cardiac cycle, the pressure drop is about 32 mmHg.In the dispersed particle, fluid flow model, the pressure drop at the peak systole in the first and the fifth cardiac cycles is significantly different due to particle deposition patterns in the stenotic area.Higher level of particle motion makes more pressure drop.At t = 0.65 s, a pressure drop is about 55 mmHg and 40 mmHg at t = 3.86 s, respectively.
The level of particle motion varies with time.The highest level is at the peak systole of the first cardiac cycle t = 0.65 s.To show the effect of particle motion on the blood pressure and the wall shear stress, the results at the peak systole t = 3.86 s obtained from the turbulent flow model and the turbulent dispersed particle-fluid flow model are compared in Figures 9-14   fluid as shown in Figure 10 has significant effect on the wall shear stress.High wall shear stress occurs in the area with high particle concentration (particle cluster).

Discussion and Conclusion
This paper presents the mathematical model of the dispersed bioparticle-blood flow in the left coronary artery (LCA) with its branches including the LAD and LCX.The combination of the mean continuity equation, the Reynoldsaveraged Navier-Stokes equations, and the Menter's SST k- models is employed to investigate the turbulence flow of blood, the non-Newtonian incompressible fluid.Describing the movement of dispersed particle phase, the Newtonian equations are used to examine the translation and rotational motion of bioparticles.Running the simulation of two-phase flow inside LCA, the 3D computational domain together with initial and boundary conditions are necessary.The dispersed phase flow, the pressure distribution and the wall shear stress distribution are analyzed corresponding to three different stenosis intensities of 25%, 50%, and 75% at the proximal LAD.The results demonstrate a significant effect of turbulence blood flow on the particle trajectories and a high impact of particle motion with the particle-particle and the particle-wall collisions on the blood flow pattern.The coronary artery with critical 75% stenosis generates a sudden drop of pressure with high wall shear stress around the stenosis cite.Higher degree of stenosis gives higher values of the pressure drop and wall shear stresses.Pressure distribution and wall shear stresses are plotted when the flow is at a maximum at the peak systole in the first and the fifth cardiac cycle.Pressure drop is significantly different due to particle deposition patterns in the stenotic area at the peak systole.Higher level of particle motion makes more pressure drop and has significant effect on the wall shear stress that occurring in the area with higher particle concentration.

( 23 ) 2 . 3 .
Movement of Dispersed Particle Phase.For the dispersed particle phase, we assume that all particles are spherical, and heat and mass transfer are omitted.All particles are treated as point masses and inert type.The translation and rotational motions of the th particle,  = 1, ...,   , are governed by the following Newtonian equations:

Figure 2 :
Figure 2: Computational domain of validation study.

Figure 3 :
Figure 3: Blood pressure obtained from the model with no particle (solid line) and the model with particles (dotted line) in a normal curved tube at two different peak systoles: (a)  = 0.65 ; (b)  = 3.86s.

Figure 4 :Figure 5 :
Figure 4: Blood pressure obtained from the model with no particle (solid line) and the model with particles (dotted line) in the 75% stenosis tube at two different peak systoles: (a) t = 0.65 s; (b) t = 3.86s.

Table 2 :− 3
Model parameters.Particle mass flow rate ( ṁ ) 5 .5 0x1 0 −8 kg s −1 Particle time step (Δ) 1 × 10 −6 s model and the turbulent dispersed particle-fluid flow model demonstrate the effect of turbulence blood flow on the particle trajectories and the impact of particle motion with International Journal of Differential Equations Pulsatile pressure at the LCA inflow surface

Figure 6 :
Figure 6: Stenosis conditions at the proximal LAD, two investigated lines along the LCA connecting to the LAD and the LCX, and pulsatile pressure with two investigated times at the peak systole in the first and the fifth cardiac cycles.
(a) A part of LCA (b) Streamline of blood flow (c) Velocity field of blood flow (d) Particle deposition patterns

Figure 7 :
Figure 7: Dispersed phase flow in the LCA connecting to the LAD and LCX.
. The results indicate that the coronary artery with critical 75% stenosis generates a sudden drop of pressure with high

Figure 15 :
Figure 15: The system of coronary arteries including the base of the aorta and the normal left and the normal right arteries.

Figure 16 :
Figure 16: The system of coronary arteries including the base of the aorta and the 75% stenosed left artery and the normal right artery.