Characterizing pulsatile blood flow in a specific carotid bifurcation: insights into hemodynamics and rheology models

1.   Introduction

The carotid artery is one of the most important arteries in the human body [1]; its primary function is to supply blood to the brain through the basilar trunk. Its complex structure has many curved areas, particularly the carotid sinus. Furthermore, at the carotid branches, there is a siphon that forms a highly curved spiral. Under real physiological conditions, hemodynamic factors (blood pressure, blood flow, flow velocity, etc.) in this complex anatomy contribute to the appearance of atherosclerosis [2]–[4]. The latter is strongly linked to stenotic obstruction (reduction of the section of the blood vessel). In addition to differential pressure, bifurcation geometry, and arterial wall properties, several studies confirm that hemodynamic factors such as flow distribution, recirculation, low and oscillatory wall shear stress play a significant role in the development and progression of atherosclerotic plaques and other arterial lesions [5]–[7]. Therefore, thoroughly understanding the carotid artery's blood flow seems crucial.

Accessing biological flows poses a significant challenge in the study of hemodynamic phenomena. In recent years, computational fluid dynamics (CFD) has proven to be an effective solution to study similar problems and understand hemodynamics using mathematical tools [8],[9]. Numerical simulation of constitutive models can accurately describe the rheological blood response under physiological flow conditions [10]. It is recognized as an invaluable tool for interpretation and analysis of the functionality of the circulatory system in physiological and pathological situations [11],[12].

In this context, the present work studies blood flow in a specific carotid bifurcation. The flow is considered incompressible, laminar in the first case, and turbulent modeled with the k-ϵ model in the second case. The carotid wall is assumed to be rigid. Two different forms of physiological pulses from the literature were applied as boundary conditions to generate blood movement. These pulses are summarized in two pressure-velocity states: P₁-V₁ (theoretical) and P₂-V₂ (actual). Three viscosity models were used to model the non-Newtonian blood flow: Cross, Carreau, and Quemada. The resolution of the governing equations system is carried out using the Fluent solver based on the finite volume method. The calculations are performed over a total time of 2.4 s with a period of 0.8 s.

2.   Mathematical modeling

2.1.   Problem description

The study investigates the blood flow characteristics in a specific rigid carotid bifurcation. (Figure 1, Table 1). It consists of two main parts, separately examining laminar and turbulent flow regimes. The laminar flow case considers incompressible blood with a 1050 kg/m³ density. In the turbulent flow case, the k-ϵ model is employed.

The rheological behavior of blood is analyzed by comparing the Newtonian case, where the blood viscosity is constant at 0.0035 Pa.s, with the non-Newtonian instances described by the Cross (1), Carreau (2) and Quemada (3) models [13]. These models capture the non-Newtonian characteristics of blood flow.

Physiological pulsatile flow is modeled using two pairs of periodic boundary conditions: P₁-V₁ (theoretical) and P₂-V₂ (real). The velocity conditions are imposed at the inlet, while the pressure conditions are imposed at the outlets of the carotid artery. The Reynolds number at the inlet ranges from 431.15 to 605.81, indicating the flow regime.

The calculations are performed for a period of 0.8 s and a total time of 2.4 s, which covers three complete periods of the physiological flow.

µ_∞ = 0.00345 Pa.s, µ₀ = 0.0364 Pa.s, λ = 0.380, n = 1.45 for the Cross model (1);

µ_∞ = 0.00345 Pa.s, µ₀ = 0.056 Pa.s, λ = 3.313, n = 0.3568 for the Carreau model (2);

µ_c = 0.002982 Pa.s, τ_y = 0.02876 Pa.s et λ = 4.020 s for the Quemada model (3).

2.2.   Governing equations

The mathematical system for transient, incompressible and laminar blood flow in a rigid artery includes the continuity (4) and the Navier-Stokes equations (5):

For turbulent flows, the use of turbulent models is needed, and some modifications are included to the momentum equation (6):

