{'fraction_non_alphanumeric': 0.025415764311690434, 'fraction_numerical': 0.008785191844073061, 'mean_word_length': 4.8257557724351345, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Haemodynamic analysis using multiphase flow dynamics in tubular lesions

Background and Objective: The role of red blood cell dynamics is emphasised in certain cardiovascular diseases, and thus needs to be closely studied. A multiphase model of blood flow allows the resolution of locally varying density of red blood cells within a complex vessel geometrical domain, and haemodynamic consequences of such build up. Methods: A novel computational fluid dynamics solver for simulating multiphase flows is used for modelling blood flow using level set for a sharp interface representation. Single-phase simulations and reduced order models are used for pressure comparisons. The new solver is used for numerically studying AHA type B lesions. The impact of hematocrit and degree of stenosis on the haemodynamics of coronary arteries is investigated. Results: The comparisons with single-phase flow simulations indicate differences in pressure when considering red blood cell aggregation. Multiphase simulations provide slightly lower pressure drop for the same stenosis severity compared to the single-phase simulations. Secondary flow patterns and the interactions between the two phases leads to the red blood cell aggregation at the end of the diastole cycle, which significantly changes the red blood cell distribution, the shear stresses and velocity in tubular lesions. Conclusions: Neither pressure drop nor mean velocity are not strongly changed in the multiphase modelling, but particle buildup significantly changes which is only revealed by the multiphase approach.


Introduction
Blood is a multi-component fluid primarily consisting of dispersed red blood cells (RBCs), white blood cells (WBCs) and platelets, in plasma. Normal levels of the volume fraction of RBCs (hematocrit) range from 35% to 50% [1].
Consequently, the rheological behaviour of RBCs strongly affect blood properties which significantly change when the hematocrit level varies for instance in hemodilution-increased plasma/anemia-reduced RBC [2]. Particle aggregation strongly depends on hematocrit and can lead to large fluctuations of shear rate locally which effects blood viscosity [3]. This non-Newtonian shear-thinning rheological behaviour of blood can have a significant impact on blood flow [4].
This manifests that a single fluid approach with uniform properties for numerically investigating cardiovascular diseases such as atherosclerosis might fail to model the underpinning phenomena due to particulate motion. for the volume fraction of each phase is solved. These Eulerian approaches offer the opportunity of calculating the forces exerted on the fluids by solving for momentum for each phase allowing for studying inter-phase mass and momentum transfer. Some studies choose to couple the Eulerian-Eulerian framework with implicit methods for interface capturing, leading to a more accurate calculation of the volume fraction for each phase [5]. These are usually based on the volume of fluid (VOF) method [6]. VOF is one of the most commonly used implicit methods which combines simplicity and mass conservation by default.
Given a sufficient grid resolution, the VOF approach models the larger scales and bulk fluid-fluid interfaces or for cases where interdispersion between two fluids can be neglected [7]. The correct calculation of the interface curvature and conservation of mass for the RBC component while maintaining numerical stability in regions where volume fraction exhibits sharp changes are extremely important for accurate multiphase CFD models for numerically investigating the behaviour of RBC and plasma in cardiovascular diseases such as atherosclerosis. Increasing the accuracy of interface capturing for the larger scales is a key point for increasing overall accuracy of the multiphase analysis. For addressing these challenges here, we deploy more advanced implicit methods rather than the traditional VOF, while guaranteeing that mass is conserved.
In this paper we propose an improved numerical approach for the calculation of the volume fraction of each phase for modelling the RBC distribution. For this purpose, we employ a coupled VOF and level set method that we have recently developed and validated for different multiphase flow cases in [8] and has been proven to give a more accurate interface description compared to other state-of-the-art VOF-based approaches. Blood is modelled as a mixture and the interface of the two phases is captured using one of the most attractive features of the level set e.g a sharp interface representation. The formulation leads to a multi-scale approach that combines the advantages of both level set for the larger scales, and Eulerian-Eulerian two-fluid approach which can be used to model the characteristics of the mixture into the under-resolved (sub-grid) scales.
The presented methodology is tested for the study of partially occluded arteries due to a stenotic obstruction of eccentric shape. Although these lesions account for the majority of lesions among patient population [9,10] they are given little attention in literature from numerical works. Contrary to symmetric stenosis models, asymmetry not only offers a more realistic representation of diseased arteries, but the generated secondary flow patterns alter the flow downstream at the end of the stenosis, resulting in a fundamentally different haemodynamics compared to the idealised axisymmetric stenoses [11]. In addi-tion to eccentricity, the length of a lesion, L s is another important factor that also affects blood flow [12]. Although some studies for discrete lesions with L s ≤ 10mm -type A in the American College of Cardiology/American Heart Association (ACC/AHA) classification of coronary lesions, have been studied [13,14], there are no detailed studies for tubular coronary lesions (type B, with 10mm ≤ L s ≤ 20mm) which involve a moderate risk, higher than in discrete lesions [15]. To the best knowledge of the authors, the effects of hematocrit and stenosis on the RBCs transport and its effect on the blood flow in tubular coronary lesions (with 10mm ≤ L s ≤ 20mm, L s = lesion length) have not been previously reported elaborately. In this study, the blood flow patterns and particulate aggregation in coronary arteries with a tubular lesion are numerically studied using the developed multiphase solver.
Owing to the complexity of the multiphase flow phenomena and their dependence on the lesion geometry, apart from the eccentricity and length of a lesion, the hematocrit is also studied for investigating its role on the local haemodynamic parameters such as the velocity and the WSS on the endothelial cell membrane in the coronary artery. Although the effects of the hematocrit in the viscosity have been previously studied [16,17], many open questions still remain regarding the local haemodynamic factors, such as the wall shear stress, red blood cell buildup and velocity [18][19][20][21][22]. These important parameters which vary alongside the curved vessel and during the cardiac cycle [23,24] are the primary focus of this research here. First, the impact of hematocrit and stenosis severity on the calculated pressure drop ∆p is studied. Calculating ∆p along a lesion is extremely important since it is typically used for characterising the significance of the stenosis through the calculation of Fractional Flow Reserve (FFR). FFR is a haemodynamic index, based on invasive measurement of the pressure ratio distal and proximal to a given lesion which is employed for guiding revascularisation. Since in the particular case of tubular lesions, there are no numerical studies that evaluate the effects of the RBC levels on ∆p, the velocity and red blood cell distribution, we include some tests together with comparisons with single phase and reduced order models.
Comparisons with single-phase flow demonstrate the differences in pressure when considering the mixture. Single-phase simulations over-predict the pressure drop for the same stenosis severity compared to the multiphase simulations.
The model here generally predicted higher RBC buildup in the area of maximum curvature and the end of the stenosis compared to the upstream location. WSS region was found to vary with respect to the local RBC buildup, in the area of maximum curvature of the vessel. Secondary flow patterns and vessel geometry significantly influence the RBC buildup, primarily because of the prolonged particulate residence times. These observations could have significant implications towards better understanding of the haemodynamic parameters that might have an impact on the development of cardiovascular diseases such as atherosclerosis in tubular coronary lesions.

Materials and Methods
The open-source code OpenFOAM [25] is used for the simulations presented. The numerical implementation here has used as its basis the solver multiphaseEulerFoam solver [26] and enables simultaneous modelling of complex multiphase flows with different flow regimes (i.e. both dispersed and segregated flows). The solver uses a coupled level set (LS) and VOF interface capturing method to model the large fluid-fluid interfaces, whereas the Euler multi-fluid approach is used for modelling the dispersed flow. The two approaches are coupled by a simplified switching term, which is defined for each phase pair [26].
The solver is a combination of two multiphase solvers i.e. a coupled VOF-LS and an Eulerian-Eulerian solver with additional sub-models to model the forces involved. In addition, the solver includes models for phase change (e.g. cavitation) which is deactivated here. Although the methodology can be used to simulate multiple fluids, here we choose to present it with two fluids. This is an Eulerian-Eulerian two-fluid model which regards RBCs as suspended particles dispersed in plasma, which has been shown that it accurately describes the blood flow [4]. Following other researchers, the effect of other components such as platelets was considered here to be insignificant, since they exhibit volume fraction values which are less than 1% [18]. The multi-fluid model equations for incompressible, isothermal flow are given by a set of mass and momentum equations for each phase k, where the density, phase fraction, and velocity for phase k are given by ρ k , α, u k . The ρ k α k g is the gravitational force and the other two interfacial forces in Eq.1 are the drag force F D,k and the surface tension force F s,k . For selected phase pairs where the VOF method is used to resolve the interface, an additional compression term is added to their phase transport equation. The method here includes the level set method for capturing the interface of each pair of fluids.
The drag term F D,k is given by where the subscripts c and d denote the continuous (plasma) and dispersed phase (RBC) values respectively and K is Among the various models for calculating the drag coefficient C D available in OpenFOAM, the commonly used model of Schiller and Naumann [27] is employed here. In this model, the drag coefficient is a function of the Reynolds number Re according to where the Reynolds number is calculated based on the relative velocity of the two phases, the dispersed phase diameter and the kinematic viscosity ν according The lift force acts on the RBCs due to the developed velocity gradients in the plasma which has its own velocity field. The lift forces are calculated as A constant droplet diameter size, defined independently for each phase is used in this framework with the RBCs considered as spherical particles with a diameter of 8 µm [18].

Marker function for level set
The level set function is used here to describe the interface dynamics of the fluids. For simplicity, the case of two fluids is presented. The level set function ψ is defined to be a signed distance function calculated from the distance at the interface, Γ that separates the two fluids A given point x can have either a positive or a negative distance which defines whether it belongs in one fluid or the other. The interface that separates the two fluids is then the set of points that have a zero-distance (zero-level). A Heaviside function H(ψ, t) is used to calculate the smoothed volume fraction that depends on the values of ψ > 0. Here, the following expression is chosen [28], In order to conserve mass, the equation for H needs to be satisfied Then H(ψ, t) resembles a smoothed volume fraction α c of one fluid and 1-H(ψ, t) for the other fluid. The total mass with respect to time is then V ρ i H(ψ)dV .
The level set function according to Eq.8 is a distance function that is defined wherever an interface exists. The distance function can be advected using Eq.11 is solved first and the distance function is advected. The level set is then re-initialised for allowing ψ to remain a signed distance function. For this, the Eikonal equation |∇ψ| = 1 is solved.
The solutions of Eq. (12) have the same zero-level as ψ and are distance functions. Denoting withψ d e.g.ψ d = ψ n+1 d the new value of ψ after solving Eq. (12) an additional step for satisfying mass conservation is performed. [8] The correction denoted here withψ is added toψ d so that the following is satisfied for an appropriate ψ [29] The following reaction-diffusion equation is solved where λ is a parameter that is a function of the cell size and is generally taken The method avoids the assumption made by similar approaches for the interface position that is not located at the iso-surface contour α-0.5 as in [30][31][32][33].
Compared to other previous works, the level set function is not reconstructed only from the advected volume fraction [31,33] but neither is based on specifically designed schemes for re-initialisation [33,34] that allows for simpler and faster calculations.

Switching between segregated and dispersed flow regimes
First the transport equation for obtaining an initial approximation of the liquid volume fraction is solved. The equation includes an interface compression term u c and is given as where the velocity u c is applied normally to the interface to compress the volume fraction field and maintain a sharp interface. The α c (1 − α c ) term ensures the term is only active in the interface region. The value for the artificial interface compression velocity is given by The ∇αc |∇αc| term is the interface unit normal vector for the direction of the applied compression velocity. The coefficient C α is used for controlling the interfacial compression, and when considered to be ∈ [0, 1] it is simply a binary coefficient which switches interface sharpening on/off when it is equal to 1/0 respectively. When C α = 0 for a given pair of phase, no interface compression is imposed resulting in phase dispersion. When C α = 1, sharp interface capturing is applied and the flow is treated with VOF.
The switch between approaches is done dynamically imposing the interface sharpening in those areas of the domain that the C α disctates. This allows for a unified method which allows simulation of complex flows with any combination of regimes ranging from fully dispersed to fully segregated. The switch is done based on evaluating the magnitude of the gradient of the volume fraction (γ), assuming that when the gradient of the volume fraction becomes smaller than a defined cut-off value (γ cut−of f ) phase dispersion occurs and C α = 0. The parameter γ is then defined by [26] Thus, when γ ≥ γ cut−of f the value of C α is set to 1. The value of γ cut−of f is set to 0.4.

Results
The numerical tests here included an eccentric stenosis with a stenosis severity S D that varied from 50% to 70%. Two cases of different hematocrit levels are used, one with RBC=35% and another one with RBC=45%. The length of the stenosis L s = 14mm and the non-constricted diameter of the artery was D = 4mm, shown in Fig.1. For the inlet velocity, a pulsating unsteady blood flow with a period of T = 1s was imposed and is shown in Fig.1. The flow upstream the stenosis was assumed to be fully developed downstream with a Neumann boundary condition considered at the outlet. A zero slip velocity was employed for the RBCs and plasma at the lesion walls. The inlet profiles for velocity and RBC were considered to be uniform across the inlet cross section similar to other similar works [18].    Results are compared with the reduced-order models and the single phase flow approach. The shear stresses are shown with respect to time at the point of maximum wall shear stress at x=L2 (on the right).

Discussion
Although the drag and pressure forces are the ones that primarily influence the motion of the RBCs The pressure obtained with higher RBC was closer to the sigle-phase results. Since both CFD approaches use the same numerical algorithm for solving Navier-Stokes (presure-based, PIMPLE algorithm and same numerical schemes for the discritisation), this difference is attributed to the presence of RBC and the different forces included and the calculated shear stresses [42]. The results here are in agreement with the experiments of Freidoonimehr et al., [43] who studied lesions with much shorter stenoses (L s /D ≤ 1) than the ones considered here.
The effect of hematocrit is more evident for the lesions with a higher stenosis severity as shown in Fig.3. For all cases, the pressure drop increases for higher S D which is expected and in accordance to previous studies and the ROM [36,37]. In all simulations the flow pattern changes due to the varying cross previous studies [44].
The radial velocity profiles in Fig.4 were parabolic for all cases and the effect of the hematocrit was shown to have a lesser effect on the radial velocity compared to the stenosis severity. The magnitude of the axial velocity were higher for both lesions with S D = 70% with a peak velocity magnitude of 1.6 m/s for RBC=45% and 0.95 m/s for S D = 70% for RBC=45% during the systole cycle. Our results of the velocity profile are qualitatively in agreement with the ones in Jung et al. [18] who also observed similar flat profiles at the radial direction of the vessel.
At the L2 position the velocity profile has a parabolic-shape profile. This shape differs from the profiles reported in the work of [45] who observed a plugshaped flow for axisymmetric lesions of various stenosis severities including the ones studied here. The reason for this difference in the profile with the results here, is that in the case of the axisymmetric idealised lesions, the RBCs move towards the centreline of the vessel whereas in this study, due to eccentricity, the red blood cells exhibit a core with high RBC levels at the centreline, but also at regions near the wall at the lower part of the vessel.
The shear stresses observed with the single-phase and the two-phase approaches have the same trend with the inlet velocity. For the same stenosis, the single phase approach predicted higher values for shear stress (Fig.4-right).
The WSS when the cross-sectional areas is smaller becomes higher, having high shear rates near the vessel wall. The non-Newtonian shear thinning is responsible for this difference, and is related to a lower shear rate in the case the flow is considered to be multiphase. The approach described here treats the RBCs (which mostly contribute to the haemodynamic forces) as a shear-thinning fluid with a viscosity which is a function of the volume fraction [46].
In this study it is assumed that the red blood cells are deformable and non-spherical. Following previous studies the dispersed phase consists of rigid spherical particles without considering the effects of the normal stress effects [47]. The stress tensor for the RBCs is, where τ ij,d is the partial stress tensor for the RBC with a corresponding shear rate S ij,d . The constant k 2 accounts for the shear-thinning effect of RBCs, calculated by [48],γ is the RBC shear rate, µ ∞ is an asymptotic viscosity for the case of an infinite shear rate, µ 0 is the viscosity for the case of zero shear rate and k is a parameter that describes the shear thinning aspects of the RBCs.
The shear stresses at the L2 point in Fig.4 were also higher for higher RBC concentration as shown in Fig.6 for the same stenosis, and higher than those of the single-phase model. Depending on the region considered inside the vessel, WSS varied with the L3 giving lower WSS in the downflow region in Fig.7 compared to L2. In this downflow position, the high shear stresses regions formed in the areas of higher axial velocities and low RBC (shown in Fig.6).
The low viscosity in those areas might be also a parameter that contributes to the formation of these low-RBC regions as observed in [18]. For higher initial RBC levels from 35% to 45%, the relative RBC volume fraction generally increased in the downflow location L3 but WSS decreased near the wall, while preserving higher values in the location L2 for both cases for different stenosis severity.
The curvature-induced secondary flow profiles exhibited vortices in the flow as shown in Fig.5, with a recirculation pattern that varied during the cardiac cycle. As shown in Fig.4 high velocity gradients were observed at the near wall regions with a zero gradient closer to the centre. Large secondary recirculating flows were observed in all cases here with similar patterns for the different stenosis severity and initial RBC level (Fig.5). These recirculation zones might contribute to large pressure fluctuations locally, and lead the migration of the RBCs towards locations other than the centreline [49]. At the upstream position of the stenosis (L1) a lesser RBC migration is observed with similar relative profiles for all the vessel models studied here.
On the other hand, for all the vessels, low RBC volume fraction values were generally obtained at the downstream location (L3) compared to L2, where lower WSS led to a RBC migration.
Areas with lower WSS such as areas with flow recirculation might potentially form atherosclerotic plaque [50]. The RBC patterns at the downstream position were less uniform than the ones in the other locations, due the impact of turbulence, the impact of which, increases in more stenosed vessels, and was also reported in [51].
Local fluctuations in the red blood cell distribution at the maximum curvature and downstream locations together with the changes in viscocity might be a factor that could potentially lead to a limited intake of oxygen at the arterial endothelial surface causing hypoxia [52]. Furthermore, nitric oxide (NO) is captured by red blood cells and scavenged inside through hemoglobin. Hence, the changes in the RBC distribution due to the particulated deposition might lead to changes in the NO transport by the red blood cell hemoglobin which has been reported to cause atherogenesis and atherosclerosis [53,54].
The WSS is a significant haemodynamic parameter that is usually considered to study the platelet and vascular behaviour in blood flow through vessels. Areas with high WSS might have increased risk for diseases such as thrombosis. Areas with low or oscillating WSS which are more evident at the lesion models with stenosis severity of 70% were observed as shown in Fig.7.
The oscillatory patterns together with the secondary flow that occurs inside the stenosis and especially at regions with high curvature, have a significant contribution to the RBC deposition in the lesion. These changes downstream the maximum curvature position can potentially lead to a reduced WSS. In the areas of lower WSS, the endothelium might be more easily compromised by other particulates (monocytes and proteins) which might cause thrombosis and atherogenesis [55,56].

Conclusion
The paper here presents a novel study for the impact of hematocrit for calculating pressure drop through eccentric tubular lesions considering the most important fluid forces such as lift, drag, virtual, gravity, which has not been previously presented in the literature. Our results showed that pressure drop is not strongly affected from single phase, and neither is mean velocity thus confirming single phase modelling is adequate for estimating these variables.
This implies the validity of the majority of previous studies. On the other hand, Options for allowing nonspherical particles using the shape factor concept [18] could provide more insights for the haemodynamic parameters studied here.

Conflict of Interest
The authors declare no conflict of interest.