Blood flow CFD simulation on a cerebral artery of a stroke patient

The purpose of this paper is to conduct a numerical simulation of the stroke patient's cerebral arteries and investigate the flow parameters due to the presence of stenosis. The computational fluid dynamics (CFD) simulations are based on simplified and realistic cerebral artery models. The seven simplified models (benchmarks) include straight cylindrical vessels with idealized stenosis with variable d/D (0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1). The realistic model of the cerebral artery is based on magnetic resonance imaging (MRI) for patient-specific cerebral arteries. The simulation for the realistic model of the cerebral artery is performed at boundary conditions measured by ultrasonography of the input and the output flow profiles (velocity and pressure). The obtained CFD results of the benchmarks are validated with actual data from the literature. Furthermore, a previous vascular contraction is assumed to be exist and the effect of this contraction area ratio on the blood flow regime is discussed and highlighted. Furthermore, CFD results show that a certain vascular contraction area critically affects the blood flow which shows increasing the wall shear stress WSS at the stenosis site. An increase in the blood velocity and vortex appears after the contraction zone, this lead to vessel occlusion and strokes. The pressure drop across the arterial contraction is reduced when the area ratio d/D is increased. In some cases, the vortex can prevent blood flow from crossing, this leads to vessel occlusion especially at low d/D The WSS near the contraction area is high. Increasing the WSS can cause embolism that leads to lead to vessel occlusion. The pressure drop across the arterial contraction is reduced when the area ratio d/D is increased. In some cases, the vortex can prevent blood flow from crossing, this leads to vessel occlusion especially at low d/D The WSS near the contraction area is high. Increasing the WSS can cause embolism that leads to lead to vessel occlusion.


