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

Authors: Mongkol Kaewbumrung; Somsak Orankitjaroen; Pichit Boonkrong; Buraskorn Nuntadilok; Benchawan Wiwatanapataphee
Journal: International Journal of Differential Equations
 (2018)

---

## 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.


---

## Body

## 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. Figure1 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.Figure 1
Tight, severe stenosis (95%) of the proximal LAD in a patient with Wellens’ warning.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) (Figure16) 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:(1)∂ui∂Xi=0,(2)ρf∂ui∂t+ρfuj∂ui∂Xj=-∂∂Xip+∂∂Xjμ∂ui∂Xj+Fiextwhere ui and Fiext denote, respectively, velocity component of fluid and the external force, p is the fluid pressure, ρf is the fluid density, and μ is the viscosity of the non-Newtonian fluid based on the Carreau model; i.e.,(3)μ=μ∞+μ0-μ∞1+λ2∑i,jSij2an-1/a,where μ0 and μ∞ represent zero shear viscosity and infinite shear viscosity, a is constant shape parameter, n is consistent index, λ is time constant, and Sij=(∂ui/∂Xj+∂uj/∂Xi)/2.


### 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 SSTk–ω model [21] as follows:(4)∂u-i∂Xi=0,(5)ρfDu-iDt=-∂p-∂Xi+∂τ^ij∂Xj+Fiext,(6)ρfDkDt=-τij∂u-i∂Xj-ρfβ∗ωk+∂∂Xjμ+μtσk∂k∂Xj,(7)ρfDωDt=-αμtτij∂u-i∂Xj-ρfβω2+∂∂Xjμ+μtσω∂ω∂Xj+2ρf1-f1σω,2ω∂k∂Xj∂ω∂Xj,with the operator(8)DDt≡∂∂t+u-j∂∂Xj,where u-i denote mean velocity component of fluid and μ is the fluid viscosity based on the Carreau model (3). In (5),(9)τ^ij=2μSij-13δij∂uk∂Xk+τij,where τij is the Reynolds stress term. According to the concept of Reynolds decomposition, the individual instantaneous velocity component ui,i=1,2,3 (the dependent variables) of the continuity equation are decomposed into Navier-Stokes equations. Consequently, the mean part u-i and the fluctuation part ui′ are defined as(10)ui=u-i+ui′,where ui′,(i=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.In (9), τij representing the nonlinear convective transport due to turbulent velocity fluctuations is defined by(11)τij=2μtS-ij-13δij∂u-k∂Xk-23ρfkδij,where uk¯  is component of the mean velocity, S-ij=∂u-i/∂Xj+∂u-j/∂Xi/2, and μt is eddy viscosity defined by(12)μt=ρfkω1max⁡1/α∗,S-f2/a1ω,for α∗=α∞∗α0∗+Ret/Rk1+Ret/Rk-1; Ret=ρfkμω-1 and Rk=6, α0∗=βi/3, βi=0.072, a1=0.31, and S-≡2S-ijS-ij1/2 is the modulus of the mean rate-of-strain tensor.Two blending functionsf1 in (7) and f2 in (12) are defined by(13)f1=tanh⁡ϕ14;ϕ1=min⁡max⁡kβ∞∗ωy,500μρy2ω,4ρkσω,2Dω+y2;(14)f2=tanh⁡ϕ22;ϕ2=max⁡2kβ∞∗ωy,500μρy2ωwith Dω+=max⁡2ρσω,2-1ω-1∂k/∂Xj∂ω/∂Xj,10-10 and y is the distance to the next surface. In (6),(15)β∗=β∞∗4/15+Ret/Rβ41+Ret/Rβ4,where β∞∗=0.09 and Rβ = 8. Other four parameters α,β,σk, and σω in the SST k-ω model are determined by(16)α=f1α-1+1-f1α-2,σk=f1σk,1+1-f1σk,2,(17)β=f1β1+1-f1β2,σω=f1σω,1+1-f1σω,2,which are the relationships that transform the constants of the original k-ω model into the constants of the SST k-ω model.The external forceFiext 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 m˙p and the time step Δt; i.e.,(18)Fiext=∑p=1np18μCDRev24ρp2ui-uipm˙pΔt,where ρp is the particle density and Reν is the dimensionless relative particle Reynolds number defined by(19)Rev=ρfdpu-vpμThe velocity component of the pth particle, vip, is the arithmetic average nodal values on the face based on Green–Gauss node-based method defined by(20)vip=1Nface∑p=1Nfacev~ip,where Nface represents number of nodes on face and v~ip is nodal value on the face.The drag coefficient,CD, is based on the spherical drag law [22, 23]; i.e.,(21)CD=α1+α2Rev+α3Rev3for(22)α1,α2,α3=0,24,00.0≤Rev<0.13.690,22.73,0.09030.1≤Rev<11.222,29.1667,-3.88891.0≤Rev<100.6167,46.50,-116.6710≤Rev<1000.3644,98.33,-2778100≤Rev<10000.357,148.62,-475001000≤Rev<50000.46,-490.546,5787005000≤Rev<100000.5191,-1662.5,5416700Rev≥10000.To be more specific, the remaining constant model parameters introduced in this subsection are declared as follows:(23)α-1=59,β1=340,σk,1=0.85,σω,1=0.5,α-2=0.44,β2=0.0828,σk,2=1,σω,2=0.856.


