Biomedical Simulation of Non-Newtonian Fluid Dynamics in Cardiovascular Systems: A Finite Volume Method Approach to Pulsatile Flow and Atherosclerosis Analysis

In the realm of fluid dynamics, the study of non-Newtonian fluids has garnered significant attention due to their complex flow behavior, which deviates from the simple linear relationship between shear stress and shear rate observed in Newtonian fluids. Non-Newtonian fluids, such as blood, exhibit a range of behaviors including shear-thinning, shear-thickening, and viscoelasticity, making their analysis crucial for understanding various biological and industrial processes. The complexity of non-Newtonian fluids lies in their viscosity, which is not a constant but rather a function of the shear rate. This dependency poses challenges for accurately predicting flow behavior, especially in intricate geometries like those found in the human cardiovascular system. One of the prevalent models used to describe the viscosity of non-Newtonian fluids is the Carreau-Yasuda model, which effectively captures the shear-thinning nature of fluids like blood. The Carreau-Yasuda model provides a more accurate representation of blood viscosity by incorporating parameters that adjust the fluid's behavior across different shear rates, thus enabling a more realistic simulation of blood flow in physiological conditions [1].

The effect of viscosity on blood flow is profound, influencing both macroscopic and microscopic dynamics within the cardiovascular system. Blood, being a shear-thinning fluid, experiences a decrease in viscosity with increasing shear rates, which is particularly relevant in the context of pulsatile flow induced by the heartbeat [2, 3]. The pulsatile nature of blood flow means that the velocity and shear rates within arteries are continuously fluctuating, leading to complex temporal and spatial variations in viscosity. These variations play a critical role in determining the hemodynamic forces exerted on the arterial walls, which are crucial for understanding the progression of diseases such as atherosclerosis. Atherosclerosis, characterized by the buildup of plaques within the arterial walls, leads to geometric constrictions that further complicate the flow dynamics. The interplay between the non-Newtonian viscosity of blood and the altered geometry due to atherosclerosis necessitates sophisticated modeling approaches to accurately capture the resultant flow patterns and stress distributions. The Carreau-Yasuda model is particularly advantageous in this regard, as it provides a robust framework for modeling the shear-thinning behavior of blood. The model is defined by the Eq. (1).

where, $\mu$ is the viscosity $\dot{\gamma}$ is the shear rate, $\mu_0$, $\mu_{\infty}$ are the zero-shear and infinite-shear viscosities respectively, $\lambda$ is a time constant is a parameter that defines the transition region between the Newtonian plateau and power-law region and is the power-law index. This model effectively bridges the low and high shear rate viscosities, capturing the essential characteristics of blood flow under various physiological conditions. By incorporating this model into numerical simulations, researchers can gain deeper insights into the flow behavior in health and disease [4-6].

To simulate the complex flow of non-Newtonian fluids in cardiovascular systems, the Finite Volume Method (FVM) is a powerful numerical technique. FVM is widely used in computational fluid dynamics (CFD) due to its flexibility in handling complex geometries and conservation properties. Unlike other numerical methods, FVM discretizes the domain into small control volumes and applies conservation laws (mass, momentum, and energy) to each control volume [7-9]. This approach ensures that the fluxes entering and leaving a control volume are balanced, which is particularly advantageous for simulating flow in irregular geometries such as those found in atherosclerotic arteries. The ability of FVM to handle complex boundaries and ensure conservation makes it an ideal choice for simulating the pulsatile and shear-dependent flow of blood in the cardiovascular system. By leveraging FVM, researchers can discretize the arterial geometry affected by atherosclerosis and accurately capture the local variations in flow velocity, pressure, and shear stress [10]. The benefits of using the Finite Volume Method to simulate non-Newtonian fluid dynamics in cardiovascular systems extend beyond accurate flow predictions. Such simulations provide valuable insights into the hemodynamic factors that contribute to the initiation and progression of cardiovascular diseases. For instance, understanding the distribution of wall shear stress and its temporal variations can help identify regions susceptible to atherosclerotic plaque formation [11]. Furthermore, these simulations can aid in the design of medical interventions, such as stents, by predicting the impact of different geometries and materials on blood flow [12, 13]. The integration of physiological and biomathematical knowledge through advanced numerical simulations also enhances our understanding of the fundamental mechanisms driving blood flow and its interaction with arterial walls. This comprehensive approach bridges the gap between theoretical models and clinical applications, enabling the development of more effective diagnostic and therapeutic strategies for cardiovascular diseases [14, 15].