The turbulent blood flow is modeled using the k-ϵ model. The turbulent kinetic energy k, as well as its dissipation rate ϵ, are obtained from the following transport equations:

µ_t is the turbulent viscosity and it is given as follow: µt=ρCµk2ϵ

G_k is the generation of turbulence kinetic energy due to mean velocity gradients; σ_k = 1 and σ_ϵ = 1.3 are the turbulent Prandtl numbers for k and ϵ, respectively; C_1ϵ = 1.44 and C_2ϵ = 1.92 are constants.

2.3.   Boundary conditions

In our study, several pressure and velocity profiles have been considered to solve the above system. Taking into account the Womersley approximation [14], the pulsating velocity profile V₁ (Figure 2a) [15], given by equation (9), will be used as the first condition at the inlet of the carotid bifurcation.

Where: A = 2/3 and ω=2πT .

The second velocity condition, V₂, is fitted using a Gaussian function (10) as described by Kleinstreuer [16]:

The coefficients values a_i b_i and c_i are given in the Table 2 below:

The pressure conditions have been used to model the pulsating pressure at the carotid outlets. The first condition P₁ (11) is the one used by P. Kumar Mandal et al [17],[18]. The pressure gradient was taken according to Burton [19] as follows [20],[21]:

Where A₀ is the constant amplitude of the pressure gradient, A₁ is the amplitude of the pulsatile component giving rise to the systolic and diastolic pressures.

The second condition, P₂, is the one used by Vasava et al. [22], it is given as an eighth-degree polynomial correlation (12) developed from data provided by Conlon et al. [23].

The C_i coefficients are summarized in Table 3:

The discussed conditions (V₁, V₂, P₁ and P₂) are shown in Figure 2:

3.   Resolution method

Our study employed the finite volume method [24] implemented in the Fluent solver to solve the governing equations. The pressure-velocity coupling was achieved using the SIMPLE algorithm. Spatial discretization involved using a second-order scheme for pressure and a second-order upwind scheme for momentum.

The boundary conditions mentioned earlier, specifically the P₁-V₁ and P₂-V₂ pairs, were applied as boundary conditions in the solver using User-Defined Function (UDF) files. These files defined the prescribed velocity and pressure conditions at the respective boundaries.

Since the phenomenon being investigated is transient, we adopted an implicit first-order scheme for time discretization. The calculations were conducted for a total duration of 2.4 seconds, equivalent to three periods of the pulsatile flow. A time step size of 0.004 seconds was employed during the simulations.

Convergence of the solution was determined based on the criterion that the velocity components fall below a threshold value of 10⁻⁶.

As presented below, we utilized an undisclosed formula (13) to evaluate the wall shear stress at the carotid walls.

4.   Mesh size effect and validation

The discretization of the physical domain directly impacts the computational process and the accuracy of the results. Thus, the analysis of the mesh effect is a precious step in any simulation procedure, aiming to ensure the independence of the obtained results from the chosen mesh. With this objective in mind, a thorough investigation of the mesh effect on the phenomena under study was conducted. Due to the intricate nature of the carotid structure, an unstructured mesh composed of tetrahedral elements proves to be more suitable for addressing our specific problem. By Figure 3, a mesh of 145068 elements was carefully chosen to accurately solve the blood flow within the designated domain, as indicated in Table 4.

To establish our solver's reliability, we compared our results with those obtained by K. Mamuna and K. Funazakia [25], who investigated blood flow in a rigid and elastic stenotic artery. Their study considered the flow incompressible, non-Newtonian using the Cross model, and turbulent utilizing the k-omega model. We specifically focused on comparing the diastole wall shear stress (at t = 0.5945 s) for the 55% stenotic artery. The results comparison, as depicted in Figure 4, reveals a slight disparity between the two curves, indicating a significant level of agreement between our findings and those of K. Mamuna and K. Funazakia [25].

5.   Results