### 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 thepth particle, p = 1,...,np, are governed by the following Newtonian equations:(24)dxpdt=vp,(25)mpdvpdt=Fpf+Fpc(26)IpdΩpdt=ρp2dp25CωΥ,where xp denotes the particle position vector, vp and Ωp are, respectively, Lagrangian velocity and angular velocity of the p particle, mp is the particle mass, and Ip the moment of rotational inertia given by Ip=π/60ρpdp5 for the sphere particle with diameter of dp and density of ρp. The right side term of (26) is the resulting torque depending on the rotational drag coefficient Cω and the relative particle-fluid angular velocity Υ defined by(27)Υ=12∇×u-Ωp;Fpf are the sum of various forces including drag force FD, Saffman force FS, and the Magnus or rotational lift force FRL acting on the particle by carrier fluid:(28)FD=u-vpτp,(29)Fs=2κυρfdijρpdpdlkdkl1/4u-vp,(30)FRL=ρf2AipCRLVΥV×Υ,where τp=ρpd2p18μ-1·24CDRev-1, κ = 2.594, and djj is the deformation tensor along the path of particles defined by(31)dij=12∂vi∂xj+∂vj∂xiRegarding (30), Aip is the ith projected particle surface area, V is the relative fluid-particle velocity, and CRL denotes the rotational lift coefficient. Taking into account the effect of the rotational Reynolds number Reω and the particle Reynold number Reν, the coefficients CRL and Cω are assigned as(32)CRL=0.45+ReωRev-0.45exp⁡-0.05684Reω0.4Rev0.3,(33)Cω=6.45Reω+32.1Reω;Reω=ρfΩdp24μThe term Fpc in (25) represents the contact forces due to the particle-particle collision and the particle-wall collision; i.e., Fpc=Fpq+Fpw. Based on the Spring-Dashpot Collision Law [24], Fpq is given by(34)Fpq=Kδ-γvpq·epqepq,whereKis the elastic collision coefficient, δ is the overlap of any two particles, γ is the damping coefficient, vpq is the relative velocity, and epq is a unit vector. For the position vectors xp, xq and the radii rp and rq of the particles p and q, we have(35)δ=xq-xp-rp+rq,(36)vpq=vq-vp,(37)epq=xq-xpxq-xpThe damping coefficient γ in (34) depending on the mass loss in the collision process mpq and the collision time scale tc are defined by(38)γ=2mpqln⁡ηtc,(39)mpq=mpmqmp+mq,(40)tc=π2+ln2⁡ηmpqKwhere η,mp, and mq denote the damper restriction coefficient, the mass of particle p, and the mass of particle q, respectively. The contact force due to the particle-wall collision Fpw is calculated in the same way as Fpq.


