Computational evaluation of smoothed particle hydrodynamics for implementing blood flow modelling through CT reconstructed arteries

Abstract

Simulation of blood flow in a stenosed artery using Smoothed Particle Hydrodynamics (SPH) is a new research field, which is a particle-based method and different from the traditional continuum modelling technique such as Computational Fluid Dynamics (CFD). Both techniques harness parallel computing to process hemodynamics of cardiovascular structures. The objective of this study is to develop and test a new robust method for comparison of arterial flow velocity contours by SPH with the well-established CFD technique, and the implementation of SPH in computed tomography (CT) reconstructed arteries. The new method was developed based on three-dimensional (3D) straight and curved arterial models of millimeter range with a 25% stenosis in the middle section. In this study, we employed 1,000 to 13,000 particles to study how the number of particles influences SPH versus CFD deviation for blood-flow velocity distribution. Because further increasing the particle density has a diminishing effect on this deviation, we have determined a critical particle density of 1.45 particles/mm² based on Reynolds number (Re = 200) at the inlet for an arterial flow simulation. Using this critical value of particle density can avoid unnecessarily big computational expenses that have no further effect on simulation accuracy. We have particularly shown that the SPH method has a big potential to be used in the virtual surgery system, such as to simulate the interaction between blood flow and the CT reconstructed vessels, especially those with stenosis or plaque when encountering vasculopathy, and for employing the simulation results output in clinical surgical procedures.

Keywords

Blood flow simulation smoothed particle hydrodynamics stenosed artery critical density computed tomography parallel computing

1 Introduction

Stenosis in arteries is the result of atherosclerosis, and is one of the most common cardiovascular problem that can lead to the malfunction of the vascular system [1, 2]. The progression of plaque can result in the redistribution of wall shear stress inside the blood vessel wall, resulting in plaque growth and rupture occurrence. Before the rupture occurs, the stenosis can cause flow separation and oscillation, which lead to higher pressure drop and flow resistance [3, 4]. The accumulation of cholesterol plaque in the coronary artery can occlude the artery to cause myocardial ischemia. Therefore, it is significant to determine the relationship between arterial hemodynamics and plaques formation by using advanced simulation techniques for the study of this disease.

There are several studies on simulation methods for determination of cardiovascular flow. In the 1800s, Claude Navier and George Stokes, formulated the famous Navier-Stokes equations that describe the dynamics of fluids [6]. Historically, Computational Fluid Dynamics (CFD) has been implemented for fluid analysis, and became even more widely used since the rapid development of computers in the 1950s. In addition to the Navier-Stokes momentum conservation equations, two additional equations are needed to simulate fluid flows: (1) continuity equation describing mass conservation, and (2) state equation describing energy conservation. In the CFD literature, many methods have been proposed to simulate fluid flows [7]. Most CFD software, such as ANSYS, was developed based on fixed grid (2D) or volume (3D) in space to compute fluid flow velocity and pressure. We have studied fixed grids in space, and shown how the fluid quantities measured at these grids change in time. Using this method, we can determine the velocity data, pressure data and Wall Shear Stress(WSS) data.

In recent years, there have been significantly more researches on arterial hemodynamics. Many researchers have studied pathological mechanisms by analyzing arterial hemodynamics [5, 6]. Annich et al. researched on catechin treatment for improving cerebrovascular flow-mediated dilation and learning abilities due to atherosclerosis [5]. Callaghan et al. studied the stress distribution in cervical internal carotid artery, to investigate the role of the carotid sinus in the reduction of arterial wall stresses due to head movements [6]. Several studies have investigated the correlation of cerebrovascular atherosclerosis with Alzheimer’s disease [7 –9]. Alex et al. measured cerebral blood flow to research on the potential correlation between cerebral blood flow volumes, pulse pressure, and cognitive measurements [7]. Mark et al. investigated the frequency and severity of atherosclerotic plaques in the circle of Willis in Alzheimer’s disease [8]. Rovshan introduced the Dean number, which indicates the influence of curvature on the resistance to blood flow, to explain the specific correlation of Wills atherosclerosis correlates with Alzheimer’s disease [9].

Various other studies have carried out detail investigation of the hemodynamic parameters of blood flow in various arterial models [10 –18]. For examples, Bhaganagar et al. carried out numerical investigation of the effect of plaque morphology on the flow characteristics in a diseased coronary artery, by using realistic plaque morphology [11]. Chaichana et al. investigated the hemodynamic changes in various types of coronary stenosis in a patient-specific left coronary artery bifurcation [12]. Liu et al. conducted a numerical simulation of flow in curved coronary arteries with progressive amounts of stenosis, by using fluid-structure interaction modelling, to examine the influence of the variability in curvilinear arterial wall geometry on the blood flow characteristics [13]. Torii et al. carried out fluid-structure interaction analysis, to study physiological velocity and pressure waveforms in a patient-specific right coronary artery [14]. Wong et al carried out mathematical modelling of blood flow through an artery with multiple stenosis and post-stenotic dilatations, to examine the effects of variability in arterial wall geometry on the blood flow resistance [15]. Then, Wu et al. studied transient blood flow in elastic coronary arteries with varying degrees of stenosis and dilatations, in order to determine the haemodynamic parameters associated with stenosis and dilatations [16]. Furthermore, Siauw et al. employed both Newtonian and non-Newtonian models to simulate unsteady flow through a hypothetical stenotic geometry, in order to compare the difference of hemodynamic parameters between Newtonian and non-Newtonian models [17].