Figure 1 illustrates a conceptual model of an artery affected by atherosclerosis. The diagram depicts the lumen (inner channel) radius, representing the open space through which blood flows, and the outer wall radius, representing the artery's external boundary. The constricted region, spanning 3 mm in length and narrowing to a minimum lumen radius of 0.5 mm, demonstrates the reduced blood flow caused by plaque buildup. The flow pathways (red arrows and dots) highlight the impact of the narrowing on blood dynamics, including velocity and shear stress distribution. This representation is crucial for understanding how arterial constrictions affect blood flow, potentially leading to conditions such as ischemia or plaque rupture. The study of non-Newtonian fluid dynamics, particularly in the context of blood flow in atherosclerotic arteries, is essential for advancing our understanding of cardiovascular health and disease. The Carreau-Yasuda model provides a robust framework for capturing the shear-thinning behavior of blood, while the Finite Volume Method offers a powerful tool for simulating the complex geometries and pulsatile flow conditions characteristic of the cardiovascular system [17, 18]. By combining these advanced modeling techniques, researchers can achieve accurate and detailed simulations that not only enhance our understanding of hemodynamics but also inform the development of better clinical interventions. The integration of physiological and biomathematical knowledge through these simulations holds the potential to significantly improve cardiovascular health outcomes, making this area of research both scientifically intriguing and clinically relevant [19-22].

The novelty of this research lies in its comprehensive approach to simulating non-Newtonian fluid dynamics within cardiovascular systems, specifically addressing the complex interplay between pulsatile blood flow and atherosclerotic artery constrictions using advanced numerical methods. By integrating the Carreau-Yasuda model to accurately capture the shear-thinning behavior of blood and employing the Finite Volume Method (FVM) to discretize intricate arterial geometries, this study provides unprecedented insights into the hemodynamic environment under realistic physiological conditions. Unlike previous studies that often simplified either the fluid properties or the geometrical complexities, this research combines both aspects, allowing for a more precise analysis of the temporal and spatial variations in shear stress and pressure due to pulsatile flow. The focus on the impact of pulsatility in the presence of atherosclerosis offers a deeper understanding of how cyclic variations in flow dynamics influence plaque progression and endothelial function, which is critical for developing more effective diagnostic tools and therapeutic strategies.

This holistic approach not only enhances the theoretical understanding of cardiovascular hemodynamics but also has significant implications for clinical applications, making it a pioneering effort in the field of cardiovascular fluid dynamics.

The organization of this paper is structured to provide a systematic and comprehensive exploration of non-Newtonian fluid dynamics in cardiovascular systems, focusing on the interaction between pulsatile blood flow and atherosclerotic constrictions [23, 24]. The paper begins with an Introduction section, offering background information on non-Newtonian fluids, the significance of viscosity variations, the Carreau-Yasuda model, and the relevance of these topics to cardiovascular health. The Literature Review follows, summarizing previous research efforts in the field, highlighting gaps that this study aims to address, and establishing the novelty of the current research. The Methodology section details the numerical techniques employed, including the Finite Volume Method (FVM) for discretizing the complex arterial geometry and implementing pulsatile boundary conditions, as well as the integration of the Carreau-Yasuda model to represent blood viscosity [25, 26]. The Modeling and Simulation section presents the development of the 2D atherosclerotic artery model, the application of the FVM, and the simulation setup to capture pulsatile flow dynamics. The Results and Discussion section showcases the simulation outcomes, analyzing the spatial and temporal variations in shear stress and pressure, and discussing their implications for atherosclerosis progression and endothelial function [27, 28]. This is followed by the Validation and Comparison section, where the results are compared with experimental data and previous studies to establish the accuracy and reliability of the simulations. The Implications for Clinical Practice section explores how the findings can inform the development of diagnostic tools and therapeutic strategies. This structured approach ensures a thorough examination of the research questions and facilitates a clear understanding of the study's contributions to the field of cardiovascular fluid dynamics.