### 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 rateQ(t) and the pressure p(t) are calculated by the following Fourier series:(41)Qt=Q-+∑n=14anQcos⁡nwt+bnQsin⁡nwt,(42)pt=p-+∑n=14anpcos⁡nwt+bnpsin⁡nwt,where w denotes the angular frequency defined by w=2π/T with a cardiac period T, and all values of the parameters are given by Wiwatanapataphee et al. [16]. On the inflow boundary, we have the pulsatile velocity as(43)u=u-=Unt;vp=Unt=0.2s,for(44)Unt=QtAwhere Q(t) is the pulsatile flow rate and A is the inflow surface area. The turbulence kinetic energy and the specific dissipation rate at the inflow are(45)k=32Untς2,(46)ω=0.1643k1.5Lwith the percentage of turbulence intensity ς=0.16(ρUn(t)D/μ)1/8 and turbulence length scale L=0.07D for D=1.13A1/2.On the outflow boundariesΓLCA including ΓLCA1, ΓLCA2, and ΓLCA3, the boundary conditions are set to σ·n=p(t)n and the normal pressure gradient field is corresponding to ∂p/∂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.


## 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:(1)∂ui∂Xi=0,(2)ρf∂ui∂t+ρfuj∂ui∂Xj=-∂∂Xip+∂∂Xjμ∂ui∂Xj+Fiextwhere ui and Fiext denote, respectively, velocity component of fluid and the external force, p is the fluid pressure, ρf is the fluid density, and μ is the viscosity of the non-Newtonian fluid based on the Carreau model; i.e.,(3)μ=μ∞+μ0-μ∞1+λ2∑i,jSij2an-1/a,where μ0 and μ∞ represent zero shear viscosity and infinite shear viscosity, a is constant shape parameter, n is consistent index, λ is time constant, and Sij=(∂ui/∂Xj+∂uj/∂Xi)/2.


## 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 SSTk–ω model [21] as follows:(4)∂u-i∂Xi=0,(5)ρfDu-iDt=-∂p-∂Xi+∂τ^ij∂Xj+Fiext,(6)ρfDkDt=-τij∂u-i∂Xj-ρfβ∗ωk+∂∂Xjμ+μtσk∂k∂Xj,(7)ρfDωDt=-αμtτij∂u-i∂Xj-ρfβω2+∂∂Xjμ+μtσω∂ω∂Xj+2ρf1-f1σω,2ω∂k∂Xj∂ω∂Xj,with the operator(8)DDt≡∂∂t+u-j∂∂Xj,where u-i denote mean velocity component of fluid and μ is the fluid viscosity based on the Carreau model (3). In (5),(9)τ^ij=2μSij-13δij∂uk∂Xk+τij,where τij is the Reynolds stress term. According to the concept of Reynolds decomposition, the individual instantaneous velocity component ui,i=1,2,3 (the dependent variables) of the continuity equation are decomposed into Navier-Stokes equations. Consequently, the mean part u-i and the fluctuation part ui′ are defined as(10)ui=u-i+ui′,where ui′,(i=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.In (9), τij representing the nonlinear convective transport due to turbulent velocity fluctuations is defined by(11)τij=2μtS-ij-13δij∂u-k∂Xk-23ρfkδij,where uk¯  is component of the mean velocity, S-ij=∂u-i/∂Xj+∂u-j/∂Xi/2, and μt is eddy viscosity defined by(12)μt=ρfkω1max⁡1/α∗,S-f2/a1ω,for α∗=α∞∗α0∗+Ret/Rk1+Ret/Rk-1; Ret=ρfkμω-1 and Rk=6, α0∗=βi/3, βi=0.072, a1=0.31, and S-≡2S-ijS-ij1/2 is the modulus of the mean rate-of-strain tensor.Two blending functionsf1 in (7) and f2 in (12) are defined by(13)f1=tanh⁡ϕ14;ϕ1=min⁡max⁡kβ∞∗ωy,500μρy2ω,4ρkσω,2Dω+y2;(14)f2=tanh⁡ϕ22;ϕ2=max⁡2kβ∞∗ωy,500μρy2ωwith Dω+=max⁡2ρσω,2-1ω-1∂k/∂Xj∂ω/∂Xj,10-10 and y is the distance to the next surface. In (6),(15)β∗=β∞∗4/15+Ret/Rβ41+Ret/Rβ4,where β∞∗=0.09 and Rβ = 8. Other four parameters α,β,σk, and σω in the SST k-ω model are determined by(16)α=f1α-1+1-f1α-2,σk=f1σk,1+1-f1σk,2,(17)β=f1β1+1-f1β2,σω=f1σω,1+1-f1σω,2,which are the relationships that transform the constants of the original k-ω model into the constants of the SST k-ω model.The external forceFiext 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 m˙p and the time step Δt; i.e.,(18)Fiext=∑p=1np18μCDRev24ρp2ui-uipm˙pΔt,where ρp is the particle density and Reν is the dimensionless relative particle Reynolds number defined by(19)Rev=ρfdpu-vpμThe velocity component of the pth particle, vip, is the arithmetic average nodal values on the face based on Green–Gauss node-based method defined by(20)vip=1Nface∑p=1Nfacev~ip,where Nface represents number of nodes on face and v~ip is nodal value on the face.The drag coefficient,CD, is based on the spherical drag law [22, 23]; i.e.,(21)CD=α1+α2Rev+α3Rev3for(22)α1,α2,α3=0,24,00.0≤Rev<0.13.690,22.73,0.09030.1≤Rev<11.222,29.1667,-3.88891.0≤Rev<100.6167,46.50,-116.6710≤Rev<1000.3644,98.33,-2778100≤Rev<10000.357,148.62,-475001000≤Rev<50000.46,-490.546,5787005000≤Rev<100000.5191,-1662.5,5416700Rev≥10000.To be more specific, the remaining constant model parameters introduced in this subsection are declared as follows:(23)α-1=59,β1=340,σk,1=0.85,σω,1=0.5,α-2=0.44,β2=0.0828,σk,2=1,σω,2=0.856.