However, continuum models, such as those mentioned above, may be excessively time-consuming. A large number of methods have been proposed in CFD literatures to simulate virtual prototypes on computer. However, this has involved a number of tedious and time-consuming tasks [19, 20]. These methods include geometry translation, geometry cleanup, manual meshing processes and re-applying e.g. boundary-and volume-conditions whenever the product geometry changes. Therefore, time-consumption limits the application of CFD to real-time system, although it has ability to analysis simulation results. As computer technology has advanced rapidly over the past decades, discrete particle based methods such as Smoothed Particle Hydrodynamics (SPH) have experienced increasing utilization in arterial flow simulation [21]. In fact, meshless methods have arisen and grown in popularity, as they can be applied to problems that are highly nonlinear in arbitrarily complex geometries [22] and difficult for grid-based methods. As the meshless techniques become more researched, SPH is proving to be popular among computational hemodynamics researchers. More detailed reviews on standard SPH theory can be found in literature [23, 24].

The SPH theory possesses several benefits over traditional grid-based techniques. Firstly, SPH guarantees conservation of mass without extra computation, e.g. computing “divergence-free” equation [25, 26], since the particles themselves represent mass. Secondly, SPH computes pressure from weighted contributions of neighboring particles rather than by solving linear systems of equations. Moreover, for interactive applications [26, 27] such as virtual surgery, particle based fluid simulation methods like SPH are commonly preferred to grid-based fluid representations. This is because the fluid is able to flow freely in the scene and only perform computation where necessary, without the need to define a finite grid that is costly in terms of memory and computation. Particle methods are also more convenient to integrate into existing physics systems [28], as particles can collide against the scene geometry just like other rigid objects, without the need to voxelize the scene geometry [29, 30]. SPH also has advantages for modelling multiphase flow over grid-based method [31].

However, the drawback of SPH technology is that these simulations are mostly implemented for visualisation and cannot quantify the hemodynamics in System International (SI) units for measurement reference. In other words, they can only be used for aesthetic presentation of the cardiovascular flow, and lack the important quantifiable simulated units for analysis. This has motivated us to perform the research on calibrating SPH in reference to the CFD results that can be assumed to be ground truth data.

In this paper, we describe the SPH method and velocity data acquisition of our experiment in Section 2. Then, the vessel wall geometry, the parameters used in SPH and CFD and the framework of our comparative experiments are described in Section 3. The resulting velocity distribution contours based on SPH method are compared with the CFD results in Section 4, and we discuss therein the reasons for the deviation and associated complications. Finally, we conclude our study and propose future work in Section 5.

2 Theory of smoothed particle hydrodynamics

2.1 Theoretical formulations

We note again that the SPH is a particle-based method, namely the fluid is represented by particles and every particle has its quantities, such as mass, density and velocity, etc. Furthermore, SPH is an interpolation method. Field quantities defined at discrete particle locations can be evaluated anywhere in space. In order to do this, the SPH method distributes quantities in the local neighborhood of each particle, using symmetrical smoothing kernels. Therefore, a scalar quantity A is interpolated at location x_i by a weighted sum of contributions from all neighborhood particles, as shown in Fig. 1. Therein, the fluid is represented by spheres and the smoothing length h determines the resolution of simulations and controls the number of neighboring particles (the green solid spheres in circle) which are contributing to the evaluation of the quantity at any given point. The corresponding particle function is depicted as: $A_{S} (x_{i}) = \sum_{j} m_{j} \frac{A_{j}}{ρ_{j}} W (| x_{i} - x_{j} |, h),$ (1) where j iterates over all neighbor particles in the scope of h, m_j is the mass of particle h, x_j and ρ_j and A_j are the position density and arbitrary attributes of the field quantity at x_j. The ‘s’ indicates a smoothed quantity.

The scalar function W (x, h) is called the smoothing kernel with radius h. Since we only use kernels with finite support, so the ‘h’ is the radius of support in our formulation. Furthermore, we require W to be even such that W (x, h) = W (- x, h), and normalized smallintW (x) dx = 1. Another useful distributed particle function is the poly6 kernel [32]:

$W_{poly 6} (r, h) = \frac{315}{64 π h^{9}} {\begin{matrix} {(h^{2} - r^{2})}^{3} & 0 \leq r \leq h \\ 0 & otherwise \end{matrix},$ (2)

We now have all the ingredients to compute a smooth density field from the individual positions and messes of the neighborhood particles, as follows: $ρ (x_{i}) = \sum_{j} m_{j} \frac{ρ_{j}}{ρ_{j}} W (| x_{i} - x_{j} |, h) = \sum_{j} m_{j} W_{poly 6} (| x_{i} - x_{j} |, h),$ (3)

In the fluid flow equations, the derivatives of field quantities need to be computed. A good property of this formulation is how the gradient of such field quantities affects the smoothing kernel. In this regard, the gradient and the Laplacian of A are calculated by using the following formulae: $\nabla A_{S} (x) = \sum_{j} m_{j} \frac{A_{i}}{ρ_{j}} \nabla W (| x_{i} - x_{j} |, h),$ (4) $\nabla^{2} A_{S} (x_{i}) = \sum_{j} m_{j} \frac{A_{i}}{ρ_{j}} \nabla^{2} W (| x_{i} - x_{j} |, h) .$ (5)

2.2 Modelling fluid flow with SPH

In the Eulerian (grid based) formulation, the fluids are described by a series of fields, such as the velocity field v , density field ρ and pressure field p . The variation of these quantities over time is given by two equations. The first equation assures conservation of mass: $\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ v) = 0,$ (6) while the Navier-Stokes (N-S) equation assures conservation of momentum: $ρ (\frac{\partial v}{\partial t} + v \cdot \nabla v) = - \nabla p + ρ F + μ \nabla^{2} v,$ (7) where, F is an external force and μ is the viscosity of the fluid.

Differing from the Eulerian method, the particle system (instead of stationary grid) has mainly two advantages [26]. On the one hand, because the number of particles is constant and each particle has a constant mass, the mass conservation property is guaranteed; hence, Equation (6) can be completely omitted. On the other hand, the terms ∂v/∂t + v · ∇ v on the left hand side of Equation (7) can be replaced by the substantial derivative Dv/Dt [26]. The substantial derivative of the velocity field is simply the time derivative of the velocity of the particles, since the particles move with the fluid, which means that the term v · ∇ v is not needed for particle systems.

Based on this assessment, there are three forces on the right hand of the N-S equation, for modeling the influence of pressure (- ∇ p), external forces (ρF) and viscosity (μ∇²v). Equation (7) can now be rewritten as: $\frac{\partial v}{\partial t} = - \frac{1}{ρ} \nabla p + F + \frac{1}{ρ} μ \nabla^{2} v,$ (8) and

$a_{i} = \frac{\partial v_{i}}{\partial t}, v_{i} = \frac{\partial x_{i}}{\partial t},$ (9) where a_i is the acceleration of particle i.

By combining Equations (8) and (9), we can evaluate the velocities of particles if we can evaluate the forces, represented by the three terms on the RHS of Equation (8). Once the velocities are calculated, the displacements can be computed by integrating these velocities over time.

Particularly, the term F in Equation (8) represents types of forces. Firstly, for low Reynolds number flows, local variations in the pressure gradient forcing the fluid motion can be very small, compared with the hydrostatic pressure gradient. Another part of F is the body force derived from the hydrostatic pressure gradient, which drives flow [33]. The other parts of F are the forces exerted by the vessel wall on the particles contacting with the vessel wall: (i) constraint force to prevent fluid particles from moving across the vessel, and (ii) friction force applied to the particles.

