Comparison of CFD and FSI Simulations of Blood Flow in Stenotic Coronary Arteries

Abstract

Cardiovascular diseases are amongst the main causes of death worldwide, and the main underlying pathological process is atherosclerosis. Over the years, fatty materials are accumulated in the arterial which consequently hinders the blood flow. Due to the great mortality rate of this disease, hemodynamic studies within stenotic arteries have been of great clinical interest, and computational methods have played an important role. Commonly, computational fluid dynamics methods, where only the blood flow behavior is considered, however, the study of both blood and artery walls’ interaction is of foremost importance. In this regard, in the present study, both computational fluid dynamics and fluid-structure interaction modeling analysis were performed in order to evaluate if the arterial wall compliance affects considerably the hemodynamic results obtained in idealized stenotic coronary models. From the overall results, it was observed that the influence of wall compliance was noteworthy on wall shear stress distribution, but its effect on the time-averaged wall shear stress and on the oscillatory shear index was minor.

1. Introduction

Over the last year, COVID-19 has been the most widely spoken and researched disease worldwide and, inevitably, other existing pathologies were moved to the background. Amongst these, cardiovascular diseases (CVDs) should be highlighted because they still are by far the major contributor to global mortality [1]. CVDs are mainly caused by a blockage that prevents blood from flowing properly, known as atherosclerosis [2, 3]. This is a complex disease that affects medium and large size arteries and consists of a build-up of fatty deposits on the inner walls of the blood vessels, hampering the blood flow through the body [4]. The effect of atherosclerosis can be exacerbated by other diseases that affect blood circulation. Particularly, it has been shown that patients with COVID-19 are prone to develop blood clots on both arteries and veins [5], and thus atherosclerosis may be an even more important factor to global mortality.

Given the prevalence of this disease, atherosclerosis has been intensely studied through both cardiovascular modeling and experimental procedures, as reviewed elsewhere [6, 7, 8, 9]. Nevertheless, with the growing trend of greater computer power, computational approaches have become a valuable, cheaper, and efficient alternative for numerous researchers to predict blood flow behavior [10, 11, 12, 13].

There are two main approaches for simulating blood flow, computational fluid dynamics (CFD) and fluid-structure interaction (FSI). In the first one, the arterial wall is assumed as rigid, while in the second one, arteries are considered elastic and the interaction between the blood and the arterial walls are included in the simulation [14]. Although CFD has been widely applied in the study of blood flow under pathological conditions in virtue of lower computational cost [15, 16, 17, 18, 19, 20], since FSI provides a more realistic simulation of the human vasculature behavior, it has received increasing interest [21, 22]. Nonetheless, this approach requires significantly more computational effort, and the foremost difficulty is stability and convergence [23, 24].

In order to evaluate if the differences between CFD and FSI results are significant, some researchers have investigated and compared both. A commonly mentioned work was developed by Torii et al. [25]. The authors studied the effects of wall compliance on a stenotic patient-specific coronary artery and found noticeable differences in the instantaneous wall shear stress (WSS) produced by the FSI and rigid wall models. However, the effects of wall compliance on time-averaged WSS (TAWSS) and oscillatory shear index (OSI) were negligible. Malvè et al. [26] performed a similar study in the left coronary artery bifurcation and the conclusions regarding the WSS agreed with the previous study, but they observed significant differences in the TAWSS, especially on its spatial distribution. More recently, another similar study was performed on carotid bifurcation [10]. The authors found that the rigid model overestimates the flow velocity and WSS, but its influence on the TAWSS is minimal.

The differences between the CFD and FSI simulations have thus been the subject of several studies, however, as demonstrated, there is still a debate as to whether it is really necessary to use the most realistic approach, namely in computing WSS dependent variables, and sometimes the findings are contradictory. In this regard, this work presents the comparison of the results of both CFD and FSI simulations in an idealized stenotic coronary artery, with a degree of stenosis of 50%.

2. Methods

The first step necessary to study blood flow consisted in obtaining the three-dimensional geometry of the lumen of an idealized coronary artery as depicted in Figure 1. This geometry was previously adopted in other works by these authors [27, 28].

Figure 1.

Geometry and mesh of the coronary artery for both solid and fluid domain with 50% of stenosis.

In the construction of the solid domain, a thickness of 0.8 mm was considered, according to an in vivo study [29].

Moreover, the fluid domain and the solid domain were discretized in 428,800 and 15,480 hexahedral elements, respectively, making a total of 444,280 elements for the FSI simulations. To ensure the quality of the mesh, the skewness and orthogonal quality were evaluated. The average skewness and orthogonal quality were approximately 8.3 $e^{-2}$ and 0.98, respectively. These parameters prove that the mesh is reliable for the study presented [30] and it is worth mentioning that a previous mesh study was carried out for CFD simulations [16], which was then adapted for this study.