Introduction
Hemodynamic analysis of the patient's arteries is important to understand and diagnose the diseases of the vascular system (saccular aneurysm, fusiform aneurysm, angiosclerosis, and angiostenosis). There are many studies have been conducted to investigate some of the blood flow parameters such as flow recirculation, wall shear stress, and wall pressure distribution. These parameters are strongly related to cerebrovascular diseases [1][2][3]. The computational fluid dynamics (CFD) technics for the hemodynamics analysis are done using two geometrybased approaches, they are benchmarks (simple models), and realistic models (based on real parts of the vascular system) [4][5][6][7][8]. Klepaczko et al. (2014) conduct a simulation for the blood vessels (angiography imaging) of Magnetic Resonance. The study cases used one for known modeling of MRA imaging which was used for validation and ensuring the adapted methodology. A blood flow simulation MR signal computation kernel is used in the second model. The model presented has been validated experimentally using the flow phantoms. The findings support the development of the technique for quantifiable validation of MR angiography image processing algorithms [9]. Perinajová et al. (2021) study the blood flow turbulent and wall shear stress (WSS) in aortic stenosis using simulations. The study is conducted in the beginning on a flow phantom (Vessel with 180 bending with narrowing) with unsteady flow (pulsating flow). MRI using a 4D-flow MRI technique is used, besides the CFD method. The results show an excellent agreement between the experimental work done using the flow phantoms and the CFD simulations [10]. Zhang et al. (2018) perform a CFD simulation based on Digital Subtraction Angiography (DSA) for a 65-yearold patient with asymptomatic carotid contraction. The boundary condition used is simplified, and the CFD simulation was done in DSA data to calculate the Carotid arterial functional assessment (CAFA) index. The CFD study conducted shows a good correlation with the proposed idea which is the CFD simulation is a cheap, efficient, and time-saving solution for studying the blood flow disorder created by the internal carotid artery (ICA) stenosis [11]. Chen et al. (2020) conducted a study on a stroke patient's blood flow in the cerebral artery. The purpose of the study is to find a time-saving solution for the complex models of the cerebral arteries which makes the numerical simulation a good fit for pre-and/or post-surgery planning, and that's what the results showed [12].
Annually, 15 million people suffer a cerebrovascular accident (stroke) worldwide. Of these, 33.33% die, and 33.33% are left permanently disabled. Cerebrovascular accident is not common among the youth, when it does occur, the main reason is high blood pressure [13]. Nowadays, understanding cerebrovascular diseases, and finding their relationship with the blood flow parameters (flow recirculation, wall shear stress, wall pressure distribution, Vortex, etc.) is a major issue. Zouggari (2018) perform a numerical study to assess the carotid atherosclerosis severity. The purpose is to investigate plaque appearance and its effect on the WSS distribution. The investigation shows that WSS values are reduced after the zone that has a plaque and leads to the presence of the flow recirculation [14]. Xie et al. (2018) study the flow disorders (such as the flow recirculation) presented due to the stenosis in the eccentric coronary. The numerical study was done using several benchmarks (Models of idealized eccentric coronary stenosis) each one simulating one of the morphological parameters (such as diameter stenosis (DS)). Models with sharper stenosis shapes have a higher maximum shear rate and a larger recirculation zone. As DS grows, the recirculation zone length and maximum shear rate both increase significantly [15]. Chaichana et al. (2011) conduct a numerical simulation to study the effect of hemodynamics due to changes in the left coronary artery angle. The study consists of 12 models (8 benchmarks and 4 realistic). The benchmarks have different angles of the artery varied from (15°-120°). The models with wide angles were observed to have a recirculation flow at the beginning of the branches. The WSS investigation shows a low distribution at left bifurcations for the models with wide angles. The pressure contours show a reduction when the flow changed from the left main stem to the bifurcated regions [16]. Lu and Zhang (2019) investigate the effect of physiological parameters on vascular wall stress (WSS). A numerical study has been adopted to estimate the effect of four physiological parameters (including blood pressure, blood flow rate, vascular elasticity, and blood viscosity) on vascular wall stress (WSS) in the right coronary artery. The WSS is affected significantly by the blood flow rate, the blood flow rate has the main role in the WSS variation under normal physiological conditions. In contrast, the blood pressure has a minor influence on WSS [17]. Luo et al. (2019) perform a CFD simulation on a patient's artery with stenosis to study the blood flow problems to understand the hemodynamics. The simulation performed in the beginning on a simplified geometry (benchmark) to validate the methodology and the boundary condition adapted for a straight cylindrical vessel with an idealized stenosis with d/D = 0.4, while the patient-specific cerebral arteries with two inlets and 15 outlets. The second part of the simulation conducted on a cerebral artery geometry based on a real case of a patient with artery stenosis with measured boundary condition such as velocity inlet and pressure outlet profiles. The results show the intracranial stenosis appears to be causing distal pressure loss [18].
In the present paper, we focus to include the entire cerebral artery with patient-specific geometry and parameters. In addition, the effect of a certain previous stroke with different contraction area is included. The remainder of this paper is organized as follows: In Sect. 2, the computational domain and the dimensions of the geometry used for the simulation are discussed in detail. Furthermore, the problem description and the system of incompressible Navier-Stokes equations with a fully implicit finite element discretization are presented. Moreover, the carreau model is used to characterize the behavior of the blood which has shear thinning properties and fluids with shear thickening properties [14]. In Sect. 3, the CFD results are first validated with the results found in literature review [18]. In addition, a deep study for the blood flow and hemodynamic parameters are conducted on benchmarks, where a straight pipe with certain vascular area contraction is represented to simulate the straight cylindrical. The solution is upgraded to consider an actual vascular stroke case in patient-specific cerebral arteries. In Sect. 4, some concluding remarks are given.

