Influence of Blood Rheology and Turbulence Models in the Numerical Simulation of Aneurysms

1. Introduction

An aneurysm is a vascular malformation, than can be classified according to its location (cerebral, aortic); and according its shape (saccular, fusiform, and mycotic). The study of blood flow interaction with aneurysms has gained attention from physicists and engineers in recent years [1]. Computational fluid dynamics (CFD) has allowed the use of pacient specific aneurysm geometries [2], substituting the numerical simulation and experimental measurements on idealized models [3].

Furthermore, CFD solutions allow the evaluation of hemodynamics parameters, like wall shear stresses (WSS) and oscillatory shear index (OSI); as potential variables to describe the origin an evolution of aneurysms [4]; and in recent years these parametes have been linked with the aneurysm rupture risk in combination with morphological parameters (size, shape and location) [5]. Therefore, it is of great importance to assure that the numerical results adequately describe the biological part. Selecting laminar or turbulent flow and the rheological model of the blood are factors that have a great impact on the numerical solutions.

The Reynolds number of blood flow in intracranial aneurysms is aroung 400 and a mean value of 600 for abdominal aneurysms [6]; therefore, according to the classical pipe theory the flow is laminar. Lee et al., [7] performed numerical simulations considering laminar flow to evaluate the size and shape of the aneurysm. Han et al., [8] also employed laminar flow in patient specif cerebral aneurysms for rapid evaluation and treatment planning. However, Poelma et al., [9] found that the flow in aneurysms is turbulent due the oscillatory nature of the blood flow. MoMoradicheghamahi et al., [10] found that turbulent models overestimated the time-averaged wall shear stress in stenosed carotid arteries, estimating less damage in the wall than the real value.

On the other hand, it is well known that blood is a non-Newtonian fluid and its rheological behaviour depends on the diameter of the vessel, the Fåhraeus–Lindqvist effect [11]. According to the characteristic length of abdominal aneurysm, the blood can be considered as a Newtonian fluid. For intracranial aneurysms, the mean arteries diameter is around 1 [mm], and the rheology is in the limit between Newtonian and non-Newtonian. Abbasian et al., [12], studied the influence of rheological models in atherosclerotic coronary arteries, considering the most common models used in numerical simulations are Casson, Carreau, and power law. They found that Carreau, modified Casson, and Quemada viscosity blood models were the most accurate; however, the results shown that blood flow can be considered laminar in some sections with high WSS. Abugattas et al., [13] analyzed the influence of the power-law, Cross, and Carreau–Yasuda models in the carotid artery flow pattern, their results shown that the power-law estimates lower WSS in comparison to the other models.

A third effect that is relevant in the blood flow simulation in aneurysms is the elastic properties of the wall. Most of the numerical simulations employed the rigid wall considerations ([14,15]); however, the cerebral arteries are deformable, therefore the FSI must be considered ([16,17]) . According to Moradicheghamahi et al., [10] the elastic wall has a major effect on the CFD analysis than the rheology or laminar/turbulent flow.

The use of turbulent and rheological models to numerically describe the flow in intracranial aneurysms is still under debate, as well as the rigid or elastic artery wall. In the present work a comparison of a turbulent model ($k-ϵ$ model) and laminar model is done for three aneurysms, two aortic aneurysms (AAA1 and AAA2) and a cerebral aneurysm (IA). Four rheological models are evaluated: newtonian, power-law, Carreau and Casson. Finally, the rigid wall solutions are compared with the aneurysms with elastic wall. Wall shear stress (WSS), time average wall shear stress (TAWSS), oscillatory shear index (OSI), relative residence time (RRT), and gradient oscillatory number (GON) were calculated for each case, according to the definitions given by [18], and the results were compared with the laminar case.

2. Materials and Methods

To assess the influence of blood rheology, turbulence models, and wall elasticity on the fluid dynamics in aortic and intracranial aneurysms numerical simulations were performed. Patient-specific cases were employed to rebuilt the geometries.

2.1. Governing equations