The subsequent section presents the results obtained from simulating pulsed blood flow through a rigid carotid bifurcation, considering two sets of pressure-velocity boundary conditions (P₁-V₁ and P₂-V₂) and three viscosity models (Cross, Carreau, and Quemada). Three radial sections were created to visualize the evolution of velocity profiles within the studied geometry (Figure 5). To ensure a fully developed flow, the results are presented for the second heart cycle, specifically within the interval of [0.8 s, 1.6 s]. For each set of boundary conditions, pressure, velocity, and wall shear stress variations are presented at the systolic peak and early diastole moments. The corresponding timing for each of these moments is summarized in Table 5.

5.1.   Time evolutions

Figures 6, 7, 8, and 9 illustrate the temporal evolution of pressure and velocity in the CCA (common carotid artery) for all viscosity models and boundary conditions mentioned in section 2.1. From Figures 6a, 7a, 8a, and 9a, it can be observed that the pressure in the carotid bifurcation remains unaffected by the various parameters investigated, including the flow regime, viscosity model, and type of boundary conditions. However, when examining velocity variations, it becomes evident that the impact of the viscosity model is more pronounced in laminar flow than turbulent flow.

Despite the slight differences observed in the curves (Figures 6a and 7a), it is apparent that the Carreau model yields the highest velocity value (37 m/s). The Cross model exhibits a lower weight (36.5 m/s). In turbulent flow, the simulation results using the P₁-V₁ boundary condition (Figure 8b) demonstrate that the velocity differences between the Newtonian, Cross, and Carreau cases are negligible, with all curves reaching a maximum velocity of 38 m/s. However, the Quemada model predicts the highest velocity value (38.5 m/s). When considering the second boundary condition pair, P₂-V₂ (Figure 9b), it is noteworthy that the Cross model yields the highest velocity value, while the Carreau model corresponds to the lowest velocity value.

5.2.   Pressure

Figures 10 to 17 depict the variations in flow field pressure during the systolic peak and early diastole under different conditions, including flow regime (laminar and turbulent), viscosity models (Newtonian, Cross, Carreau, and Quemada), and boundary conditions (P₁-V₁, P₂-V₂). Across all the figures, a common observation is that the pressure rapidly decreases in the CCA (Common Carotid Artery), regardless of these parameters. The unique anatomy of the carotid bifurcation, particularly the branching of the CCA and the curvature of the ECA (External Carotid Artery), introduces a radial pressure gradient and gives rise to significant secondary flows (Figures 24, 25, 26, and 27), unlike typical internal flows characterized by axial pressure gradients.

Furthermore, the pressure in the ICA (Internal Carotid Artery) is higher than in the ECA. However, the branch point of the bifurcation exhibits the highest-pressure values due to the abrupt change in velocity. The pressure distribution in this region varies depending on the viscosity model and the pair of boundary conditions. Figures 11 and 13 show that the P₂-V₂ boundary condition yields the highest-pressure values.

The transition from laminar to turbulent blood flow modeling in the carotid bifurcation is essential. The intricate anatomy of this artery theoretically promotes the development of chaotic movements, particularly in the ICA. Figures 24, 25, 26, and 27 vividly demonstrate the presence of secondary solid flows in the carotid sinus.

By examining Figures 14, 15, 16, and 17, it is evident that the pressure difference between the laminar and turbulent regimes is minimal during the systolic peak. However, during early diastole, the turbulent flow exhibits higher pressure, potentially damaging endothelial cells [26],[27].

5.3.   Velocity

The Figures presented in this section depict the velocity profiles in different areas of the geometry (S1, S2, and S3) for each case studied. In section S1, the velocity profile exhibits a quasi-parabolic shape (Figures 18 and 19). Notably, the asymmetry of the carotid artery alters the position of the velocity maxima, which is shifted towards the ICA (internal carotid artery).

In the branching zone, as shown in contours (Figures 20, 21, 22, and 23), it is observed that a significant amount of fluid is dragged towards the inner wall of the ICA [10],[28]. This region is of great physiological importance and exhibits more complex physical phenomena than those observed in the CCA (common carotid artery).