Continuity equation (Mass Conservation) is:

For incompressible flow, this simplifies to:

For an incompressible, non-Newtonian fluid, the Navier-Stokes equations in two dimensions (2D) are [29]:

where, $\mathrm{u}=(u, v)$ is the velocity vector with and being the velocity components in the and directions, respectively.

For non-Newtonian fluids, the stress tensor is given by [30, 31]:

where, $\mu(\dot{\gamma})$ is the dynamic viscosity, which is a function of the shear rate $\dot{\gamma}$, $-\mathrm{D}=\frac{1}{2}\left(\nabla u+(\nabla u)^T\right)$ is the rate-of-deformation tensor.

In the presence of atherosclerosis, the geometry of the artery is no longer a simple straight channel. Instead, there is a constriction that affects the flow characteristics. This constriction can be modeled as a reduction in the cross-sectional area of the artery in the constricted region. This impacts the continuity and momentum equations by introducing spatial variations in the flow properties. Pulsatile flow, such as the flow induced by the heartbeat, introduces time-dependent boundary conditions. This affects the momentum equations by adding a time-varying component to the velocity field at the inlet. The continuity equation for incompressible flow remains the same but must be satisfied in the modified geometry.

The Navier-Stokes equations in the presence of atherosclerosis and pulsatile flow are given by x and y directions respectively:

Eq. (7) and Eq. (8) model blood flow’s complex behavior in arteries affected by atherosclerosis, where constrictions and pulsatile flow create time-varying and spatially-dependent velocity and pressure fields. This helps to understand shear stress distribution and hemodynamics in diseased arteries.

The shear rate in the presence of atherosclerosis must be calculated considering the modified velocity gradients due to the constriction.

The viscosity $\mu(\dot{\gamma})$ is then given by the power-law model.

The pulsatile nature of the blood flow can be modeled by a time-varying inlet velocity [32, 33].

where, $U_0$ is the mean inlet velocity, $\alpha$ is the amplitude of the pulsation, f is the frequency of the pulsation, and t is time. This time-varying inlet velocity affects the boundary conditions of the Navier-Stokes equations, making the problem time-dependent [34, 35].

There are to criteria to ensure that the fluid flows in a stable condition [36, 37].

where, $C_{\max }$ the number of counts that are usually less than 1 for explicit scheme stability.

where, $D_{\max }$ the maximum limit to ensure stability, usually less than 0.5.

The algorithm for simulating non-Newtonian fluid flow in arteries was somewhat basic, with a focus on simple geometry definitions, straightforward boundary conditions, and a limited explanation of fluid properties [38]. It lacked considerations for complex artery geometries, like realistic atherosclerosis modeling, and did not fully account for the pulsatile nature of blood flow. The approach to grid creation was also standard, without utilizing advanced techniques like adaptive mesh refinement. Additionally, there was minimal detail on post-processing or visualization, which are critical for interpreting the results of a simulation in practical medical or engineering contexts.

Based on the algorithm, it became more comprehensive and capable of handling realistic simulations of blood flow in arteries affected by atherosclerosis. The enhancements included more sophisticated geometry definitions, such as using medical imaging data to model arterial constriction accurately. The boundary conditions were refined to include pulsatile inlet flow and pressure outlet conditions, which better reflect the natural behavior of blood flow (Table 1).

The grid generation process was made more efficient by incorporating techniques like adaptive mesh refinement, which allows for greater resolution in critical areas like the constriction. The fluid properties were also enhanced by introducing a shear-thinning model with regularization, ensuring stable and realistic simulations of non-Newtonian blood flow. Finally, the post-processing step was expanded to include detailed analysis of wall shear stress, pressure drops, and flow visualization, making the algorithm more suitable for medical and engineering applications in understanding arterial blockages and blood dynamics.

The Navier-Stokes equation modified for the blood flow (non-Newtonian fluid) can be written as:

Incorporating the displacement function $\varepsilon$ involves introducing a term that characterizes the alteration in position or motion of the vessel wall, which in turn impacts the flow of blood. Let $\varepsilon$ be a position-dependent and time-dependent movement function. To make this adjustment, simply incorporate the word into the speed component $u$. Let's consider the new velocity component as $\mathrm{u}+\varepsilon$. The Navier-Stokes equation modified with the function to be:

Expanding this equation further will yields as follows:

Since $\varepsilon$ is a small transition function, we can ignore the term small square $\varepsilon . \nabla \varepsilon$. So, the equation can be simplified to:

The impact of the addition of functions on the results of simulation and analysis of blood flow in the cardiovascular system can have several important aspects. The function, which denotes a displacement or deformation of the vessel wall, will result in an alteration of the blood flow velocity profile. This is because the displacement of the wall will impact the boundary condition of the blood flow. When the wall of the vessels shifts or displaces, it also affects the blood flow near the wall, causing a change in the distribution of speed inside the vascular intersection. The blood arteries might experience pressure variations because to the deformation of the wall caused by $\varepsilon$. The pressure will intensify in the region where the vessel wall constricts and diminishes in the regions where the wall widens. Understanding the hemodynamic response to changes in the elasticity of blood vessel walls requires a thorough comprehension of this pressure fluctuation. The inclusion of will impact the allocation of tension on the walls of the blood vessels. Understanding the relationship between wall tension and plaque development is crucial in studying atherosclerosis, as regions experiencing elevated wall tension are more prone to plaque formation. This model has the potential to offer an enhanced understanding of regions with a high susceptibility to atherosclerotic plaque development by providing a more precise depiction of stress distribution.

The momentum equation in the x direction is:

Discrimination using an explicit scheme:

From Eq. (17), we evaluate to solve the as follows:

The momentum equation in the y direction is:

Discrete using an explicit scheme:

From (19) we evaluate to solve the as follows:

The continuity equation is:

Discrimination using an explicit scheme:

To ensure the stability of the explicit scheme at Eq. (12) and Eq. (13) we know that the parameters.

• 1. Courant-Friedrichs-Lewy Condition

where, the we know that:

So that with a time step less than or equal to 0.01 seconds, we can run a simulation with guaranteed numerical stability.

We used computational fluid dynamics (CFD) to simulate blood flow through a constricted artery, aiming to understand the flow behavior in both normal and pathological conditions. By utilizing CFD, we were able to model the complex dynamics of blood flow, including the velocity distribution, wall shear stress, and the presence of recirculation zones that typically occur downstream of a constriction due to atherosclerosis. The artery was modeled as a 2D domain with a length of 3 mm and an outer radius of 1 mm. A constriction representing stenosis was applied in the central region, reducing the artery’s radius by 50% over a 1 mm segment. The inlet and outlet regions on either side of the constriction maintained the full radius. This simplified geometry provided a controlled environment to study how blood flow is altered by stenosis. The blood flow was simulated using the incompressible Navier-Stokes equations, which are fundamental to fluid dynamics and capture the conservation of mass and momentum:

Blood was modeled with a density of 1060 kg/m³ and a dynamic viscosity that varies with shear rate. , This configuration accurately represents the shear-thinning nature of blood under different flow conditions.

The meshing shown in Figure 2 geometry is crucial for studying blood flow in constricted arteries, as it allows the simulation to focus on the dynamics of how blood accelerates through the constriction and how pressure and wall shear stress are affected. This kind of meshing is often used in computational simulations to ensure accurate results near-critical regions, such as around the constriction where flow separation and recirculation may occur.

Figure 3 compares the velocity magnitude and shear stress distributions for Newtonian and Non-Newtonian blood flow in an artery. In both velocity plots (top row), the non-Newtonian fluid maintains higher velocities in the core region due to the shear-thinning behavior, while the Newtonian fluid exhibits a smoother decrease in velocity along the artery. The shear stress plots (bottom row) reveal significantly higher wall shear stress for the non-Newtonian fluid, particularly downstream, where the fluid experiences reduced viscosity in regions of high shear. This indicates that modeling blood as a non-Newtonian fluid better captures the dynamics near the artery walls, which is crucial for understanding the mechanical stresses that may contribute to cardiovascular conditions such as plaque rupture. Before this comparison, it is important to understand that blood behaves more like a non-Newtonian fluid under physiological conditions. The assumption of Newtonian behavior is a simplification often used in theoretical studies but may not accurately capture the behavior of blood, especially in constricted arteries.