With the geometry defined, a CFD simulation implemented in OpenFOAM [19] software was used to solve the mass and momentum conservation equation,Eq. 1 and Eq. 2 respectively, under the assumptions of blood as a homogeneous incompressible fluid and taking into account that the fluid flow is isothermal. Therefore, the Navier-Stokes equations can be written as:

$$\frac{\partial }{\partial t}\left(\rho \right)+\frac{\partial }{\partial x_i}\left(\rho u_i\right)=0$$

(1)

$$\frac{\partial }{\partial t}\rho u_i+u_j\frac{\partial }{\partial x_j}\rho u_i=-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}$$

(2)

Where $u_i$ is the velocity vector, $\rho$ is density, t is time, and p is pressure field. To evaluate the influence of blood rheology, the Newtonian, Carreau, Casson, and Power-law rheological models were used to describe the behaviour of blood viscosity ($\mu$) in two abdominal aortic aneurysms (AAA) and an intracranial aneurysm (IA).

Newtonian fluid is characterized by a linear relationship between shear stress and shear rate. The relation is ${\tau }_{ij}=-\mu {\dot{\gamma }}_{ij}$, where ${\tau }_{ij}$ is the stress tensor, $\mu$ is the fluid viscosity, and ${\dot{\gamma }}_{ij}=\nabla u+\left({\nabla u}^T\right)$ is the rate-of-strain tensor. On the other hand, the simplest non-Newtonian model is the power-law, where the relation between the shear stress and shear rate is non-linear, defined as:

$${\mu }_{eff}=K{\dot{\gamma }}_{ij}^{(n-1)}$$

(3)

Where n is the power-law index, and K is the consistency index. For $n>1$, Eq. 3 describes a dilatant fluid; meanwhile, for $n<1$ a pseudoplastic or shear-thinning fluid is obtained. In shear-thinning fluids, like blood, the shear rate increases while the apparent viscosity decreases. The drawback of this model is that it can overestimate or underestimate the viscosity at very high or very low values of strain rates, respectively, which may not have a physical sense.

In order to describe properly the shear-thinning blood behaviour, two models had been defined. The first one is the Cassson model, given by:

$${\mu }_{eff}={\left(\sqrt{m}+\sqrt{\frac{{\tau }_c}{{\dot{\gamma }}_{ij}}}\right)}^2$$

(4)

Where ${\tau }_c$ is interpreted as the yield stress and m is a viscosity consistency. The second one is the Carreau model, which provides a smooth transition between the viscosity limit values, given by the equation:

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

(5)

${\mu }_{\infty }$ and ${\mu }_0$ are the minimum and maximum viscosity limit values obtained experimentally. ${\mu }_0$ is the viscosity value when the blood is not subjected to shear stress, and ${\mu }_{\infty }$ is the minimum viscosity value reached after a large strain-rate (>100) is applied. Commonly, the ${\mu }_{\infty }$ value is used as the Newtonian constant value. The non-Newtonian behavior is described by the time $\lambda$ and the exponent n, where $n<1$. The parameters used in the rheological models are shown in the Table 1.

2.2. Transport equations for the standard $k-ϵ$ model.

The $k-ϵ$ model was used to study the influence of turbulence on the hemodynamic in aneurysms. The model solves the equations of kinetic energy (6) and the dissipation of kinetic energy (7) to determine the turbulent viscosity (8).

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left[(\mu +\frac{{\mu }_t}{{\sigma }_k})\frac{\partial k}{\partial x_j}\right]+P_k-\rho ϵ$$

(6)

$$\frac{\partial }{\partial t}\left(\rho ϵ\right)+\frac{\partial }{\partial x_i}\left(\rho ϵu_i\right)=\frac{\partial }{\partial x_j}\left[(\mu +\frac{{\mu }_t}{{\sigma }_{ϵ}})\frac{\partial ϵ}{\partial x_j}\right]+C_{1ϵ}\frac{ϵ}{k}{P}_k-C_{2ϵ}\rho \frac{{ϵ}^2}{k}$$