Typically, in circular pipes, decreasing the flow section results in an increase in velocity values. However, in the case of carotid bifurcation, this relationship no longer holds. The velocity values in the ICA (more extensive section) are higher than those in the ECA (smaller section). This is attributed to the increase in the CCA section for a constant flow rate and the curvature of the ECA, which introduces an opposite pressure gradient and relatively lower flow in the ECA compared to the ICA. A compensatory increase in the ICA velocity occurs to maintain flow conservation, explaining the apparent differences in velocity values between the two branches.

The shape of the branching zone plays a significant role in the deformation of flow and the interference of streamlines, leading to complex blood movement in the carotid sinus. Examining the contours in Figures 24, 25, 26, and 27, as well as the curves presented in Figures 28 and 31, it is evident that the velocity profiles are aligned towards the internal wall of the ICA (internal carotid artery). This profile shape is closely associated with a secondary flow near the outer wall of the carotid sinus and the formation of a strong vortex in the positive direction within this acceleration region.

However, despite the vortex motion, it is essential to note that this region cannot be regarded as a recirculation zone since the fluid is not trapped within it. The contours also demonstrate that velocity variations along this branch are not dependent on the flow cross-section. It can be observed that the velocity increases from the carotid sinus towards the exit of the ICA, with values higher than those surveyed in the ECA (external carotid artery). This acceleration near the wall results from locally large pressure gradients induced by the effects of curvature.

Considering the ECA (external carotid artery), the curvature of the vessel contributes to generating a radial pressure gradient, as observed in Figures 10–17. This gradient leads to a secondary flow and the emergence of two counter-rotating vortices directed toward the center of the artery, as depicted in Figures 35, 36, 38, and 39.

Pathologically, these phenomena observed in the carotid sinus are closely associated with developing atheromatous plaques, particularly in atherosclerosis's early and moderate stages. It is worth noting that the ECA is typically affected by lipid deposits, primarily in the advanced stages of the disease.

The dimensions of the physical domain highly influence the rheological behavior of blood in the carotid artery. In the CCA (common carotid artery), the distinction between Newtonian and non-Newtonian behavior of blood is minimal. As depicted in Figures 18 and 19, the curves for the Newtonian case and the non-Newtonian instances described by the Cross, Carreau, and Quemada models overlap, and the velocity values are nearly identical at the selected presentation times of systolic peak and early diastole.

However, in the ECA (external carotid artery) (Figures 34 and 37), a noticeable difference is observed compared to the ICA (internal carotid artery) (Figures 28 and 31). The diameter of the vessel directly impacts the rheological behavior of blood. Moreover, the figures indicate that the non-Newtonian behavior of blood becomes more prominent at the P₂-V₂ systolic peak. Among the non-Newtonian models, the Carreau model exhibits the highest velocity values, followed by the Quemada and Cross models. The Newtonian case demonstrates the lowest velocity values.

The following figures depict the velocity variations obtained from turbulent blood flow modeling using the k-epsilon (k-ϵ) model. In the CCA and ICA (Figures 40, 41, 50, and 53), the viscosity model and type of boundary conditions have a negligible effect on the velocity curves. However, a slight difference can be observed in the velocity curves of the ECA for the P₂-V₂ boundary condition pair (Figure 59). The velocity obtained using the Carreau model is slightly higher than that of the Quemada model, followed by the Cross model. The Cross model yields the same velocity as the Newtonian approach. This difference in velocity is more noticeable in laminar flow than in turbulent flow.

Regarding axial velocity, the contours presented in Figures (42, 43, 44, and 45) appear to have similar distributions, with only the velocity values changing with the presentation moments and the pairs of boundary conditions (P₁-V₁ and P₂-V₂).

5.4.   Wall shear stress

The figures presented in this section depict the distribution of wall shear stress along the carotid artery during the systolic peak and early diastole for all conducted simulations. In general, the wall shear stress exhibits variations, reaching a maximum value of approximately 10 Pa during the systolic peak and 8 Pa during early diastole.