Computational model
The blood flow is simulated in a pipe of 3.5 mm diameter and length of 61.6 mm with a vascular contraction ratio range from 0 to 60% in the first case (benchmark). This contraction is assumed to be a previous in-born medical issue with the patient. Furthermore, different inlet mass flow rates are trained to figure out the blood flow regimes and recirculation zone after this previous medical stroke contraction. In the second case, a real life situation for patient-specific cerebral arteries [12] is considered. The computations for both cases are performed with similar CFD model adjustment to present the blood hemodynamics in the considered geometry.

Governing equations
Steady-state three-dimensional CFD computations were performed for the considered geometry shown in Fig. 1. Differential equations governing the blood flow are given as in [19,20]. The blood motion was governed by Navier-Stokes equations in the laminar flow condition as: where u, v, w are the blood flow velocities in x, y, z coordinates, respectively. The blood flow pressure defined as P while is the blood density. It is assumed that blood can be described as incompressible Newtonian fluid with viscosity model (Carreau model) and density of 1063 kg/m 3 . The wall of the blood vessels supposed to be rigid without slipping as shown in Tables 1 and 2 for the blood properties and Carreau model [14]. It is noted that Carreau fluid model treats the flowing fluid (blood in the present case)   is a type of generalized Newtonian fluid. That is mean, at low shear rate Carreau fluids behaves as a Newtonian fluid and at high shear rate as power law fluid. So, it is recommended for blood artery simulation and accurate prediction [21]. It is noted that shear-thinning and wall shear stress with resistive impedance and their growth are successfully captured by using Carreau fluid model. These assumptions account for the rigid wall model prescribed in the present simulation. The carreau flow model equations can be found as follows [20].
where H(t) is the temperature dependence, known as the Arrhenius law.
where is the ratio of the activation energy to the thermodynamic constant and T is a reference temperature for which H(t) = 1 . The T 0 is the temperature shift, is set to 0 by default, and corresponds to the lowest temperature that is thermodynamically acceptable. The time constant is defined by .

Numerical scheme
In the present simulation, the governing equations are solved in the three-dimension domain using the commercial package Ansys Fluent 19.2 [20]. The flow is assumed to be laminar as the Reynolds number is less than 2000 in all flow cases. Also, a second-order upwind scheme are used for all conservation equations to save computational time. The coupled scheme is used for pressure-velocity coupling and the pressure is calculated with the second-order scheme. The coupled algorithm solves the momentum and pressure-based continuity equations together which makes the solution robust and more efficient over the segregated solution scheme. To accomplish robust convergence, various pseudo time values are tested for the fluid and solid zones of the computational domain. The selected pseudo time scale factors for fluid and solid zones are 0.7 and 1, respectively. The pseudo-transient fashion is recommended for the complex blood flow. A grid independence study is performed to grantee the computational results and the generated mesh is shown in Fig. 2. The total number of the used tetrahedral mesh is 897,865 for the straight cylindrical vessel with an idealized stenosis (case 1) and 1,292,839 for patient-specific cerebral arteries (case 2), the number and the size of the grids are previously recommended by a previous experimental/CFD study [18]. The mesh sensitivity study was imported from the previous validated study of the model results confidence [18]. The mesh is refined near the tube wall and y + is approximately equal to 1. Moreover, standard wall treatment is used for a smooth transition between the boundary layers. Blood flow simulation is highly influenced by the wall effect. Therefore, y + has been tried to be close to 1 even if the flow is considered to be laminar.
Boundary conditions specify the flow on the physical domain boundaries. The velocity inlet boundary condition is applied at the vascular inlet. The pressure outlet boundary condition is specified as a constant value and is equal to the pressure set at the vascular outlet. The isolated wall of the duct has been set as hydraulically smooth walls with non-slip boundary conditions. It is noted that the artery wall is considered as a rigid and the blood/wall interaction is not considered in the present simulation but will be accounted in the future study. The hybrid initialization is applied for the present

