Abstract
A mathematical model of dispersed bioparticle-blood flow through the stenosed coronary artery under the pulsatile boundary conditions is proposed. Blood is assumed to be an incompressible non-Newtonian fluid and its flow is considered as turbulence described by the Reynolds-averaged Navier-Stokes equations. Bioparticles are assumed to be spherical shape with the same density as blood, and their translation and rotational motions are governed by Newtonian equations. Impact of particle movement on the blood velocity, the pressure distribution, and the wall shear stress distribution in three different severity degrees of stenosis including 25%, 50%, and 75% are investigated through the numerical simulation using ANSYS 18.2. Increasing degree of stenosis severity results in higher values of the pressure drop and wall shear stresses. The higher level of bioparticle motion directly varies with the pressure drop and wall shear stress. The area of coronary artery with higher density of bioparticles also presents the higher wall shear stress.
1. 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–12]. The flow phenomena includes asymmetric flow, unsteady laminar-to-turbulent flow [13–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–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 fluid-structural 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 two-way 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 three-dimensional 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.
2. Mathematical Model
The computational domain is modelled by using two coordinate systems, i.e., an Eulerian frame Ω(X1;X2;X3) for the fluid flow and a Lagrangian frame ΩL(x1;x2;x3) 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.
2.1. 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 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 and represent zero shear viscosity and infinite shear viscosity, is constant shape parameter, is consistent index, is time constant, and .
2.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 operatorwhere 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 (the dependent variables) of the continuity equation are decomposed into Navier-Stokes equations. Consequently, the mean part and the fluctuation part are defined aswhere 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.
In (9), representing the nonlinear convective transport due to turbulent velocity fluctuations is defined bywhere is component of the mean velocity, , and is eddy viscosity defined byfor ; and , , , , and is the modulus of the mean rate-of-strain tensor.
Two blending functions in (7) and in (12) are defined bywith and is the distance to the next surface. In (6),where and = 8. Other four parameters , and in the SST - model are determined bywhich are the relationships that transform the constants of the original - model into the constants of the SST - model.
The external force 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 is the dimensionless relative particle Reynolds number defined byThe velocity component of the particle, , is the arithmetic average nodal values on the face based on Green–Gauss node-based method defined bywhere 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.,forTo be more specific, the remaining constant model parameters introduced in this subsection are declared as follows:
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, = , are governed by the following Newtonian equations:where denotes the particle position vector, and are, respectively, Lagrangian velocity and angular velocity of the particle, is the particle mass, and the moment of rotational inertia given by 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 are the sum of various forces including drag force , Saffman force , and the Magnus or rotational lift force acting on the particle by carrier fluid:where , = 2.594, and is the deformation tensor along the path of particles defined byRegarding (30), is the th projected particle surface area, is the relative fluid-particle velocity, and denotes the rotational lift coefficient. Taking into account the effect of the rotational Reynolds number and the particle Reynold number , the coefficients and are assigned asThe term in (25) represents the contact forces due to the particle-particle collision and the particle-wall collision; i.e., . Based on the Spring-Dashpot Collision Law [24], is given bywhere K is the elastic collision coefficient, is the overlap of any two particles, is the damping coefficient, is the relative velocity, and is a unit vector. For the position vectors , and the radii and of the particles and , we haveThe damping coefficient in (34) depending on the mass loss in the collision process and the collision time scale are defined bywhere , 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 is calculated in the same way as .
2.4. 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 are calculated by the following Fourier series:where denotes the angular frequency defined by 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 asforwhere is the pulsatile flow rate and is the inflow surface area. The turbulence kinetic energy and the specific dissipation rate at the inflow arewith the percentage of turbulence intensity and turbulence length scale for .
On the outflow boundaries including , , and , the boundary conditions are set to and the normal pressure gradient field is corresponding to . 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.
3. 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.
3.1. 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 Reynolds-averaged 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.
(a)
(b)
(a)
(b)
3.2. 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.
(a) LAD with stenosis
(b) The LCA-LAD axial line
(c) The LCA-LCX axial line
(d) Pulsatile pressure at the LCA inflow surface
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 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 the particle-particle and the particle-wall collisions on the blood flow pattern.
(a) A part of LCA
(b) Streamline of blood flow
(c) Velocity field of blood flow
(d) Particle deposition patterns
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.
(a) Cutting plane
(b) Blood velocity in each plane
Distribution of systolic pressure along two axial lines including the LCA-LAD line in as shown in Figure 6(b) and the LCA-LCX line as shown in Figure 6(c) is investigated in the first and the fifth cardiac cycle. Regarding the fluid flow model and the dispersed particle-fluid flow model, blood pressure at the peak systole at t = 0.65 s and t = 3.86 s along the two axial lines with three different degrees of stenosis severity are plotted in Figures 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. The results indicate that 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 cycles (t = 0.65 s and t = 3.86 s). The dispersed particle-fluid flow model gives high variation of pressure and wall shear stress, especially in the LCA to the proximal LAD. Figures 11 and 12 describe the effect of particle motion on the pressure distribution in two investigated lines, the LCA-LAD axial line and the LCA-LCX axial line as shown in Figures 6(b) and 6(c).
(a) Model with no particle
(b) Model with particle motion
(a) Pressure distribution with particle deposition patterns
(b) Wall shear stress distribution with particle deposition patterns
(a) 25%
(b) 50%
(c) 75%
(a) 25%
(b) 50%
(c) 75%
(a) 25%
(b) 50%
(c) 75%
(a) 25%
(b) 50%
(c) 75%
Wall shear stress along the first half of the LCA connecting to the LAD is higher as shown in Figure 13. At the first cardiac cycle t = 0.65 s, the maximum wall shear stress in the proximal 75% stenosis of the LAD is about 225 Pa in the fluid flow model and 275 Pa in the dispersed particle-fluid flow model. In addition, at the fifth cardiac cycle t = 3.86 s, this maximum wall shear stress is about 240 Pa in the fluid flow model and 250 Pa and the dispersed particle-fluid flow model. In Figure 14, presents variation of wall shear stress along the LCA connecting to the normal LCX in the first and the fifth cardiac cycle. It indicates that particle motion in the carrier 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).
4. 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 Reynolds-averaged 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Supplementary Materials
The symptoms of coronary heart disease depend on the severity of the blockage in the 2 main coronary arteries, the left main, and the right coronary arteries. Understanding blood flow around the blockage is thus necessary for a bypass surgery. Two supplementary figures show the coronary system of human arteries that was constructed using 1000 images of computed tomography scans of the human coronary system. The system consists of the base of the aorta connecting with the normal right coronary artery (RCA) and the left coronary (LCA) with appearance of a LAD stenosis located at 5 mm from the aorta-LCA connection. One supplementary figure presents the normal coronary system. Another one is the system with LAD stenosis of 75%. In this study, we considered only the left main coronary artery with two branches, the left anterior descending artery (LAD), and the circumflex artery (LCX), to investigate the effect of the severity of coronary artery stenosis. (Supplementary Materials)