(7)

$${\mu }_t=\rho C_{\mu }\frac{k^2}{ϵ}$$

(8)

Where k is the kinetic energy, $ϵ$ is the dissipation of k, $P_k$ is the generation of kinetic turbulence by the mean velocity gradients; ${\sigma }_k$, ${\sigma }_{ϵ}$, $C_{1ϵ}$ and $C_{2ϵ}$ are constants. The constant values have the following default values: ${\sigma }_k=1.0$, ${\sigma }_{ϵ}=1.3$, $C_{1ϵ}=1.44$, $C_{2ϵ}=1.92$$C_{\mu }=0.09$, according to [10].

2.3. Fluid-Solid Interaction (FSI)

In order to study the effect of the rigid artery wall simplification on the results, a FSI simulation was developed in AAA2 where the linear elastic model was chosen to model the elastic behavior of the artery wall. Bazilevs et al. [24] and Torii et al. [25] concluded that the linear elastic model is a good approximation to model the wall elasticity, after compare differents wall models.

Mechanical properties of the wall tissue are considered uniform, and were based on Isaksen et al. [26]: ${\rho }^s=1000$$kg/m^3$ , $E=1$$MPa$, ${\nu }^s=0.45$, and wall thickness of 2 mm.

The solid momentum equation with the linear elastic consideration is written as follows:

$${\rho }^s\frac{{\partial }^2\delta }{\partial t^2}-\nabla ·[(2{\mu }^s+{\lambda }^s)\nabla \delta ]=\nabla ·[-({\mu }^s+{\lambda }^s)\nabla \delta +{\mu }^s{\nabla }^T\delta +{\lambda }^{s}tr(\nabla \delta )I]$$

(9)

Where ${\rho }^s$ is density, $\delta$ is displacement tensor, ${\mu }^s$ and ${\lambda }^s$ are the Lame constants, and I is the second order tensor.

2.4. Boundary conditions

For the reological and turbulence analysis, the artery wall was assumed to be rigid, and the no-slip boundary condition was applied. A pulsatile flow was defined at the inlet, using the waveform from Banerje et al. [27], with a mean Reynolds value of 600 in AAA, and 400 in IA. In the outlet branches, a pressure value of 0 $Pa$ was imposed. Density was consider constant with 1056 $kg/m^3$ value. Viscosity ($\mu$) depends on the rheological model. Six cardiac cycles were used to do the transient simulations, as the condition for the existance of the periodic solution, and the results of the last cycle were analyzed.

For FSI case, the inlet and outlet boundary conditions of the rigid wall case were employed. The displacement of the fluid and solid interfaces must be the same (Eq. 10), and the traction in both domains must be in equilibrium (Eq. 11). Furthermore, the fluid obeys the non-slip condition on the solid wall. The solid sides are fixed, zero displacement, and a pressure of 0 Pa was imposed on the external wall.

$${\delta }_s={\delta }_f$$

(10)

$${\sigma }_s·{\widehat{n}}_s={\sigma }_f·{\widehat{n}}_f$$

(11)

The time step was variable to maintain a Courant number less than 0.1, for numerical stability. A second-order scheme was used to discretize the governing equations.

2.5. 3D Model

The patient-specific aneurysms geometries were extracted from digital subtraction angiography (DSA) images, using the open-source software 3D Slicer. To isolate the area of the aneurysm and the parent vessel, the threshold tool was used to distinguish blood from the rest of the tissues. The artery and the tissue around the artery were removed, the 3D model was exported to STL format, which was used to mesh the analyzed geometry in the simulation. The aortic aneurysms (AAA1 and AAA2) and cerebral aneurysm are shown in Figure 1

2.6. Grid generation

Mesh generation of the computational models was performed using the OpenFOAM’s tool snappyHexMesh. The mesh was hexahedral dominant, and a prismatic boundary layer was implemented to capture accurately the influence of the wall and the WSS.