## 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 thepth particle, p = 1,...,np, are governed by the following Newtonian equations:(24)dxpdt=vp,(25)mpdvpdt=Fpf+Fpc(26)IpdΩpdt=ρp2dp25CωΥ,where xp denotes the particle position vector, vp and Ωp are, respectively, Lagrangian velocity and angular velocity of the p particle, mp is the particle mass, and Ip the moment of rotational inertia given by Ip=π/60ρpdp5 for the sphere particle with diameter of dp and density of ρp. The right side term of (26) is the resulting torque depending on the rotational drag coefficient Cω and the relative particle-fluid angular velocity Υ defined by(27)Υ=12∇×u-Ωp;Fpf are the sum of various forces including drag force FD, Saffman force FS, and the Magnus or rotational lift force FRL acting on the particle by carrier fluid:(28)FD=u-vpτp,(29)Fs=2κυρfdijρpdpdlkdkl1/4u-vp,(30)FRL=ρf2AipCRLVΥV×Υ,where τp=ρpd2p18μ-1·24CDRev-1, κ = 2.594, and djj is the deformation tensor along the path of particles defined by(31)dij=12∂vi∂xj+∂vj∂xiRegarding (30), Aip is the ith projected particle surface area, V is the relative fluid-particle velocity, and CRL denotes the rotational lift coefficient. Taking into account the effect of the rotational Reynolds number Reω and the particle Reynold number Reν, the coefficients CRL and Cω are assigned as(32)CRL=0.45+ReωRev-0.45exp⁡-0.05684Reω0.4Rev0.3,(33)Cω=6.45Reω+32.1Reω;Reω=ρfΩdp24μThe term Fpc in (25) represents the contact forces due to the particle-particle collision and the particle-wall collision; i.e., Fpc=Fpq+Fpw. Based on the Spring-Dashpot Collision Law [24], Fpq is given by(34)Fpq=Kδ-γvpq·epqepq,whereKis the elastic collision coefficient, δ is the overlap of any two particles, γ is the damping coefficient, vpq is the relative velocity, and epq is a unit vector. For the position vectors xp, xq and the radii rp and rq of the particles p and q, we have(35)δ=xq-xp-rp+rq,(36)vpq=vq-vp,(37)epq=xq-xpxq-xpThe damping coefficient γ in (34) depending on the mass loss in the collision process mpq and the collision time scale tc are defined by(38)γ=2mpqln⁡ηtc,(39)mpq=mpmqmp+mq,(40)tc=π2+ln2⁡ηmpqKwhere η,mp, and mq denote the damper restriction coefficient, the mass of particle p, and the mass of particle q, respectively. The contact force due to the particle-wall collision Fpw is calculated in the same way as Fpq.