With regards to the specific expressions for forces due to pressure gradient and viscosity, and the kernel functions, we mainly refer to Matthias Müller [26], for defining them as follows: $F_{i}^{Pressure} = - \frac{1}{ρ_{i}} \nabla p (x_{i}) = \frac{1}{ρ_{i}} \sum_{j} m_{j} \frac{p_{i} + p_{j}}{2 ρ_{j}} \nabla W_{spiky} (| x_{i} - x_{j} |, h),$ (10) $F_{i}^{Viscosity} = \frac{1}{ρ_{i}} μ \nabla^{2} v (x_{i}) = \frac{1}{ρ_{i}} μ \sum_{j} m_{j} \frac{v_{j} - v_{i}}{2 ρ_{j}} \nabla^{2} W_{viscosity} (| x_{i} - x_{j} |, h),$ (11) $p_{z} = k (\frac{ρ_{z}}{ρ_{0}} - 1), z represents i and j .$ (12) $W_{spiky} (r, h) = \frac{15}{π h^{6}} {\begin{matrix} {(h - r)}^{3} & 0 \leq r \leq h \\ 0 & otherwise \end{matrix},$ (13) $W_{viscosity} (r, h) = \frac{15}{2 π h^{3}} {\begin{matrix} - \frac{r^{3}}{2 h^{3}} + \frac{r^{2}}{h^{2}} + \frac{h}{2 r} - 1 & 0 \leq r \leq h \\ 0 & otherwise \end{matrix},$ (14)

The fluid molecules are subjected to attractive and repulsion forces, which are equal in all directions and balance each other. In contrast, the forces acting on the molecules at the free surface are unbalanced, which results in surface tension. For realistic simulation of blood flow, we use an alternative cohesion force [34], which is defined as: $F_{i}^{Cohesion} = - γ m_{i} m_{j} C (| x_{i} - x_{j} |) \frac{x_{i} - x_{j}}{| x_{i} - x_{j} |},$ (15) where γ is the surface tension coefficient, and C is a spline function created for 3D SPH simulation, which can be given as: $C (r, h) = \frac{32}{π h^{9}} {\begin{matrix} {(h - r)}^{3} r^{3} & \frac{h}{2} < r \leq h \\ 2 {(h - r)}^{3} r^{3} - \frac{h^{6}}{64} & 0 < r \leq \frac{h}{2} \\ 0 & otherwise \end{matrix},$ (16)

The cohesion term is similar to the Lennard-Jones potential interaction between a pair of molecules. However, this cohesion term stops increasing as the particles being closer, which avoids very stiff forces and resultant stability issues.

2.3 No-slip boundary conditions

The solid wall that the fluid is in contact with is considered to be impermeable. No fluid element is allowed to cross the wall boundary. Thus, the normal component of velocity has to be zero, namely: $\vec{u} \cdot \hat{n} = 0,$ (17) which pertains to the static vessel wall. However, for the elastic vessel, the normal component of the fluid velocity needs to match the normal component of the elastic solid’s velocity, as: $\vec{u} \cdot \hat{n} = {\vec{u}}_{solid} \cdot \hat{n},$ (18)

In both these equations, $\hat{n}$ is the vector normal to the solid boundary. This is the no-slip condition, since we restrict the normal component of fluid velocity, and only keep the tangential component of fluid velocity, for the fluid to slip past freely along the tangential direction of the solid wall boundary. The flow chart for the implementation of SPH method is shown in Fig. 2.

3 Experimental methods

3.1 Model reconstruction

In order to validate the measured results of the SPH method, we selected the 3D vessels computational model for modeling the flow in the regions of stenosis channel (as shown in Figs. 3 and 4, wherein both the outline size and mesh model are shown), and computed the flow velocity distribution with the CFD method using Ansys software.

The vessel models reconstructed from the CT/MRI data are typically represented by triangulated surface elements. The geometrical model of a 3D artery with one stenosis is shown in Fig. 3(a). The non-diseased sections have a constant diameter of D = 4.4 mm and the d = 2.4 mm. The axial geometry of an atherosclerotic lesion artery is characterized by the non-occluded wall height (D-d) and diameter D. The cylindrical domain has a total length of L = 50 mm, wherein L2 = 5 mm from the inlet segment. The diseased stenotic segment has a length of L3 = 7.5 mm, and the stenotic segment is at a distance L4 = 5 mm from the outlet segment. The vessel’s inlet segment is L1 = 10 mm and the outlet segment is L5 = 22.5 mm; these inlet and outlet segments L1 and L5 are meant to eliminate the inlet and outlet effects. There are 1308 surface mesh elements as depicted in Fig. 4(a), and 28,782 volume elements as depicted in Fig. 4(b).

Figure 3(b) presents a curved arterial segment in which α = 60°, R = 11.456, the diameter of vessel lumen D = 4.4 mm, the thickness of vessel wall is 0.8 mm, the straight length of inlet L1 = 10 mm, the outside length of curved segments L2 = 7.5 mm and L4 = 7.5 mm, the straight length of outlet L5 = 20 mm. Also, d is the diameter of vessel lumen in the narrowest plane and equal to 2.4 mm, D-d is the height of plaque, L3 = 15 mm is the length of the plaque. The model is represented by 3290 surface mesh elements as depicted in Fig. 4(c), and 347,000 volume elements shown as Fig. 4(d). The grid independence test was carried out to determine the optimal number of theseelements.

3.2 Parameter settings for boundary conditions

In the simulation, the parameters used in SPH and CFD are shown in Table 1, which are set according to human body blood parameters – dynamic viscosity, density, velocity at inlet [16, 35], which is corresponding to the Reynolds number Re = 200 considered to be a laminar flow [36]. In the present study, we focus on steady blood flow in vessels, namely that the velocity at the inlet is constant and the outlet pressure is set to be at 0 Pa. The flow models are processed using GPU-based parallel computing framework.

3.3 Data acquisition in the SPH method

In order to compare velocity contours of blood flow with the CFD method, we need to acquire the velocity data of particles. Considering that the fluid is granular and less uniform in the SPH method simulation, implying that the particles may not fill the vascular cavity entirely, we take measures to weaken that effect and increase precision when acquiring the velocity data. On the one hand, after getting the positions and velocities of the particles in a plane, we put the particles into a uniform grid according to their positions (resolution is 4 grid/mm²); we then calculate the average velocity value for every grid with the particles in it. On the other hand, we acquire the velocity data for each grid 10–15 times and use the average as its velocity. The process of data acquisition is shown as Fig. 5.

3.4 Validation of SPH by using CFD

CFD methods are generally regarded to be effective ways for obtaining solutions to flow problems. Analysis based on CFD that utilizes the anatomical reconstructed arterial geometry forms a new branch of study termed as computational hemodynamics. The arterial flow field has been extensively studied by CFD simulation over many years [37, 38]. Hence, we can assume the computed results employing CFD as the properties of true flow, and evaluate SPH by comparing against it.

It is important to establish a measure for percentage deviation of the examined velocity measurement matrix V_E. The deviation matrix θ can be defined as the ratio of the deviation map V_D to the velocity of true flow V_T, as follows: $θ = \frac{V_{D}}{V_{T}} \times 100 %,$ (19) $V_{D} = | V_{E} - V_{T} |,$ (20) $\bar{θ} = \frac{1}{n^{2}} \sum_{i = 1}^{n} \sum_{j = 1}^{n} θ_{ij},$ (21) where n is the width of deviation map in pixels. Ideally, if the system can generate flow velocity fields that correspond close enough to those that are measured by true flow, the deviation tends to approximate zero, and θ also tends to zero.

Based on velocity contour maps of straight and curved arteries, (i) the sum of percentage deviations θ, for all pixels in the contour map, is computed, and (ii) the average of all deviations represented by every pixel $\bar{θ}$ is then determined, and used as our indicator of how different the SPH predicted flow is from its ground truth.

The framework of our study is depicted by Fig. 6. After the geometry of the vessel structure (usually triangular surface mesh) is obtained, for CFD a grid/voxelization feature is adopted for the model. Next, a suitable algorithm solves the Navier-Stokes equations of motion, to get the value of every grid/voxelization which represents fluid field. For SPH, there is no need to voxelize the model, since the particles represent fluid itself and the attributed variables (of velocity, position, etc.) with them. We have thereby solved the resolution stated in section 2.2, and depicted the entire simulation process of blood flow through the vessels.

4 Results and discussion

4.1 Flow of particles through stenosed straight and curved arteries by SPH

The flow of particles through stenosed arteries is created by using our in-house SPH based system (developed using Microsoft Visual Studio 2010) that has GPU-accelerated interactive SPH particle simulation and rendering capability. Based on our SPH framework, we can simulate more than 10,000 particles at 3.4 G Hz on a GPU with a GeForce GTS 450. We can achieve a final simulation and average render frame rate of more than 15 FPS. We have developed time-dependent flow photos for both the straight and curved arteries, as shown in Fig. 7 (a) and (b) respectively.

4.2 Comparison of Velocities based on SPH and CFD simulation methods

The comparison of flow ve20-02-2017locity contours inside the stenosed vessels by SPH and CFD is illustrated in Figs. 8 & 9 for the straight vessel with one stenosis, and Figs. 10 & 11 for the vessel with segment bend and stenosis. We present the velocity contours of coronal and sagittal sections. Their coordinate systems are the same as in Figs. 4 and 5, and the A, B, C sagittal planes are shown in Fig. 4. The resolution of our contour maps is 4 grids/mm².

As the results in Figs. 8–11 demonstrate, the SPH method produces a close agreement with the CFD method for simulation of the velocity distribution in straight and bent stenoic vessel models. Figures 8 and 11 show the different plane velocity contours in the coronal and sagittal planes of flow through (i) a straight vessel (cylinder cavity) with a stenosis, and (ii) a bent vessel with a stenosis.

In Fig. 8, we can see that the velocity distribution with both the methods corresponds to laminar flow, wherein the location of the site at which the velocity maximizes is the centreline of stenosis in the axial direction, and the minimum velocity is acquired near the wall. In the lateral direction, the maximum velocity is closely approximated at the minimum lumen of stenosis throat and also in an area downstream of the stenosis.

Figure 10 shows the velocity contours of flow through a bent vessel with a stenosis. Therein, the velocity distributions are observed to be highly skewed toward the outer vessel walls; such a skewed flow structure is primarily caused by the misalignment of the mean axis of arteries [39]. While the blood flow in the straight vessel is laminar, the symmetry is broken in the case of the bent vessel with stenosis. Especially in its outlet segment, the velocity distributions are observed to be highly skewed toward the outer vessel walls; such a skewed flow structure is primarily caused by the misalignment of the mean axis of arteries [39].

4.3 Accuracy of the SPH versus CFD

As we can observe from above depicted results, the SPH velocity contours have some deviation from those obtained by the CFD method. There are some main reasons for the deviations. Firstly, the CFD is a grid-based method, which calculates the properties of the simulation at a fixed set of points in space, which means that every grid in space has property values. On the other hand, the SPH is a particle-based method, and calculates the properties of a set of particles as they move through space. In the other words, the entire space may be filled with fewer particles, especially when the geometric construction changes suddenly past the stenotic section. Hence, we could increase the number of particles to decrease the deviation.

Figure 12 shows the performance of our system versus the number of particles used. It is to note that we implemented 10 trials, and adopted the average value of the trials as the examined result. The percentage deviation is inversely proportional to the number of particles. As such, it is worthwhile to use more particles to improve the accuracy of flow. However, there is a weak deviation effect when particle number exceeds 10,000, which corresponds to the corresponding particle density of 1,454,100 particles per square meter in the stenosis section. Therefore, we determine this to be the critical density in our arterial flow simulation.

In our simulation flow study, we have also observed that the lowest $\bar{θ}$ occurs in the narrowest stenosis section of straight models, while the poorer accuracies occur in outlet sections of straight and curved models as depicted in Fig. 13. Because the particles become more tight and dense when passing through the narrowest stenosis section, it means that there are more particles per unit area, which then leads to a distribution closer to true flow by the CFD method. However, the particles become looser and more scattered after flowing into the outlet segment.

The lowest percentage deviation occurs in the stenosis section, due to the highest particle density, namely the density based on particle counts per cross-sectional area of flow is inversely proportional to the area through which the particles flow. However, the outlet section has a higher deviation than the inlet section with the same area; this implies that the particle density distribution may be uneven, as can be observed from Figs. 9(c) and 11(c). In our study, we have observed that while particles are initially arranged in order, they disorder and loosen close to the outlet section, which introduces random force deviations [40].

Secondly, as we described in section 2.4, although we grid the plane to acquire data, the distribution of particles in the vessel models is also non-uniform and particularly close to the border (as illustrated by Figs. 8–11 (a). We have used the force exerted by the vessel wall to prevent fluid particles from moving across the vessel, which has introduced an artificial term to correct the movement of the particle. This correction is non-physical and results in a rough edge, as shown in Figs. 8–11 (a) This causes a large proportion of difference (deviation), as shown by the solid black circles in Figs. 8–11 (c). The conventional low order SPH cannot produce high accuracy results, especially in the area near the boundary of the fluid domain [41]. In this paper, we have focused on the particle number that has an effect on deviation, to obtain the particle critical density that can avoid unnecessary huge computational expenses and further effect the simulation accuracy. For our SPH method, it is possible for us to increase the number of particles, in order to produce more accurate flow information, and we have in fact obtained a critical value of the particle number. Further, it is more important and effective to study the factors influencing its accuracy, as carried out by some researchers [41, 42].

Although the simulation time increases with the increasing number of particles, which is depicted in Fig. 14, it is still much less than CFD method’s costs. In our experiments, the computational time of the simplest straight and bent vessel using CFD are 6.7 minutes and 8.4 minutes respectively. Considering many complex and irregular artery models that may exist, it is extremely difficult to simulate and compute fluid flowing in real-time due to the expensive time cost with CFD simulation. In general, the velocity contours obtained by SPH are consistent with those obtained by CFD. Therefore, we can state that while the CFD simulation method is more suited for the pure analysis of data, the SPH method is better suited for real time simulation like for virtual surgery simulation to provide the information for surgical implementation. We have obtained tolerable results balancing between time and accuracy.

The most significant advantage of SPH is that we can show the simulation process of blood flow in a 3D meshed vessel in real-time. This is particularly important in the use of SPH method in presurgical simulation, as depicted by us in Fig. 15, which provides a real-time presentation of blood flow through a realistic artery model employed in our virtual surgery system.

5 Conclusion

To the best of our knowledge, this paper is the first to compare the results of SPH simulation with CFD simulation in stenotic vessel models. It has been verified by us that the SPH method enables a fairly accurate determination of velocity distribution in stenotic arterial models; a more thorough analysis and calibration is still to follow. The results of velocity distribution between our SPH and CFD simulations study reveal a few deviations, particularly at the outlet segment of the vessels. In this regard, we have successfully studied the criticality of particles number on simulation accuracy for achieving resolution independence.

The SPH is a mesh-free method, adopting a discrete set of particles without topological connectivity (i.e. a grid) for representing the material continuum, which is the most essential reason to cause deviation when validation by using CFD. Nevertheless, the SPH method has a significant advantage for arbitrarily complex geometries (such as an anatomical vessel model reconstructed from CT/MRI data, whose curvature and radius are varied and complicated.) and for interactive systems, which are difficult to simulate by grid-based CFD methods. We can also state that it is possible to usefully apply SPH to an interactive system like virtual surgery blood flow simulation, and get both the scenario and hemodynamic parameters results. Therefore, it is worthwhile to conduct further research on simulation of blood flow by using the SPH method.

As for future work, Reynolds number may have effect on the deviation, which is also worth studying.. Furthermore, we plan to improve the accuracy of interpolation methods and the kernel function in order to solve the situation of particle disorder arrangements in SPH, so that we can provide more accurate specific information about blood flow simulation in our virtual interactive surgery system or other clinical remedial measures for arteries.

Conflict of interest

We confirm that all authors of this manuscript have no conflicts of interest to declare.

Footnotes

Acknowledgments

This work was supported in part by National Natural Science Foundation of China (No. 61672510), Shenzhen Science and Technology Program (No. JSGG20150602143414338, No. JCYJ20160331191401141), and Guangdong Science and Technology Program (No.2016A020220016).

References

Nerem

R.M.

, Vascular fluid mechanics, the arterial wall, and atherosclerosis, Journal of Biomechanical Engineering114(3) (1992), 274–282.

Worthley

S.G.

, Omar-Farouque

H.M.

, Helft

and Meredith

I.T.

, Coronary artery imaging in the new millennium, Heart Lung Circ11(1) (2002), 19–25.

Wong

K.K.L.

, Mazumdar

J.N.

, Pincombe

, Worthley

S.G.

, Sanders

and Abbott

, Theoretical modeling of micro-scale biological phenomena in human coronary arteries, Medical & Biological Engineering & Computing44(1) (2006), 971–982.

M.X.

, Beech-Brandt

J.J.

, John

L.R.

, Hoskins

P.R.

and Easson

W.J.

. Numerical analysis of pulsatile blood flow and vessel wall mechanics in different degrees of stenoses, Journal of Biomechanics40(16) (2007), 3715–3724.

Annick

, Virginie

, Nathalie

T.T.

, Élisabeth

, Priscilla

, Edward

, Frédéric

, Marc-Antoine

, Louis

, Edith

, Guylaine

and Eric

, Catechin treatment improves cerebrovascular flow-mediated dilation and learning abilities in atherosclerotic mice, Am J Physiol Heart Circ Physiol300(3) (2011), H1032–H1043.

Callaghan

F.M.

, Soellinger

, Baumgartner

R.W.

, Poulikakos

, Boesiger

and Kurtcuoglu

, The role of the carotid sinus in the reduction of arterial wall stresses due to head movements-potential implications for cervical artery dissection, Journal of Biomechanics42(6) (2009), 755–761.

Alex

E.R.

, Josef

P.D.

, Michael

, Kewei

, James

G.P.

, Sharmeen

, Christine

, Chera

L.M.

, Pradeep

, Hua

, Jesse

M.H.

, Tyler

A.K.

, Douglas

G.W.

, Jane

C.K.

, Marek

, Marwan

N.S.

and Thomas

G.B.

, Cerebral blood flow in Alzheimer’s disease, Vasc Health Risk Manag8 (2012), 599–611.

Mark

, Sharon

X.X.

, Mitchel

A.K.

, Jon

B.T.

and David

A.W.

, Cerebrovascular atherosclerosis correlates with Alzheimer pathology in neurodegenerative dementias, Brain135(12) (2012), 3749–3756.

Rovshan

M.I.

, Circle of Willis atherosclerosis, Alzheimer’s disease and the Dean number, World J Cardiol5(10) (2013), 394–396.

10.

Lee

, Shim

and Chang

H.J.

, Intensity-vesselness Gaussian mixture model (IVGMM) for 2D+t segmentation of coronary arteries for X-ray angiography image sequences, Journal of X-ray Science and Technology23(5) (2015), 579–592.

11.

Bhaganagar

, Veeramachaneni

and Moreno

, Significance of plaque morphology in modifying flow characteristics in a diseased coronary artery: Numerical simulation using plaque measurements from intravascular ultrasound imaging, Applied Mathematical Modelling37(7) (2013), 5381–5393.

12.

Chaichana

, Sun

and Jewkes

, Hemodynamic impacts of various types of stenosis in the left coronary artery bifurcation: A patient-specific analysis, Phys Med29(5) (2013), 447–452.

13.

Liu

, Wu

, Huang

, Wu

, Zhang

, Wong

K.K.L.

and Ghista

D.N.

, Numerical simulation of flow in curved coronary arteries with progressive amounts of stenosis using fluid-structure interaction modelling, Journal of Medical Imaging and Health Informatics4(4) (2014), 605–611.

14.

Torii

, Wood

N.B.

, Hadjiloizou

, Dowsey

A.W.

, Wright

A.R.

, Hughes

A.D.

, Davies

, Francis

D.P.

, Mayet

and Yang

G.Z.

, Fluid–structure interaction analysis of a patient-specific right coronary artery with physiological velocity and pressure waveforms, Communications in Numerical Methods in Engineering25(5) (2009), 565–580.

15.

Wong

K.K.L.

, Tu

, Mazumdar

J.N.

and Abbott

, Modelling of blood flow resistance for an atherosclerotic artery with multiple stenoses and poststenotic dilatations, ANZIAM Journal51 (2010), C66–C82.

16.

, Liu

, Huang

, Ghista

D.N.

and Wong

K.K.L.

, Transient blood flow in elastic coronary arteries with varying degrees of stenosis and dilatations: CFD modelling and parametric study, Computer Methods in Biomechanics and Biomedical Engineering (2014), 1–11.

17.

Siauw

W.L.

Unsteady stenosis flow prediction, a comparative study of non-Newtonian models with operator splitting scheme, Medical Engineering and Physics22(4) (2000), 265–277.

18.

N.H.

, Chen

T.B.

, Liu

K.Y.

, Hsu

S.Y.

, Twan

W.H.

, Ding

H.J.

, Hung

C.M.

, Lin

L.W.

and Huang

Y.H.

, Investigated geometrical characteristics and image density of left ventricle of multi-detector computed tomography in early coronary artery disease patients, Journal of X-ray Science and Technology24(3) (2016), 353–359.

19.

Liu

, Wen

, Chao

, Yin

, Shen

, Lai

, Lin

C.-H.

, Liu

, Sun

and Chen

, Accurate and high-resolution boundary conditions and flow fields in the first-class cabin of an MD-82 commercial airliner, Atmospheric Environment26 (2012), 33–44.

20.

Mazumdar

, Long

and Chen

, A coupled CFD and analytical model to simulate airborne contaminant transmission in cabins, Indoor and Built Environment23 (2014), 946–954.

21.

Qin

and Pang

W.M.

, Particle-based simulation of blood flow and vessel wall interactions in virtual surgery, Proceedings of Symposium on Information and Communication Technology (2010), 128–133.

22.

Adami

, Hu

X.Y.

and Adams

N.A.

, A generalized wall boundary condition for smoothed particle hydrodynamics, Journal of Computational Physics231(21) (2012), 7057–7075.

23.

Monaghan

J.J.

, Smoothed particle hydrodynamics, Rep Prog Phys68 (2005),, 1703–1759.

24.

Harada

, Koshizuka

and Kawaguchi

, Improvement in the boundary conditions of Smoothed Particle Hydrodynamics, Compute Graphics & Geometry9(3) (2007), 2–15.

25.

Jos

, Stable fluids, ACM SIGGRAPH (1999), 121–128.

26.

Müller

, Charypar

and Gross

, Particle-based fluid simulation for interactive applications, In: Proc. Of Sigraph Symposium on Computer Animation (2003), 154–159.

27.

Damien

, Fabrice

and Marie-Paule

, Interactive animation of ocean waves, In, Proceedings of the ACM SIGGRAPH symposium on Computer animation (2002), 161–166.

28.

Akinci

, Dippel

and Akinci

, Screen Space Foam Rendering, Journal of WSCG21(3) (2013), 173–182.

29.

Steinman

D.A.

, Image-based CFD modeling in realistic arterial geometries, Annals of Biomedical Engineering30(4) (2002), 483–497.

30.

Severin

, Beckert

and Herwig

, Spatial development of disturbances in plane Poiseuille flow: A direct numerical simulation using a commercial CFD code, International Journal of Heat and Mass Transfer44(22) (2001), 4359–4367.

31.

Tartakovsky

A.M.

, Ferris

K.F.

and Meakin

, Lagrangian particle model for multiphase flows, Computer Physics Communications180 (2009), 1874–1881.

32.

Müller

, Schirm

and Teschner

, Interactive blood simulation for virtual surgery based on Smoothed Particle Hydrodynamics, Technology and Health Care12(1) (2004), 25–31.

33.

Joseph

P.M.

, Patrick

O.X.

and Yi

, Modeling low Reynolds Number incompressible flows using SPH, Journal of Computational Physics136(1) (1997), 214–226.

34.

Akinci

, Akinci

and Teschner

, Versatile surface tension and adhesion for SPH fluids, ACM Transactions on Graphics32(6) (2013).

35.

Johnston

B.M.

, Johnston

P.R.

, Corney

and Kilpatrick

, Non-Newtonian blood flow in human right coronary arteries: Steady state simulations, Journal of Biomechanics37(5) (2004), 709–720.

36.

Avila

and Moxey

, The onset of turbulence in pipe flow, Science333 (2011), 192–196.

37.

Bertolotti

, Qin

and Lamontagne

, Influence of multiple stenoses on Echo-Doppler functional diagnosis of peripheral arterial disease: A numerical and experimental study, Annals of Biomedical Engineering34 (2004), 564–574.

38.

, Liu

and Huang

39.

Vetel

, garon

and Pelletier

, Lagrangian coherent structures in the human carotid artery bifurcation, Exp Fluid46(6) (2009), 1067–1079.

40.

Dehnen

and Aly

, Improving convergence in smoothed particle hydrodynamics simulations without pairing instability, Instrumentation and Methods for Astrophysics425(2) (2012), 1068–1082.

41.

Zheng

, Duan

W.Y.

and Ma

Q.W.

, Comparison of improved meshless interpolation schemes for SPH method and accuracy analysis, Marine Sci Appl9 (2010), 223–230.

42.

El-Gammal

, Khalil

E.E.

and Haridy

, Influence of smoothing length and virtual particles on SPH accuracy, International Journal of Materials, Mechanics and Manufacturing1 (2013), 166–170.