Figure 4 demonstrates that non-Newtonian fluid models better capture the complexity of blood flow, especially in regions of high shear near artery walls. The higher shear stress in the non-Newtonian model suggests that neglecting this behavior may underestimate the mechanical forces acting on the artery walls, which are critical in understanding diseases like atherosclerosis. Before this model, understanding how blood flow is affected by stenosis (narrowing of the artery due to plaque) was challenging. Clinicians needed to predict how blood would behave near and downstream of such constrictions, particularly in regions where flow separation and vortices form. Computational fluid dynamics (CFD) simulations like this model enable more detailed analysis of flow behavior, wall shear stress, and areas prone to recirculation or turbulent behavior.

Figure 4 visualizes the blood flow dynamics within a constricted artery, modeling the effects of atherosclerosis. The color gradient represents the velocity magnitude of blood flow, with higher velocities (red) concentrated in the narrowed section of the artery and lower velocities (blue) occurring in less restricted areas. Streamlines indicate the flow direction and the formation of recirculation zones or vortices downstream of the constriction, where the flow becomes disrupted. The black contours highlight areas of wall shear stress, showing elevated stress near the constricted region and lower stress in recirculation zones. Figure 4 demonstrates the complex interaction between blood velocity, wall shear stress, and recirculation zones, which are critical factors in understanding plaque progression and thrombus formation in atherosclerosis. After understanding the atherosclerosis model, it becomes evident how constricted arteries (due to atherosclerosis) affect blood flow. The increase in velocity in the constricted region can lead to high wall shear stress, while recirculation zones downstream are areas where blood flow stagnates, potentially promoting plaque buildup. These models help predict regions at risk for disease progression or intervention.

Figure 5 illustrates velocity magnitude with potential turbulence in an artery with a constriction, showcasing pulsatile flow. Before the pulsatile flow enters the artery, the fluid moves steadily, following a smooth laminar pattern. However, as pulsatile flow is introduced, the velocity fluctuates over time, mimicking the natural heartbeat. This leads to variations in flow speed and direction, with higher velocities seen near the center and slower flow near the walls. After the pulsatile flow is fully developed, the flow pattern becomes more complex, especially in the constricted area, where flow separation and potential turbulence occur. The streamlines in Figure 5 indicate disrupted flow near the walls and in the wake of the constriction, suggesting the development of recirculation zones or vortices. These turbulent regions, combined with the pulsatile nature of the flow, demonstrate how blood velocity changes dynamically throughout the cardiac cycle, potentially increasing stress on artery walls.

Figure 6 illustrates key aspects of blood flow in a constricted artery simulation. Figure 6(a) shows velocity profiles at different locations (x = 0.0005 m, 0.0015 m, 0.0025 m), with the highest velocity near the beginning of the artery (x = 0.0005 m), decreasing sharply as it moves further downstream. This indicates a deceleration of flow due to the presence of the constriction. Figure 6(b) displays the wall shear stress at time step 100, showing how shear stress peaks near the entrance of the artery and gradually decreases along its length. Higher shear stress near the inlet could result in an increased risk of endothelial damage. Figure 6(c) compares the wall shear stress at different time steps (50, 100, 150), demonstrating that while the overall trends remain similar, the shear stress magnitude slightly reduces with time, likely due to changes in the flow’s development over time.

Figure 6 collectively reveals how velocity and shear stress are affected by the geometry and flow conditions, which is critical for understanding the risk factors for cardiovascular issues in constricted arteries. The simulation results highlight that regions of high wall shear stress near artery constrictions can contribute to plaque rupture, while areas of low shear stress and flow recirculation downstream of the stenosis may promote plaque formation. The non-Newtonian behavior of blood, which causes viscosity to vary with shear rate, adds complexity to the flow, increasing the risk of thrombosis in low-shear zones. Arterial geometry, such as curves or bifurcations, further influences the flow, creating more disturbed regions. Incorporating personalized CFD simulations based on patient-specific artery models could enable tailored treatments, helping clinicians predict high-risk areas for atherosclerosis progression and improve the outcomes of interventions such as stenting, medication, or bypass surgery.