## 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 rateQ(t) and the pressure p(t) are calculated by the following Fourier series:(41)Qt=Q-+∑n=14anQcos⁡nwt+bnQsin⁡nwt,(42)pt=p-+∑n=14anpcos⁡nwt+bnpsin⁡nwt,where w denotes the angular frequency defined by w=2π/T with a cardiac period T, and all values of the parameters are given by Wiwatanapataphee et al. [16]. On the inflow boundary, we have the pulsatile velocity as(43)u=u-=Unt;vp=Unt=0.2s,for(44)Unt=QtAwhere Q(t) is the pulsatile flow rate and A is the inflow surface area. The turbulence kinetic energy and the specific dissipation rate at the inflow are(45)k=32Untς2,(46)ω=0.1643k1.5Lwith the percentage of turbulence intensity ς=0.16(ρUn(t)D/μ)1/8 and turbulence length scale L=0.07D for D=1.13A1/2.On the outflow boundariesΓLCA including ΓLCA1, ΓLCA2, and ΓLCA3, the boundary conditions are set to σ·n=p(t)n and the normal pressure gradient field is corresponding to ∂p/∂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.


## 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 SSTk-ω 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 SSTk-ω 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.Figure 2
Computational domain of validation study.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)t = 0.65 s; (b) t = 3.86s.
(a)
t = 0.65   s (b)
t = 3.86   sFigure 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.65s; (b)t= 3.86s.
(a)
t = 0.65   s (b)
t = 3.86   s


