Fluid Flow and Structural Numerical Analysis of a Cerebral Aneurysm Model

Abstract

Intracranial aneurysms (IA) are dilations of the cerebral arteries and, in most cases, have no symptoms. However, it is a very serious pathology, with a high mortality rate after rupture. Several studies have been focused only on the hemodynamics of the flow within the IA. However, besides the effect of the flow, the development and rupture of the IA are also associated with a combination of other factors such as the wall mechanical behavior. Thus, the objective of this work was to analyze, in addition to the flow behavior, the biomechanical behavior of the aneurysm wall. For this, CFD simulations were performed for different Reynolds numbers (1, 100, 500 and 1000) and for two different rheological models (Newtonian and Carreau). Subsequently, the pressure values of the fluid simulations were exported to the structural simulations in order to qualitatively observe the deformations, strains, normal stresses and shear stress generated in the channel wall. For the structural simulations, a hyperelastic constitutive model (5-parameter Mooney–Rivlin) was used. The results show that with the increase in the Reynolds number (Re), the recirculation phenomenon is more pronounced, which is not seen for Re = 1. The higher the Re, the higher the strain, displacement, normal and shear stresses values.

1. Introduction

Intracranial aneurysms (IA) are local dilations of the intracranial arteries; their occurrence in the world population varies from 0.5% to 6% [1,2]. About 0.25% of these IAs rupture and present subarachnoid hemorrhage [3]. The mortality rate from these hemorrhages reaches 60% [4]. The literature indicates how a combination of multiple factors can lead to the onset and development of IA: genetic diseases, fungal infection, hypertension, cigarette smoke inhalation, alcohol intake and just aging [5,6]. These factors result in an alteration of the original arterial wall properties and flow behavior. Factors such as weakening of the vessel wall, hemodynamic changes and inflammatory processes can also trigger the appearance of the aneurysm, which can become a stable dilatation or rupture [7,8,9]. The interaction of these factors and complex processes that lead to aneurysm formation still need to be better understood. Thus, it is essential to study the biomechanical behavior of the blood vessel together with the behavior of blood flow to determine regions more susceptible to rupture.

An important characteristic of blood vessels is the existence of soft tissues, which have a behavior called hyperelasticity [10]. Hyperelastic materials have a mechanical behavior characterized by a high deformation before reaching the breaking strength [11,12]. Due to the risks involved, it is very difficult to analyze in vivo the biomechanical behavior of blood vessels, particularly in the region of the aneurysm [13]. However, it is possible to manufacture flow phantom biomodels for experimental tests with compatible properties and geometries [14,15]. Flow phantoms basically consist of models that mimic vascular geometry for flow analysis (rigid flow phantoms) and flow analysis and arterial wall deformation (flexible flow phantoms). These biomodels were produced in silicone rubber, as they allow optical access to the flow and also reproduce the physiological compliance of the artery [16,17]. In its manufacture, materials whose biomechanical behavior is close to the real biological material were used [18]. A frequently used material is polydimethylsiloxane (PDMS) [19,20]. PDMS has been widely used in different kinds of research fields [21,22,23,24,25,26,27] due to its excellent properties and ability to replicate reliable geometries, in addition to its high flexibility, low glass transition temperature, good chemical resistance, excellent thermal stability and biocompatibility [28].

In addition to the experimental approach, in recent years, there has been an increase in numerical hemodynamic studies using the computational fluid dynamics tool CFD [29,30,31,32,33,34]. This tool has been effective in investigating the relationship of vascular hemodynamics with the mechanisms of initiation, growth and rupture [9,35,36,37]. CFD is also used to evaluate endovascular treatments [38,39] and validate new endovascular devices [40,41]. Most of these studies use Newtonian flow models; however, blood flow does not have this type of behavior, and there is no consensus on which type of viscosity model should be adopted to model blood flow in IA [37,42]. Some studies consider that blood viscosity follows a Newtonian behavior in cerebral arteries, as shear rates are considered to be higher than the range necessary for non-Newtonian properties to become effective [37,43]. Other studies report that, for medium and large arteries, blood viscosity is practically constant, which justifies, in blood flow simulations, considering it as Newtonian [44]. However, some studies claim that Newtonian viscosity results in excessive predictions of aneurysm wall shear stress, which could compromise the effectiveness of predicting the risk of rupture [37]. Given these controversies about which type of viscosity to use, it is important to compare the type of regime for each application. In this study, numerical simulations were performed using different Reynold numbers (Re), and comparisons were made between Newtonian and non-Newtonian flow (Carreau fluid model).