2.1 Mathematical formulation

Blood flow is governed by the incompressible Navier-Stokes and the continuity equations as described in Eqs. (1) and (2),

$$\nabla ·u=0$$E1

$$\rho \left(\frac{\mathrm{\partial u}}{\partial t}+u·\nabla \right)=-\nabla p+\mu {\nabla }^{2}u$$E2

where $u$ is the velocity, $p$ is the static pressure, $\rho$ is the fluid density, and $\mu$ the dynamic viscosity [30, 31].

In the present study, blood was modeled as incompressible, laminar, and non-Newtonian fluid, having a density of 1060 kg/${\mathrm{m}}^3$ [32]. Although it is commonly assumed as a Newtonian fluid, the ability of red blood cells to deform and aggregate makes blood a non-Newtonian fluid [33]. The well-known Carreau model was used to simulate the shear-thinning blood behavior, and it is defined by Eq. (3) [30, 31]:

$$\mu ={\mu }_{\infty }+\left({\mu }_0-{\mu }_{\infty }\right){\left[1+\mathit{\lambda \gamma }{̇}^2\right]}^{\frac{n-1}{2}}$$E3

where μ is the viscosity, ${\mu }_{\infty }=0.00345\mathrm{Pa}\bullet \mathrm{s}$ is the infinite shear viscosity, ${\mu }_0=0.0560\mathrm{Pa}·\mathrm{s}$ is the blood viscosity at zero shear rate, $\gamma ̇$ is the instantaneous shear rate, $\lambda =3.313\mathrm{s}$ is the time constant and $n=0.3568$ is the power-law index, as previously applied by [27].

The governing equation for the solid domain is the equilibrium equation (Eq. (4)) [30, 31]:

$${\rho }_s\frac{{\partial }^{2}u}{\partial t^2}-\nabla \stackrel{=}{\sigma }={\rho }_s\overrightarrow{b}$$E4

where ${\rho }_s$ is the solid density, $u$ represents the solid displacements, $\overrightarrow{b}$, the body forces applied on the structure, and $\stackrel{=}{\sigma }$ is the Cauchy stress tensor. For an isotropic linear elastic solid, the stress tensor is represented by Eq. (5) [30, 31]:

$$\stackrel{=}{\sigma }=2{\mu }_L\stackrel{=}{\epsilon }+{\lambda }_L\mathit{tr}\left(\stackrel{=}{\epsilon }\right)I$$E5

where ${\lambda }_L$ and ${\mu }_L$ are the first and second Lamé parameters, respectively, $\stackrel{=}{\epsilon }$, the strain tensor, $\mathit{tr}$, the trace function, and $I$, the identity matrix. For compressible materials, Lamé parameters can be written as a function of Young’s modulus, $E$, and Poisson’s coefficient, $v$, as follows.

$${\lambda }_L=\frac{\mathit{vE}}{\left(1+v\right)\left(2v-1\right)}$$E6

$${\mu }_L=\frac{E}{2\left(1+v\right)}$$E7

The arterial wall was modeled as a linear elastic, incompressible, isotropic, and homogeneous material with Young’s modulus of $3.77\mathrm{MPa}$ [34], a density of $1120\mathrm{kg}/{\mathrm{m}}^3$ [35], and a Poisson’s ratio of $0.49$[21].

The FSI simulations were performed using the Arbitrary-Lagrangian-Eulerian (ALE) methodology for the fluid flow. Taking into account that the interface between the lumen and the wall deforms, the equations governing fluid flow have to be expressed in terms of the fluid variables relative to the mesh movement. The ALE-modified Navier-Stokes momentum equation for a viscous incompressible flow is described as follows in Eq. (8) [30, 31]:

$${\rho }_f\left(\frac{\partial u}{\partial t}+\left(\left(u-u_g\right)·\nabla \right)u\right)=-\nabla p+\mu {\nabla }^{2}u$$E8

where ${\rho }_f$, $p$, $u$, and $u_g$ are the fluid density, the pressure, the fluid velocity, and the moving coordinate velocity, respectively. The term $\left(u-u_g\right)$, in the ALE formulation, is added to the conventional Navier-Stokes equations to account for the movement of the mesh.

The displacement and equilibrium forces at the interface are represented by Eqs.(9) and (10) [30, 31]:

$$u_{f,\Gamma }=u_{s,\Gamma }$$E9

$$\overrightarrow{t_{f,\Gamma }}=\overrightarrow{t_{s,\Gamma }}$$E10

where $u_{f,\Gamma }$ is the displacement of the fluid at the interface, $u_{s,\Gamma }$, the displacement of the solid at the interface, $\overrightarrow{t_{f,\Gamma }}$, the forces of the fluid on the interface and $\overrightarrow{t_{s,\Gamma }}$, the forces of the solid on the interface.