Due to the Poiseuille-like flow in the CCA, it is evident that the wall shear stress decreases rapidly along this branch. However, in the branching zone, an important observation can be made. A blue C-shaped band represents a stagnation zone, indicating a significantly low wall shear stress value. Moving downstream from this band, the wall shear stress increases again, and the maximum values are found in the inner wall of the carotid sinus near the apex, creating a region of high shear stress. Wall shear stress values are higher towards the height of the bifurcation and lower towards the outer wall.

In the external carotid artery, the wall shear stress exhibits relatively high values near the apex and lower values at the outer corner.

According to studies, variations in wall shear stress within the carotid bifurcation significantly impact the development of lipid plaques and the occurrence of atherosclerosis [29]. Several experimental studies [30]–[32] have demonstrated that areas affected by atherosclerosis correspond to regions of low wall shear stress. It is commonly observed that lipid particles tend to accumulate in the C-shaped part rather than at the apex of the bifurcation. The presence of low velocity and wall shear stress in this region of the carotid artery promotes the accumulation of lipid particles [33],[34]. Moreover, some researchers have identified flow separation and carotid sinus vortices as additional factors contributing to the formation of lipid particles [35],[36].

Comparing the figures also reveals that the impact of changing between the proposed boundary conditions (P₁-V₁ and P₂-V₂) on the wall shear stress is relatively small, with comparable values observed between the two states.

When transitioning to turbulent modeling, it becomes apparent that the CCA's wall shear stress is higher than that observed in laminar flow (Figures 66, 67, 68, and 69). Regarding the viscosity model effect, the results obtained in laminar flow are consistent with turbulent flow simulations. The Newtonian behavior aligns with the Cross model (66.a, 66.b, 67.a, 67.b, 68.a, 68.b, 69.a, and 69.b), while the Carreau model exhibits similarities to the Quemada model (66.c, 66.d, 67.c, 67.d, 68.c, 68.d, 69.c, and 69.d). Moreover, slightly higher wall shear stress values are observed when using turbulent flow and the P2-V2 boundary conditions than those obtained using the P₁-V₁ pair.

6.   Conclusions

This study looked at the dynamics of blood flow in a rigid carotid bifurcation, taking into account things like the flow regime, the viscosity model, and the boundary conditions. Utilizing the Fluent solver, we observed intriguing phenomena within the carotid bifurcation and presented the results through diverse visualization methods, including curves, contours, and vectors. The findings obtained can be summarized as follows:

In the common carotid artery (CCA), the flow exhibits a Poiseuille-like profile, independent of the flow regime, viscosity model, and boundary conditions form.

The branching zone creates a depression area, where the enlargement of the CCA cross-section and the curvature of the internal carotid artery (ICA) and external carotid artery (ECA) significantly impact velocity and wall shear stress. This leads to the formation of a C-shaped stagnation zone with low velocity and wall shear stress values, accompanied by secondary flow development in the carotid sinus. These currents are most prominent during early diastole under P₂-V₂ boundary conditions.

In the ICA, the flow is directed towards the inner wall with a higher velocity than in the ECA. The differences between the viscosity models are more noticeable in laminar flow than in turbulent flow, particularly when considering the actual boundary conditions of P₂-V₂.

The non-Newtonian behavior of blood in the ECA is more apparent during laminar flow when utilizing the P₂-V₂ boundary conditions.

It's imperative to acknowledge an inherent limitation in our analysis. We admit that we didn't factor in the variations in distal vascular resistance between the internal and external carotid arteries, which could significantly impact the formation of observed flow patterns and fluid behavior at the bifurcation point. These fluctuations in vascular resistance might introduce complexities that haven't been examined in this study. While we've considered this limitation when interpreting our findings, the fact that we haven't explicitly analyzed this variable is a limited aspect of our investigation. Subsequent research could benefit from a more comprehensive exploration of this factor to fully understand the mechanisms influencing flow patterns at the bifurcation.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.