Three cases were studied: two abdominal aortic aneurisms, (AAA1) and (AAA2), and an intracranial aneurysm (IA), Figure 1. A mesh independence test was developed in each model. In the abdominal aortic aneurysm (AAA1), three different meshes were created, with 460k (coarse), 620k (medium), and 935k (fine) elements. In order to compare the results of these three meshes, the area-average WSS of the last cardiac cycle was analyzed and plotted (Figure 2). The maximum difference of the fine and coarse mesh was 5 %, and respect with medium mesh, the maximum difference was 3 %. In the second abdominal aortic aneurysm (AAA2), the meshes studied comprised 525k (coarse), 670k (medium), and 1200k (fine) elements. The maximum difference of area-average WSS between fine and coarse meshes was 8 %, and in comparison with medium mesh the maximum difference was 3%. In the last model, the intracranial aneurysm (IA), the meshes studied comprised 200k (coarse), 300k (medium), and 2000k (fine) elements. The maximum differences comparing the fine mesh with the coarse and medium mesh were 7% and 4%, respectively. The meshes selected to do the rheological and turbulence analysis were the medium meshes in all the models, that have a percentage error less than 5%.

A wall turbulence model was employed, for the three cases $Y^{+}=2$ was the maximum value obtained, meanwhile the maximum average value was $Y^{+}=0.25$. Therefore, the turbulent boundary layer can be resolved correctly. Besides, the mesh for the solid artery wall in the FSI study was achieved by extruding the fluid wall, with a constant thickness of 2 mm, built by 5 layers with an expansion ratio of 1.2. This method allows the match between the fluid and solid meshes.

3. Results

Velocity, wall shear stress (WSS), time average wall shear stress (TAWSS), oscillatory shear index (OSI), relative residence time (RRT), and gradient oscillatory number (GON) were calculated, according to the definitions given by [18]. The WSS was determined for systole and diastole; meanwhile, the other parameters where calculated for the last cardiac cycle. In this section the most relevant Figures are shown, the repository [28] shows the figures of all cases.

The indices for IA, Figure 3, shows higher values of WSS for systole, being the dominant effect on the TAWSS. The circled areas in the figure indicate the zones of low WSS and high values of OSI, RRT and GON. The results of IA, which have an aspect ratio (AR) of 1.13, shows that the blood flow goes upstream through the parent artery (ICA), impacting the distal side of the neck, then the flow goes towards the aneurysm dome, and finally returns to the parent vessel by the proximal side of the neck, completing one vortex, as was described by Munarriz et al. [29]. The highest value of WSS was found in the distal neck because of the impinging flow, while the values of WSS in the dome depend on the velocity distribution in the sac. Meng et al. [30] proposed that the growth and rupture of an aneurysm is associated with high or low values of WSS, depending on aneurysm shape and flow behavior. The numerical results shows that the IA belongs to the slow recirculation class, according to [30]. High values of OSI and low values of WSS can trigger inflammatory-cell-mediated destructive remodeling that leads to large, thick wall aneurysm phenotype.

In the aortic aneurysms (AAA1 and AAA2), Figure 4, the flow descends through the aorta to the aneurysm zone, where multiple vortices develop; furthermore, the results exhibit some areas with recirculation and stagnation flow. The circled areas meet the condition of low WSS and high OSI and RRT. Almost all the proposed indices points areas with high probability of destructive wall remodeling, according to [30]; except GON, that presented zones of high values in different parts of the aneurysm and its influence it is not clear.

3.1. Rheology

The strain rate results, Table 2, were calculated at the artery wall, at the aneurysm wall, at the blood volume in the artery (artery body), and finally at the blood volume in the aneurysm (aneurysm body). The strain rate values at the wall for IA were significantly higher that the AAA cases. As was expected, the strain rates in the fluid were lower than the values at the wall; however, IA cases presented values an order of magnitude bigger than the AAA. The values are relevant in order to determine the rheology influence on the flow behaviour.