Proposed CFD model validation
The obtained CFD results are first validated with the literature data reported in [18]. It is clearly seen form Fig. 3 for the velocity and pressure contours that our present CFD model is able to capture the reported contours in [18]. In Fig. 3a and b, there are a small difference in the length of the velocity zone (downstream the contraction) comparing to the reference paper which appears longer than the results presented. The variation in those results is due to using different number of elements. In the present work the number of elements used is 1,292,839 and 207,467 elements for the reference. Furthermore, in the present work the simulation are done using the commercial package Ansys Fluent 19.2, and for the reference a specific algorithm is used to perform the simulation. Moreover, Table 3 presents the pressure drop comparison for straight cylindrical vessel with an idealized stenosis between the present CFD and results reported in [18]. The error difference in pressure drops between the present work and Luo, et al. [18], the results is estimated to be 1.56% which is within the acceptable range. This validation grantee the present CFD results for the straight cylindrical vessel with an idealized stenosis.

Case (1), results for the straight cylindrical vessel with an idealized stenosis
This section introduces the CFD results for the benchmark which is a straight pipe with area contraction as discussed previously. The effect of contraction area ratio on the artery inlet/outlet pressure ratio is shown in Fig. 4. Rising the area ratio d/D decreases the pressure drop across the artery contraction which means less flow energy loss   through contraction passing. This is due to blood eddies and vortex which is highly noticed at the lower contraction area. This blood vortex is due to the blood flow contraction/enlargement found in this case. This contraction vortex extracts the blood flow energy and, in some cases, prevent blood flow crossing. Lower contraction area means blood flow accumulation, which gives a sense of surgery, and a medical treatment should be introduced for the patient. This phenomenon is easily visualized in Figs. 7 and 8 for the blood flow streamlines. Moreover, increasing the blood velocity and vortex occurs after the contraction area which alter the hemodynamics of the blood flow as depicted in Figs. 5, 7, and 8, respectively. The presence of the vortex (flow recirculation) leads to flow asymmetry, which appears in the low value of d/D such as (0.4, 0.5, and 0.6) as shown in Figs. 5a-c and 7a-c [22]. The difference between inlet and outlet pressure across the domain is calculated in Fig. 4 for the straight cylindrical vessel with an idealized stenosis with variable d/D, the pressure drop produced by the cross-section reduction could be shown as the pressure contour in Fig. 6. It is obvious that decreasing the vascular contraction area ratio rises the pressure losses across it which implies more blood flow energy to overcome these losses. Figure 6g, shows that the flow is fully developed as appeared on the pressure evolution along the model centerline and the entrance length is sufficient for obtaining a developed flow. Furthermore, at d/D = 0 for the first case, the pressure difference between model outlet and inlet is calculated and compared to the Poiseuille solution and the error is found to be 2.1% which is computationally accepted. Moreover, this contraction can result in increasing the blood velocity and vortex just after it which alter the hemodynamics of the blood flow as depicted in Figs. 5, 7 and 8, respectively. Figure 9 shows that the wall shear stress (WSS) highly relies on the geometry of the vascular vessel. By decreasing the area of the stenosis (decrease the d/D value) we increase the value of the WSS at the stenosis zone. After the stenosis zone, the value of the WSS starts to reduce with a decrease in the d/D value. The decreasing of the WSS value after the zone of contraction leads presents the flow recirculation, especially at d/D ≤ 0.4, as shown in Figs. 8 and 9h. The presence of the vortex (flow recirculation) from increased turbulence in the fluid flow leads to flow asymmetry, which appears in the low value of d/D such as (0.4,

Case 2, Results for patient-specific cerebral arteries
A real Magnetic resonance imaging MRI for a full-size cerebral artery of patient-specific cerebral arteries which has a stenosis in the left middle part of the artery. The cross section of the case in Fig. 10 shows that the contraction is about 55%. The patient-specific cerebral arteries are constructed using Sim-Vascular software package shown in Fig. 11a. The computational domain has 2 inlets and 35 outlets. For this computation, we used a finite elements mesh with 1,292,839 elements and the average mesh size is around 1.89 × 10 -11 mm, as shown in the zoomed in view at Fig. 11b. Furthermore, the computed and measured velocity profiles on the artery inlet and pressure profile on the artery outlet are demonstrated in Figs. 12 and 13, respectively. In Figs. 12a and 13a, the measured velocity profiles on the artery inlet and pressure profile using ultrasonography are presented at 4 cycles (at 4 s). In Figs. 12b and 13b, the computed velocity profiles on the artery inlet and pressure profile for one cycle are presented in the red lines using polynomials equations extracted from the curve fitting, as shown in Eqs. 8 and 9. These polynomials equations are added to the Ansys fluent as User-Defined Function (UDF) to define the boundary conditions of the inlet velocity and the outlet pressure.
where the time of one cycle t vary from 0 to 1 s. These values have been used for the CFD computations to best capture the blood flow dynamics in the artery. In Fig. 14, the velocity, streamline and the pressure distribution contours are presented. The velocity streamline showed with values ranges from zero m∕s to 0.7 m∕s . The zoomed in views ensure the expected hypotheses that made based on the case 1 (the benchmark) results in Figs. 7 and 8, which the blood eddies and vortex appeared in due to lower contraction area which could lead to prevent blood flow crossing or change the blood flow symmetric distribution next to the wall. In Fig. 15, The velocity observed relatively high in the stenosis areas and the velocity in the middle of the artery is higher than the area that is near the artery wall and the distribution is not symmetric as shown in the right zoomed in view (artery Fig. 11 a The computational domain, and b a magnified view of the computational mesh used for patient-specific cerebral arteries with two inlets and 35 outlets Fig. 12 The velocity profiles imposed on the artery inlets with stenosis) compared to the section of normal artery (artery without stenosis) which show slight unsymmetrical velocity distribution, this velocity profile is created due to the blood inertia which generated higher wall stresses on the outer wall and less wall stresses on the inner wall. In Fig. 14b, the pressure is appeared with value varied from 74 into 80 mmHg, and as we expected from the simulation of the benchmark the pressure in the stenosis areas the pressure losses increase between the inlet and the outlet of the contraction zone. Figure 16 shows that the wall shear stress WSS near the stenosis is higher than other areas. In addition to a decrease of the WSS at the after-stenosis zone is observed with patient cerebral arteries. Increasing the WSS at the stenosis area in case the produced by plaque form could produce the break for that form and produce embolism which could lead to vessel occlusion and strokes.
The obtained results from simulations have an effect on the medical consultation. As result, we believe that such simulations help the medical team to have a deep

Conclusion
In this study, we have introduced a systematic imagebased CFD method to simulate patient cerebral arteries blood flow using two types of geometries the straight cylindrical vessel geometry and patient-specific cerebral arteries geometry that based on MRI of the patient. The pressure drop across the arterial contraction is reduced when the area ratio d∕D is increased, resulting in less flow energy loss during contraction passing. The value of pressure before the stenosis increased due to flow velocity decreasing. The stenosis of the cerebral arteries may result in an increase in blood velocity and vortex after the contraction zone, where the blood eddies and vortex can be seen for the lower contraction area ( d∕D ≤ 0.4). In some cases, the vortex could prevent blood flow from crossing which lead to vessel occlusion. The velocity in the stenosis areas is relatively high, and the velocity in the middle of the artery is higher than that at the area along the artery wall, and the distribution is not symmetric, as shown in the artery cross section. The wall shear stress WSS near the contraction area is high. Increasing the WSS at the stenosis site may cause the form to break and cause embolism, this can lead to vessel occlusion and strokes. The obtained results from the present work are keywords for our current and future calculations for estimating the two-phase flow with artery stroke. All the present and future findings will provide a quick alert and general view for required rapid medical treatment. This will save time and, in the meantime, save patient life specially for surgery planning such as (presurgery and post-surgery simulation).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.