The computational simulation conducted in this study provides valuable insights into the hemodynamic factors influencing the progression of atherosclerosis in arteries with stenosis. The results emphasize the critical role of wall shear stress (WSS) in determining plaque behavior, with high WSS regions near the constriction associated with endothelial damage and the potential for plaque rupture, while low WSS regions downstream promote plaque formation. These findings are consistent with previous studies showing that mechanical forces exerted on the arterial wall are pivotal in atherosclerotic plaque development and rupture risk.

The presence of recirculation zones and flow separation downstream of the stenosis was particularly notable in our simulations, where the disrupted flow led to turbulent behavior and low-shear regions. These areas create a favorable environment for lipid accumulation and plaque development, further exacerbating the narrowing of the artery. This is a critical observation, as flow disturbances of this nature are well-known precursors to plaque progression and are common in patients with severe atherosclerosis.

Incorporating non-Newtonian behavior in the simulation was another significant factor, as it allowed for a more realistic representation of blood flow dynamics. Blood viscosity decreases in high-shear regions, such as within the constricted section, leading to faster flow velocities through the stenosis. Conversely, in low-shear regions, the viscosity increases, potentially leading to thrombus formation due to stagnation and clotting risks. This highlights the importance of considering non-Newtonian properties when modeling blood flow, as assuming Newtonian behavior may underestimate the mechanical stresses acting on the artery walls, particularly in areas of stenosis.

The results also underscore the influence of arterial geometry on flow dynamics. While our model focused on a simplified constricted artery, variations in geometry, such as curved arteries or bifurcations, would likely introduce additional complexities, including more pronounced recirculation zones and localized peaks in WSS. These geometric factors, combined with flow disturbance, are key contributors to the progression of atherosclerosis and must be considered when assessing patient risk. From a clinical perspective, these simulations provide valuable information that can be used to predict high-risk areas for plaque rupture and progression. By identifying regions of high WSS and flow stagnation, clinicians could tailor treatments to target these areas more effectively. For instance, stent implantation in regions of high WSS could help stabilize vulnerable plaques, while medication aimed at reducing blood viscosity may help prevent clot formation in low-shear regions. Moreover, using personalized CFD models based on patient-specific arterial geometries could enhance diagnostic precision and treatment planning, offering a customized approach to managing cardiovascular diseases. This study demonstrates the critical role of hemodynamic forces in the progression of atherosclerosis and highlights the potential of computational simulations in enhancing clinical understanding and treatment. Future work could expand on these findings by investigating more complex arterial geometries and incorporating patient-specific data to refine predictions and optimize interventions for those at risk of severe cardiovascular events.

The computational simulation of blood flow through a constricted artery reveals the significant impact of wall shear stress, flow recirculation, and non-Newtonian fluid behavior on the progression of atherosclerosis. High wall shear stress in the stenosed region increases the risk of plaque rupture, while low shear stress and recirculation zones downstream promote plaque formation. The results underscore the importance of incorporating patient-specific arterial geometries and non-Newtonian properties of blood in simulations to accurately predict high-risk areas for atherosclerosis progression. These findings suggest that computational fluid dynamics (CFD) simulations can be a valuable tool for improving diagnosis, treatment planning, and personalized management of cardiovascular diseases. The implications of this study's findings on blood flow simulations in constricted arteries extend significantly to real-world healthcare and the management of atherosclerosis. By providing detailed insights into how wall shear stress and flow recirculation contribute to plaque formation and rupture, this research highlights the potential for using computational fluid dynamics (CFD) in clinical settings to improve early detection of high-risk regions in arteries. Through patient-specific simulations, physicians could tailor interventions such as stent placements, angioplasty, or drug treatments more effectively by targeting areas most prone to complications, thereby reducing the risk of heart attacks or strokes. Moreover, CFD models could become part of routine diagnostic protocols, allowing for personalized treatment plans that account for individual arterial geometries and flow dynamics, ultimately leading to better outcomes in managing cardiovascular diseases. In a broader public health context, this approach could help in preventative medicine, reducing the overall burden of atherosclerosis by identifying at-risk individuals earlier and applying interventions before severe symptoms manifest.