### 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 Figure5. 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.Figure 5
Computational domain of the LCA and its branches including the LAD and the 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 Figure6(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.Table 1
Mesh information.
% Stenosis Number of elements 0 % 314,150 25 % 318,897 50 % 350,265 75 % 444,764Figure 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)
LAD with stenosis (b)
The LCA-LAD axial line (c)
The LCA-LCX axial line (d)
Pulsatile pressure at the LCA inflow surfaceThe 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.Table 2
Model parameters.
Parameters value units Blood density (ρ) 1050 kg m−3 Zero shear viscosity (μ0) 0.56 g cm−1 s−1 Infinite shear viscosity (μ∞) 0.0345 g cm−1 s−1 Time constant (λ) 3.3130 s Consistency index (n) 0.3568 - Shape parameter (a) 2 - Particle diameter (dp) 100 μm Particle density (ρp) 1050 kg m−3 Particle mass flow rate (m˙p) 5.50 x 10−8 kg s−1 Particle time step (Δt) 1×10-6 sTable 3
Numerical setting.
Numerical Parameters Setting Solver Pressure - based Pressure – Velocity Coupling Coupled Multigrid AMG Under – relaxation factors 0.75 - momentum 0.75 - pressure 0.80 - turbulent kinetic energy 0.80 - specific dissipation rate 0.90 - discrete phase source Spatial discretization Second order upwind - momentum Second order upwind - turbulent kinetic energy Second order upwind - specific dissipation rate Second order upwind - transient formulation Gradient Green - Gauss node - basedFocusing on the proximal part of LCA as shown in Figure7(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.Figure 7
Dispersed phase flow in the LCA connecting to the LAD and LCX.
(a)
A part of LCA (b)
Streamline of blood flow (c)
Velocity field of blood flow (d)
Particle deposition patternsIn the fluid flow model, small turbulent flow appears in the transition area connecting the LCA with its branches as shown in Figures7(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.Figure 8
Vector plot of blood velocity at the peak systole at the timet= 0.65 s in three different planes.
(a)
Cutting plane (b)
Blood velocity in each planeDistribution of systolic pressure along two axial lines including the LCA-LAD line in as shown in Figure6(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 att= 0.65 s andt= 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 systolet= 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. Att= 0.65 s, a pressure drop is about 55 mmHg and 40 mmHg att= 3.86 s, respectively.The level of particle motion varies with time. The highest level is at the peak systole of the first cardiac cyclet= 0.65 s. To show the effect of particle motion on the blood pressure and the wall shear stress, the results at the peak systolet= 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 andt= 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).Figure 9
Systolic pressure att= 0.65 s in the artery model with 75% intensity of stenosis at the proximal LAD obtained from two different models: (a) model with no particle; (b) model with particle motion.
(a)
Model with no particle (b)
Model with particle motionFigure 10
Pressure distribution and wall shear stress distribution with particle deposition patterns in the artery model with 75% intensity of stenosis at the proximal LAD at the peak systole in the fifth cardiac cycle att= 3.86 s.
(a)
Pressure distribution with particle deposition patterns (b)
Wall shear stress distribution with particle deposition patternsFigure 11
Blood pressure at the peak systolet= 0.65 s andt= 3.86 s along the LCA-LAD axial line obtained from the model with no particle (solid line) and the model with particles (dotted line) obtained from the domain with different intensity of stenosis at the proximal LAD: (a) 25%; (b) 50%; (c) 75%.
(a)
25% (b)
50% (c)
75%Figure 12
Blood pressure at the peak systolet= 0.65 s andt= 3.86 s along the LCA-LCX axial line obtained from the model with no particle (solid line) and the model with particles (dotted line) obtained from the domain with different intensity of stenosis at the proximal LAD: (a) 25%; (b) 50%; (c) 75%.
(a)
25% (b)
50% (c)
75%Figure 13
Wall shear stress at the peak systolet= 0.65 s andt= 3.86 s along the LCA-LAD axial line obtained from the model with no particle (solid line) and the model with particles (dotted line) obtained from the domain with different intensity of stenosis at the proximal LAD: (a) 25%; (b) 50%; (c) 75%.
(a)
25% (b)
50% (c)
75%Figure 14
Wall shear stress at the peak systolet= 0.65 s andt= 3.86 s along the LCA-LCX axial line obtained from the model with no particle (solid line) and the model with particles (dotted line) obtained from the domain with different intensity of stenosis at the proximal LAD: (a) 25%; (b) 50%; (c) 75%.
(a)
25% (b)
50% (c)
75%Figure 15
The system of coronary arteries including the base of the aorta and the normal left and the normal right arteries.Figure 16
The system of coronary arteries including the base of the aorta and the 75% stenosed left artery and the normal right artery.Wall shear stress along the first half of the LCA connecting to the LAD is higher as shown in Figure13. At the first cardiac cyclet= 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 cyclet= 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).


## 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 SSTk-ω 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.Figure 2
Computational domain of validation study.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)t = 0.65 s; (b) t = 3.86s.
(a)
t = 0.65   s (b)
t = 3.86   sFigure 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.65s; (b)t= 3.86s.
(a)
t = 0.65   s (b)
t = 3.86   s