Although the numerical study using the CFD tool has achieved some advances in the understanding of IAs, they have some limitations, most of these studies use rigid vessel wall aneurysm models, and in reality, there is an interaction of blood flow with the vessel wall. Endothelial cells present in vessel walls detect blood flow through a variety of mechanical cell surfaces and intracellular sensors, any change in blood flow. Endothelial cells initiate different biological processes [45,46] that lead to a dynamic cycle of wall degeneration and subsequent remodeling. Studies have also shown that endothelial cells respond to WSS [47]. Thus, recent studies reported the importance of taking into account the interaction between blood flow and the vessel wall in the model of IA formation, growth and rupture [18,48]. However, there are few studies that precisely characterize this phenomenon of interaction; as such, in order to fill the little information reported in scientific publications, the aim of this work was to numerically determine the state of displacement and stress on the IA wall caused by different flow conditions. Thus, in order to characterize the mechanical behavior of biological tissues under different flow conditions, numerical simulations were performed using the Mooney–Rivlin constitutive model that mathematically describes hyperelastic materials [49].

In this work, numerical simulations were performed by combining both fluid flow and structural analysis [50,51] using Ansys^® software. The fluid flow analysis simulates the blood flow inside the vessels, and the resulting pressure in the vessel wall was analyzed for mechanical deformations occurring in the vessel wall as well as in the aneurysm. The existence or absence of the recirculation phenomenon in different Reynolds numbers (Re) was studied, and it was possible to visualize and analyze the development of the fluid recirculation phenomenon inside the aneurysm. Finally, the numerical characterization of the biomechanical behavior of aneurysms was carried out.

2. Problem Description

In the present study, numerical simulations of blood flow and biomechanical wall’s behavior of an in vitro channel with aneurysm, represented in Figure 1, were performed, resorting to Ansys^® software. The geometry and dimensions of the idealized IA model are based on the work developed by Parlea et al. [52].

Initially, fluid dynamics simulations were performed in order to observe hemodynamic patterns such as flow recirculation inside the aneurysm and pressure gradients for different flow velocities. The fluid flow analysis also allows obtaining the pressure distribution on the channel walls, which are later used as boundary conditions in structural simulations. Structural simulations were carried out in order to analyze the biomechanical behavior of the wall through the values of displacement, strain and stress.

2.1. Fluid Flow Simulations

For flow simulations, both continuity and Navier–Stokes equations were solved using the finite volume commercial software Fluent, Ansys 2020 R2.

The steady and incompressible regime and rigid walls were prescribed. Simulations were carried out considering and discarding the non-Newtonian properties of blood, which permits evaluating the impact of non-Newtonian characteristics of blood in the properties of studied flows. As constitutive equations, Newtonian and Carreau models were used, the last one being described by Equation (1):

$$\eta ={\eta }_{\infty }+({\eta }_{0-}{\eta }_{\infty }){\left[1+{\left(\lambda \dot{\gamma }\right)}^2\right]}^{\left(n-1\right)/2}$$

(1)

where η is the viscosity, η₀ is the zero viscosity, η_∞ is the infinite viscosity, λ is the characteristic time (the inverse of the shear rate for which pseudoplastic behavior begins), n is the power index and $\dot{\gamma }$ is the shear rate. For blood, parameters of Equation (1) assume the values η₀ = 0.056 Pa·s, η_∞ = 0.00345 Pa·s, λ = 3.313 s and n = 0.3568 [26]. In the case of the Newtonian model, blood viscosity considered was 0.00345 Pa s and the density 1050 kg/m³ [26].

2.1.1. Geometrical Domain and Mesh

Initially, the idealized geometry of the aneurysm based on other studies [11] was drawn in the SolidWorks^® CAD software, converted to Parasolid format (.x_t) and then exported to Ansys^®.

For fluid flow simulation, only the fluid domain was considered (Figure 1). The analyzed geometry is a cylindrical channel with 3.7 mm in diameter and 50 mm in length, with a central sphere representing an artery and an aneurysm idealized saccular, respectively.