The IA strain rate values were higher than the AAAs. The IA strain rate values in the total artery wall and the aneurysm wall were the order of ${10}^3$. While the strain rate values in the total artery body and only the aneurysm body were highest than ${10}^2$. AAAs values in the wall are near to 100, which is the critical strain rate value where the viscosity starts to behave like a constant Newtonian value in the no-Newtonian models. Values in the body far away from the wall are low, and differences between Newtonian and non-Newtonian models are expected.

The results of the rheological analysis in intracranial aneurysms show minimal differences between the models, considering the different indices under study with similar qualitative results between the Newtonian and non-Newtonian models. This can be explained because the strain rate is larger than 100 in the body and the aneurysm wall, as shown in the Table 2. AAAs results shown greater differences than IA in the studied indices. AAA2 was the most critical case, showing greater differences that could be attributed to low values of strain rate, even in the wall.

Figure 5 shows the Oscillatory Share Index (OSI) contours for the Newtonian and the non-Newtonian models, qualitative differences in the patterns can be observed. Below the OSI contours, (A to D), the differences of each of the non-Newtonian vs Newtonian models are shown, (E to G), some areas with differences of up to 100% can be seen. Similar differences were found in WSS, TAWSS, RRT and GON, repository.

In aortic aneurysms, far away from the wall the strain rate values are low, and the viscosity values differ from the Newtonian case. Therefore, the velocity profiles differ between newtonian and non-Newtonian models, as can be seen in the Figure 6, where the magnitude of velocity of the Newtonian and Carreau models are compared. The Newtonian model presents higher velocity values, and 25% of the plane area had a 50% difference or higher.

3.2. Turbulence

Turbulence is a complicated phenomenon, that is characterized by a caotic and instable flow. It is well known that for flow in pipes the critical Reynolds number ($Re$) is $\approx 2300$. For lower $Re$ a laminar flow can be considered, and for upper $Re$ turbulent behaviour starts. In the present study a pulsatile flow was taken into account [27], and according to the physiological conditions a mean Reynolds number ${Re}_{mean}\approx 600$ was used for AAA, with variations of ${Re}_{min}\approx 35$ to ${Re}_{max}\approx 1367$, with a Womersley number $Wo\approx 20$. For IA the following data were taken into account: ${Re}_{mean}\approx 400$, ${Re}_{min}\approx 23$, ${Re}_{max}\approx 911$, and a $Wo\approx 4$.

According to the classical theory in pipes, the blood flow in the aneurysms is laminar. However, Poelma et al., found differences in the velocity fields between the cycles under analysis Poelma et al. [9]; this changes are a characteristic of turbulent flows. The geometry of the aneurysms and the pulsatile nature of the flow are perturbations of the laminar flow and could modified the cirtical Reynolds number. Furthermore, the deceleration in the cardiac cycle is a destabilizing factor of the laminar flow; therefore Hershey and Im [31], NEREM and SEED [32], concluded that the critical Reynolds number is inversely proportional to the Womersley number.

To determine if the flow is laminar or turbulent, the results of 4 cardiac cycles were compared for both types of aneurysms. Differences between each cardiac cycle were found for aortic aneurysms, Figure 7; which is a characteristic of turbulent flows Poelma et al. [9]. On the other hand, for intracanial aneurysms no significant changes were observerd between cardiac cycles; therefore, the flow in IA can be considered as laminar.

Numerical simulations for both types of aneurysms were performed considering lamninar flow and Newtonian model, turbulent flow (using the K-epsilon model) and four rheological fluids : Newtonian, Carreau, Casson, and power-law. Figure 8 shows the WSS, TAWSS, OSI, GON, and RRT for an aortic aneurysm (AAA2) of the mentioned cases and the differences between them. The areas where the greatest differences occur between a laminar and a turbulent flow are marked in the figure (red circle). Differences exceeding 100% of the laminar value were found in large areas of the aneurysm, which overestimate the WSS, and consequently the TAWSS (lower part of the aneurysm); and promoting changes in the patterns of OSI, GON, and RTT in the central part of the aneurysm. Small differences (less than 5%) were obtained between the rheological models using the K-epsilon model, suggesting that the turbulent models are more relevant than the rheological characterization of the flow. This same behaviour was obtain for the second aortic aneurysm (AAA1). Meanwhile, for the intracranial aneurysm (IA) no differences between the rheological models nor laminar of turbulent flow were found.