## 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 Figure5. 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.Figure 5
Computational domain of the LCA and its branches including the LAD and the 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 Figure6(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.Table 1
Mesh information.
% Stenosis Number of elements 0 % 314,150 25 % 318,897 50 % 350,265 75 % 444,764Figure 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)
LAD with stenosis (b)
The LCA-LAD axial line (c)
The LCA-LCX axial line (d)
Pulsatile pressure at the LCA inflow surfaceThe 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.Table 2
Model parameters.
Parameters value units Blood density (ρ) 1050 kg m−3 Zero shear viscosity (μ0) 0.56 g cm−1 s−1 Infinite shear viscosity (μ∞) 0.0345 g cm−1 s−1 Time constant (λ) 3.3130 s Consistency index (n) 0.3568 - Shape parameter (a) 2 - Particle diameter (dp) 100 μm Particle density (ρp) 1050 kg m−3 Particle mass flow rate (m˙p) 5.50 x 10−8 kg s−1 Particle time step (Δt) 1×10-6 sTable 3
Numerical setting.
Numerical Parameters Setting Solver Pressure - based Pressure – Velocity Coupling Coupled Multigrid AMG Under – relaxation factors 0.75 - momentum 0.75 - pressure 0.80 - turbulent kinetic energy 0.80 - specific dissipation rate 0.90 - discrete phase source Spatial discretization Second order upwind - momentum Second order upwind - turbulent kinetic energy Second order upwind - specific dissipation rate Second order upwind - transient formulation Gradient Green - Gauss node - basedFocusing on the proximal part of LCA as shown in Figure7(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.Figure 7
Dispersed phase flow in the LCA connecting to the LAD and LCX.
(a)
A part of LCA (b)
Streamline of blood flow (c)
Velocity field of blood flow (d)
Particle deposition patternsIn the fluid flow model, small turbulent flow appears in the transition area connecting the LCA with its branches as shown in Figures7(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.Figure 8
Vector plot of blood velocity at the peak systole at the timet= 0.65 s in three different planes.
(a)
Cutting plane (b)
Blood velocity in each planeDistribution of systolic pressure along two axial lines including the LCA-LAD line in as shown in Figure6(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 att= 0.65 s andt= 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 systolet= 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. Att= 0.65 s, a pressure drop is about 55 mmHg and 40 mmHg att= 3.86 s, respectively.The level of particle motion varies with time. The highest level is at the peak systole of the first cardiac cyclet= 0.65 s. To show the effect of particle motion on the blood pressure and the wall shear stress, the results at the peak systolet= 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 andt= 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).Figure 9
Systolic pressure att= 0.65 s in the artery model with 75% intensity of stenosis at the proximal LAD obtained from two different models: (a) model with no particle; (b) model with particle motion.
(a)
Model with no particle (b)
Model with particle motionFigure 10
Pressure distribution and wall shear stress distribution with particle deposition patterns in the artery model with 75% intensity of stenosis at the proximal LAD at the peak systole in the fifth cardiac cycle att= 3.86 s.
(a)
Pressure distribution with particle deposition patterns (b)
Wall shear stress distribution with particle deposition patternsFigure 11
Blood pressure at the peak systolet= 0.65 s andt= 3.86 s along the LCA-LAD axial line obtained from the model with no particle (solid line) and the model with particles (dotted line) obtained from the domain with different intensity of stenosis at the proximal LAD: (a) 25%; (b) 50%; (c) 75%.
(a)
25% (b)
50% (c)
75%Figure 12
Blood pressure at the peak systolet= 0.65 s andt= 3.86 s along the LCA-LCX axial line obtained from the model with no particle (solid line) and the model with particles (dotted line) obtained from the domain with different intensity of stenosis at the proximal LAD: (a) 25%; (b) 50%; (c) 75%.
(a)
25% (b)
50% (c)
75%Figure 13
Wall shear stress at the peak systolet= 0.65 s andt= 3.86 s along the LCA-LAD axial line obtained from the model with no particle (solid line) and the model with particles (dotted line) obtained from the domain with different intensity of stenosis at the proximal LAD: (a) 25%; (b) 50%; (c) 75%.
(a)
25% (b)
50% (c)
75%Figure 14
Wall shear stress at the peak systolet= 0.65 s andt= 3.86 s along the LCA-LCX axial line obtained from the model with no particle (solid line) and the model with particles (dotted line) obtained from the domain with different intensity of stenosis at the proximal LAD: (a) 25%; (b) 50%; (c) 75%.
(a)
25% (b)
50% (c)
75%Figure 15
The system of coronary arteries including the base of the aorta and the normal left and the normal right arteries.Figure 16
The system of coronary arteries including the base of the aorta and the 75% stenosed left artery and the normal right artery.Wall shear stress along the first half of the LCA connecting to the LAD is higher as shown in Figure13. At the first cardiac cyclet= 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 cyclet= 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 SSTk–ω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.



---

*Source: 2593425-2018-08-01.xml*