A mesh consisting predominantly of hexahedral elements was generated, with a total of 170,820 elements and 174,332 nodes (Figure 2). Before choosing the mesh, a preliminary mesh test was performed. A mesh with twice the number of elements was created, keeping the same spatial distribution, and then the simulations were performed. The maximum velocities values obtained by both meshes were observed, and a difference of less than 1% was found in the results, which shows that the results are independent of the mesh used.

2.1.2. Boundary Conditions and Solver

The simulations were carried out considering atmospheric pressure as the operating pressure. As boundary conditions, four distinct mean velocities were imposed in the inlet of the channel: 0.9 mm/s, 88.8 mm/s, 444 mm/s and 888 mm/s. These velocities were chosen in order to obtain the corresponding Reynolds numbers (Re): 1, 100, 500 and 1000, respectively, for Newtonian flows. These Re values were chosen in order to verify the occurrence or not of fluid recirculation and if this phenomenon would increase according to the growth of Re values. At the structural level, observe whether the deformations were proportional to the increase in Re.

The wall of the channel was assumed rigid, and as the wall velocity is zero, the non-slip condition was applied at the wall. Additionally, the gauge pressure (0 Pa) was imposed in the outlet of the channel.

The pressure-based solver was used to solve numerical integrations, as this is the most suitable for problems with incompressible fluids. Speed–pressure coupling was performed using the PISO scheme. The spatial discretization of the pressure equation was performed with the PRESTO! approximation, while the moment equations were discretized with the QUICK scheme.

2.1.3. Model Validation

In order to validate CFD calculations, simulations were carried out for the conditions referred before in a channel similar to the one presented above but lengthier in the section before the aneurysm to guarantee a fully developed flow in that region. Two tests in that section were performed. First, velocity profiles were compared with the analytical solution for a Newtonian fully developed flow in a circular pipe which can be determined by [53]:

$$v\left(r\right)=2u\left(1-{\left(\frac{r}{R}\right)}^2\right)$$

(2)

where v is the velocity, u is the mean velocity, R is the radius of the channel and r is the radial position.

In Figure 3, it is possible to verify the good agreement between numerical results and analytic profile (mean relative error range between 2.4% and 2.6% for the three velocities), which concludes that the numerical model predicts the local properties of the studied flow well.

The second test consists of the comparison between the obtained pressure drops and predictions from the Haggen–Poiseuille equation [54]:

$$\Delta P=\frac{8\eta Lu}{R^2} \iff \frac{\Delta P}{L}=\frac{8\eta u}{R^2}$$

(3)

where ΔP represents the pressure drop and L is the channel length. Once again, a good agreement between pressure drop obtained by Equation (3) and numerical results was found, as shown in

Table 1. Pressure drop for fully developed Newtonian flows in the region before the aneurysm.

Re	(ΔP/L)Equation (3) (Pa/m)	(ΔP/L)numerical (Pa/m)	Relative Error (%)
100	716.11	713.79	0.33
500	3580.54	3575.83	0.13
1000	7161.08	7364.60	2.84

.

2.2. Structural Simulations

The structural simulations were performed using the static structure of the Ansys software. Initially, new material was inserted into the Engineering Data. The material considered was Polydimethylsiloxane (PDMS) Sylgard 184, as it is a hyperelastic material, as are the arteries. A stress/strain curve of the PDMS obtained by tensile testing was then inserted into the properties and approximated by a hyperelastic model previously defined in Ansys, Mooney–Rivlin parameter 5 model, which is a constitutive model adequate for this kind of simulation [55]. The Mooney–Rivlin constitutive law is described by Equation (4):

$$\Psi =C_{10}\left(I_1-3\right)+C_{01}\left(I_2-3\right)+C_{20}{\left(I_1-3\right)}^2+C_{11}\left(I_1-3\right)\left(I_2-3\right)+C_{02}{\left(I_2-3\right)}^2$$

(4)

where $\Psi$ is the strain energy density function; $I_1$ and $I_2$ are the strain constants; and $C_{10}$, $C_{01}$, $C_{20}$, $C_{11}$ and $C_{02}$ are the material parameters, which are determined using the data obtained experimentally by tensile test.

In Figure 4, the approximation of the mentioned curve can be seen.

The geometry for this simulation represents only the channel wall that has 0.15 mm of thickness.