2.2 Boundary conditions

Regarding the boundary conditions used in this study, at the inlet, a physiologically accurate pulsatile velocity profile was set, which is depicted in Figure 2. At the outlet, a pressure of 80 mmHg was assumed [17, 36].

Figure 2.

Velocity profile implemented in CFD and FSI simulations.

The solid and fluid wall-boundaries were defined as a fluid-structure interface, and the inlet/outlet adjacent solid boundaries were fixed in all directions.

2.3 Numerical solution

For the CFD simulations, the Ansys Fluent software was used which applies the finite-volume discretization method. In this method, the fluid domain is divided into a finite number of control volumes, the conservation equations are applied to each control volume. Then, a system of algebraic equations for the variables is obtained. For the velocity-pressure coupling, the semi-implicit method for the pressure-linked equations (SIMPLE) scheme was used [37].

In FSI simulations, the same finite-volume method is applied in the fluid domain, and the computed forces in Fluent are transferred to the solid domain, through the interface. The finite element method (FEM) is used to solve the governing equations of the solid domain. Then, the computed displacements are transferred back to the fluid domain. This two-way coupling process was repeated until the difference of the displacements and forces for the last two iterations is below 1%. A time step of 0.01 s was used for every simulation.

2.4 Hemodynamic parameters

The formation of atherosclerotic lesions has been widely associated with hemodynamic parameters, such as the wall shear stress (WSS) and its indices, time-averaged wall shear stress (TAWSS), and oscillatory shear index (OSI) [38, 39]. These have been very useful to predict and estimate disturbed flow conditions and the development of local atherosclerotic plaques [27, 40].

The spatial WSS, ${\tau }_w$, is calculated by Eq. (11), being $\stackrel{.}{\mathrm{\gamma }}$, the deformation rate, and $\mu$ the dynamic viscosity.

$${\tau }_w=\mu \frac{\partial u}{\partial y}=\mu \stackrel{.}{\gamma }$$E11

The TAWSS index allows obtaining an average temporal evaluation of the WSS exerted during a cardiac cycle ($T$) [40]. This is calculated by Eq. (12):

$$\mathit{\text{TAWSS}}=\frac{1}{T}{\int }_0^T\left|\mathit{WSS}\right|\mathit{dt}$$E12

The OSI index is the temporal fluctuation of low and high average shear stress during a cardiac cycle ($T$) and it is calculated by applying Eq. (13):

$$\mathit{OSI}=\frac{1}{2}\left(1-\frac{\left|{\int }_0^T\mathit{\text{WSSdt}}\right|}{{\int }_0^T\left|\mathit{WSS}\right|\mathit{dt}}\right)$$E13

The formulation developed in this section describes the models that couple the fluid dynamics and the mechanical interaction with the arteries’ wall which is treated as a deformable material. This methodology enables the computation of critical parameters for understanding the hemodynamics in the presence of a stenosis, such as the WSS. The advantages of this method are made evident by comparing it with a simple CFD analysis as detailed in the following section.

3. Results and discussion

3.1 Velocity profiles along the cardiac cycle

Figure 3 illustrates the velocity profiles measured along the center of the artery for both CFD and FSI simulations, in two different phases of the cardiac cycle: systole (0.4 s) and diastole (0.58 s).

Figure 3.

Axial flow velocity profiles at the center line drawn across the artery at (a) systole and (b) diastole.

Looking at the results in Figure 3, it can be observed that the estimated velocity profiles are similar for both CFD and FSI modeling. Moreover, as expected, during diastole (Figure 3b) the velocities measured are higher than during systole (Figure 3a). This happens because, during systole, the coronary arteries are compressed by the contraction of the myocardium, and so, most of the coronary flow occurs during diastole, where the flow increases. In addition, the maximum velocities in both cases are measured in the stenosis throat ($x$ = 0 mm) as observed by other investigations [15, 27].

It is also noted that the velocity is overestimated in the CFD model, particularly at and upstream of the stenosis throat. This is expected because the deformation of the elastic model provides a larger volume for the blood flow through and, as the inflow rate is equal for both models, naturally, the velocities will be higher when a rigid wall is considered. Downstream of the stenosis, at an $x$ coordinate of approximately 10 mm, the velocity is higher for the FSI case, which indicates that the pressure drop at the stenosis creates zones of low pressure, which contracts the artery, and forces the flow to accelerate. In this case, it is thus observed the effect of artery compliance, which allows a steadier supply of flow despite the variable nature of the cardiac cycle.

3.2 Velocity streamlines along the cardiac cycle