Figure 8. Hemodynamic parameters uder study for A) Laminar - Newtonian, B) $k-ϵ$ - Newtonian, C) $k-ϵ$ - Carreau, D) $k-ϵ$ - Casson, E) $k-ϵ$ - Power Law, F) Differences between Laminar-Newtonian vs $k-ϵ$-Newtonian, G) Differences between Laminar-Newtonian vs $k-ϵ$- Carreau, H) Differences Laminar-Newtonian vs $k-ϵ$-Casson, I, Differences Laminar-Newtonian vs $k-ϵ$-Power Law.

Figure 8. Hemodynamic parameters uder study for A) Laminar - Newtonian, B) $k-ϵ$ - Newtonian, C) $k-ϵ$ - Carreau, D) $k-ϵ$ - Casson, E) $k-ϵ$ - Power Law, F) Differences between Laminar-Newtonian vs $k-ϵ$-Newtonian, G) Differences between Laminar-Newtonian vs $k-ϵ$- Carreau, H) Differences Laminar-Newtonian vs $k-ϵ$-Casson, I, Differences Laminar-Newtonian vs $k-ϵ$-Power Law.

Figure 9. Newtonian value of OSI at different cycles.

Figure 9. Newtonian value of OSI at different cycles.

Figure 10. Oscillatory Share Index results of Newtonian and non-Newtonian models, comparing turbulence models. A) Newtonian laminar B) Newtonian KE C) Carreau KE D) Casson KE E)Power-Law-KE F) Newtoanian laminar-Newtonian KE Differences G) Newtoanian laminar-Carreau KE Differences H) Newtoanian laminar-Casson KE Differences I) Newtoanian laminar-PowerLaw KE Differences

Figure 10. Oscillatory Share Index results of Newtonian and non-Newtonian models, comparing turbulence models. A) Newtonian laminar B) Newtonian KE C) Carreau KE D) Casson KE E)Power-Law-KE F) Newtoanian laminar-Newtonian KE Differences G) Newtoanian laminar-Carreau KE Differences H) Newtoanian laminar-Casson KE Differences I) Newtoanian laminar-PowerLaw KE Differences

3.3. Fluid-Solid Interaction (FSI)

The movement of the elastic arterial wall under the effects of the pulsatile signal results in a different fluid dynamics behaviour in comparison to the consideration of a rigid wall. All the indices analyzed present different patterns, the one that presented the greatest difference with respect to these two considerations of the wall was the GON, where the patterns are completely different and the average of the differences in the entire area of the aneurysm is 43%, while the next index that had the greatest difference was the WSS during diastole with 25%, the third place is occupied by the index most related to intracranial aneurysms, the OSI, with a 20% difference, Figure 11.

The results of all cases are shown in the repository [28].

4. Discussion

Laminar and turbulent flow models were compared using four rheological models and rigid or elastic walls to simulate the flow in two aortic aneurysms and a intracranial aneurysm. According to Pratumwal et al., [33], the Carreau model is the best to describe the blood rheology; however, for IA the rheology has little influence on the flow behaviour. Meanwhile, for AAA the biggest differences are presented by the Carreau model with respect to the Newtonian model, for WSS (25%) and TAWSS (22%), Figure 12. For RRT the effect of the rheological model is greater than the laminar-turbulent flow, Figure 13.

No differences were found between laminar and turbulent flows for WSS and TWSS. Turbulent model underestimate the OSI, RRT, and GON in comparison with the laminar flow. Therefore, it is of great relevance the flow model to simulate the blood flow in aortic aneurysms. Finally, the most relevant effect to take into consideration for numerical simulations is the rigid-elastic wall assumption.