As the geometry for the flow simulations, this geometry was previously drawn in Solidworks and later exported to Ansys. A tetrahedral mesh with an element size of 0.5 mm was generated, with a total of 48,931 nodes and 24,370 elements. This mesh can be seen in Figure 5a. Next, the boundary conditions were defined, fixed support was applied to the sides of the model, and then the pressure imported from fluid simulations in Fluent was applied to the inner wall of the channel. Figure 5b illustrates the pressure applied to the channel wall.

Finally, the software was configured to calculate the solutions for displacement, strain, the normal and shear stress.

3. Results

In this section, the results of the fluid flow simulations for different velocities and the results of the structural simulations are presented.

3.1. Fluid Flow Analysis

Fluid flow simulations were carried out considering the positive z-axis as the main flow direction, as shown in Figure 2.

3.1.1. Flow Patterns

The hemodynamic pattern was analyzed, and similarities were found between all the studied flows (distinct inlet velocities and different rheological models). Figure 6 illustrates the flow behavior for the four Re analyzed in the central plane at x position (x = 0). It is clear that for Re = 1, the flow makes a smooth curvature, bypassing the geometry of the aneurysm wall; in this case, there is no recirculation. However, for the other Re, which has much higher flow velocities, recirculation inside the aneurysm is noted, which is more pronounced as the number of Re increases.

By observing the streamlines (Figure 7), two distinct sections can be considered concerning flow behavior. In a region close to the inlet and outlet of the channel, the typical behavior of laminar flow in a circular pipe was observed. However, when we move closer to the bulge, the flow suffers a deviation in the direction of protuberance and vortex formation inside it was observed. The last phenomenon is more pronounced for the higher velocities, as already shown.

In Figure 8, it is also possible to verify the behavior described before. Lines 1 and 3 exhibits the typical shape of velocity profiles in a circular pipe, that is, symmetric profiles about the center of the channel. The shape of these two profiles differs since the region before the aneurysm is clearly a region of developing flow, and after the aneurysm, the parabolic profile appears since flow is closer to a developed flow. This behavior is the reason for using a lengthier channel for the model validation. It should be noted that for lower velocity (88.8 mm/s), the difference presented and discussed in this paragraph was not found, and two parabolic profiles were observed.

The profile represented by line 2 in Figure 8 illustrates well the lower velocities in the aneurysm. Moreover, it is possible to verify a little shift of the maximum velocity in the bulge direction.

By comparing the obtained results for the same velocity inlet, considering and discarding the non-Newtonian properties of the blood, the differences observed in velocity’s magnitude could be neglected, which is the reason why the authors decide to present only an example of profiles and streamlines.

3.1.2. Pressure

As observed for velocity, pressure distribution was also revealed to be almost independent of non-Newtonian properties of the blood, as shown in Figure 9b. However, pressure is always higher when blood rheology is described by the Carreau model.

In Figure 9a, static pressure along the channel is presented for the Carreau model. As expected, pressure decrease in the flow direction, except in the region of the aneurysm where pressure remains almost constant due to the increase of the cross-section of the channel and velocity close to zero.

Pressure drop for all studied flows was determined and presented in

Table 2. Pressure drop for both rheological models.

Re	ΔPCarreau (Pa)	ΔPNewtonian (Pa)	ΔPratio
100	46.73	37.64	1.24
500	290.41	274.81	1.06
1000	726.12	704.84	1.03

. Analyzing the ratio between ΔP obtained, using Carreau model and Newtonian model as constitutive models, a reduction in this ratio with the increase in Re, although they are all very close to 1.

The influence of mean velocity in the ratios of pressure drop, as well as pressure, could be explained by blood viscosity obtained with the Carreau model. As shown in Figure 10, blood viscosity increase with the decrease in mean velocity and assume values higher than the Newtonian viscosity (0.00345 Pa.s), which leads to higher pressure drops.

As referred before, flow simulations were performed in order to find the blood pressure in the wall (Figure 11) and use it as a boundary condition in the wall structural analysis.

In Figure 12, static pressure obtained along the black line of Figure 11 is presented. Once again, the decrease in pressure in the main flow direction was observed, except in the region of the aneurysm. Here the pressure remains constant along with the bulge and exhibits pronounced increases/decreases in the transition of the circular channel to the aneurysm.