After evaluating the velocity profiles, velocity streamlines were created to better visualize and understand how blood flow behaves, as shown in Figure 4.

Figure 4.

Velocity streamlines obtained during systole and diastole for both CFD and FSI simulations.

The results indicate the existence of fluid recirculation downstream of the stenosis for both CFD and FSI models, but the recirculation zones are longer in FSI simulations, and this is due to the considerable vessel expansion driven by the pulsatile blood flow.

3.3 WSS and its indices

The magnitude of the WSS predicted by FSI and rigid model along the artery wall for both systole and diastole are compared in Figure 5. In this case, it is observed that the differences between compliant and rigid-wall models are remarkable. In CFD simulations, WSS values estimated in the stenosis throat are approximately twice of those obtained with FSI simulations. This indicates that the WSS distributions were substantially affected by arterial wall compliance, which is in agreement with previous research [25, 26, 41].

Figure 5.

Wall shear stress profiles obtained along artery wall at (a) systole and (b) diastole.

Taking into account that WSS-related hemodynamic parameters, such as OSI and TAWSS, play an important role in atherogenesis, these were also evaluated. Figure 6 depicts the TWASS profiles obtained in both CFD and FSI simulations.

Figure 6.

Time-averaged wall shear stress profiles obtained in CFD and FSI simulations.

The results evidence high values of TAWSS at the stenosis throat, due to flow acceleration and high-velocity gradient near the wall, and, as expected the TAWSS is overestimated by the CFD model as previously explained for WSS distributions. In spite of these observations, the overall TAWSS distributions for the FSI and rigid-wall cases are identical.

Regarding the OSI profiles depicted in Figure 7, it can be observed that for both CFD and FSI simulations, the maximum values ($\approx 0.5$) are obtained downstream of the stenosis, which indicates the presence of unsteady and oscillatory flow, commonly associated with higher susceptibility to atherosclerotic plaque development [27, 40]. Nonetheless, although the OSI profiles for the two cases look similar and unaffected by wall distensibility, in the distal region, OSI values for the FSI case are slightly higher than for the rigid wall [25]. These differences may be due to the occurrence of longer recirculation areas with the elastic model.

Figure 7.

Oscillatory shear index obtained in CFD and FSI simulations.

In general, it was found that for both rigid and compliant models the post-stenotic region presents lower TAWSS and higher OSI values, which constitute risk factors for the incidence and abnormal plaque formation [16, 27].

3.4 Displacement

Figure 8 represents the arterial wall displacement contours of the FSI simulations during systole and diastole.

Figure 8.

Displacement profile at the arterial wall obtained in CFD and FSI simulations.

In the first place, it is noted that the displacements are similar for both cycle phases. This is a consequence of the assumption of the constant outlet pressure, which is a limitation of this work. In this case, the deformations are slightly bigger in the region upstream of the stenosis, and as there is a pressure drop in the throat, the displacements are somewhat lower in the downstream region.

It is also noteworthy that the displacements of the arterial wall are approximately 1.5–2% of the vessel thickness, which is in agreement with the hypothesis that arteries become stiffer with the development of atherosclerosis. Despite this order of magnitude of the displacement values, these still bring considerable differences in the calculations of the velocity and WSS-related variables.

4. Conclusions

The study of blood flow in stenotic arteries has been made mainly assuming the artery’s wall as rigid, however, in reality, they are naturally elastic. Given that there are some inconsistencies in the literature regarding the comparison between CFD and FSI simulations, in the present work, a comparative study between FSI and CFD modeling was performed in order to investigate the influence of artery compliance on stenotic coronary artery hemodynamics and wall shear stress distribution. The main conclusions obtained from this work are:

• Comparing the rigid and compliant models, the velocities predicted differ slightly.

• The difference between WSS profiles was remarkable. The CFD simulations overestimate the WSS values, which consequently indicates that when a more realistic WSS estimation is needed, is essential to consider the effect of the atrial wall on blood flow.

• Insignificant differences were verified in the TAWSS and OSI measurements.

Acknowledgments

This work has been supported by FCT—Fundação para a Ciência Tecnologia within the R&D Units Project Scope: UIDB/00319/2020, UIDB/04077/2020, UIDB/04436/2020, and NORTE-01-0145-FEDER-030171 (PTDC/EMD-EMD/30171/2017) and EXPL/EME-EME/0732/2021, NORTE-01-0145-FEDER-029394 funded by COMPETE2020, NORTE 2020, PORTUGAL 2020. Violeta Carvalho, Diogo Lopes and João Silva would like to express their gratitude for the support given by FCT through the PhD Grants SFRH/UI/BD/151028/2021, SFRH/BD/144431/2019, and SFRH/BD/130588/2017, respectively.

Conflict of interest

The authors declare no conflict of interest.