The minimums observed in Figure 12 are in the transition from the circular duct to the bulge, where blood viscosity obtained with the Carreau model moves closer to the Newtonian viscosity, as could be verified in Figure 10.

3.2. Structural Analysis

The results of the channel wall structural simulations are presented below.

By observing Figure 13 below, it is possible to observe that the displacement along the channel occurs in a similar way; in a qualitative perspective, this effect can also be justified by the fixation support applied on the sides of the geometry during the modeling of the problem, which causes that the total displacement values are null at the ends, but do not affect the results in the aneurysm. Displacement values increase as the Reynolds number increases.

However, it is clear that at the dome of the aneurysm, there is a change in the behavior of total displacement for each Re. The total displacement becomes more and more intense on the sides of the aneurysm (areas of great effort) and is reduced at the top of the aneurysm.

Figure 14 shows the strain values. For these results, it is notorious both in the channel and in the dome of the aneurysm the increase of strain as the inflow velocity is increased. The highest strain values are concentrated in the transition zones, from channel to aneurysm and from the aneurysm to channel.

In addition to displacement and strain, normal stress and shear stress values were also computed. Next, the images of shear stress (Figure 15) and normal stress (Figure 16) for the highest Re (Re = 1000) and the graph of normal stress for all analyzed Re (Figure 17) are presented.

The behavior of the shear stress changes according to the coordinates, with the change close to the aneurysm being more noticeable (Figure 15). For shear stresses analyzed in component XY, the highest stress values are in the channel zone below the aneurysm dome. In the YZ component, the highest shear stress values are found in the channel along the curvature of the aneurysm outlet. Additionally, in component XZ, the highest values are near the neck of the aneurysm. By analyzing the mean of the values in each direction, we have the highest mean shear stress in the YZ component, which is the flow direction within the channel.

Figure 16 shows the results for normal stress. As can be seen, the normal stress is more accentuated in the Z direction. The higher stress values are close to the aneurysm entrance and exit. When we analyze the graph in Figure 17, with the mean stress values, it is clear that the normal stress in Z, for all analyzed Re, has a greater relevance (higher values) when compared to the other directions (X and Y).

4. Conclusions

The present work aimed to observe the flow behavior in an ideal intracranial aneurysm using a Newtonian blood viscosity model and a non-Newtonian model (Carreau model). In addition, the biomechanical behavior of the wall was analyzed.

By analyzing the distribution of velocities in the two rheological models, it was observed that there was no significant difference in the magnitude of velocities. Regarding the pressure distribution, despite the values being higher in the blood rheology described by the Carreau model, this also proved to be almost independent of the non-Newtonian properties of the blood. Hence, this result indicates that assuming the characteristics of blood, being Newtonian fluid, is acceptable.

In order to obtain more insights into the hemodynamic behavior, different Re was compared, from Re = 1 (without physiological relevance) to Re = 1000, in order to observe the existence or not of recirculations inside the IA. With the obtained results, it was evident that the increase in the number of Re makes the phenomenon of blood recirculation more noticeable. There is only the formation of a vortex in the center of the aneurysm, which becomes increasingly accentuated closer to the wall. Nevertheless, there is no formation of additional vortex regions as the Re number increases. For Re = 1, the phenomenon of blood recirculation does not appear due to the extremely low flow velocities at the channel inlet.

Regarding the structural analysis, it was observed that the displacement was maximum on the sides of the aneurysm and in the passage from the channel to the aneurysm and from the aneurysm to the channel, which are the areas where there is an abrupt change in geometry. Likewise, an increase in the Re leads to an increment of both normal stress and shear stress values. However, for Re = 1, the stress was extremely low, with values close to zero.

The normal stresses have shown more accentuated values at the Z direction, and the shear stresses presented the highest values at the XY plane. Hence, the normal stresses were higher at the transition region from the vessel to the aneurysm. Additionally, the shear stresses were higher at the sides of the aneurysm. Note that the representation of the type of distribution of the normal and shear stresses was similar to the tested Reynolds numbers.

Although in the present work, a simplified vessel model was used, this does not withdraw the validity of results since the form and geometry of the blood vessel vary from patient to patient. The use of a more realistic vessel, for example, constructed from 2D medical images obtained by computed tomography, gives more precise information regarding this patient. In the near future, it is intended to perform this comparison.