Influence of hematocrit-dependent blood viscosity on hemodynamics of mitral paravalvular leak

Igor Szumny; Magdalena Kopernik

doi:10.37190/abb/231653

1. Introduction

An accurate computational representation of cardiac function requires simplified yet physiologically justified assumptions regarding cardiac geometry and myocardial kinematics. Although mean anatomical parameters are commonly adopted to facilitate numerical modeling, these quantities are inherently dynamic and influenced by multiple physiological factors. Among them, hematocrit, defined as the fraction of red blood cells in blood volume, represents a particularly variable parameter and is sensitive even to mild dehydration caused by insufficient fluid intake [1]. Variations in hematocrit, while not pathological per se, may substantially affect blood rheological properties, especially in flow-restricted or high-velocity regions [2].

Paravalvular leak (PVL) provides a clinically relevant framework for investigating the interaction between blood rheology and intracardiac flow dynamics. In the setting of prosthetic mitral valves, regurgitant flow through a localized discontinuity between the valve and the surrounding tissue may reach high velocities strongly governed by leak geometry and blood viscosity [3,4]. Such jet flows can generate elevated shear stresses that may induce mechanical damage to erythrocytes and result in hemolysis [5]. During ventricular systole, blood is simultaneously ejected through the aortic valve and retrogradely through the PVL channel, producing disturbed flow patterns that deviate from physiological mitral inflow conditions. This configuration constitutes a well-defined computational framework for assessing the sensitivity of intracardiac hemodynamics to geometric and hematological factors [6].

Key parameters considered in cardiac flow simulations include velocity fields, mass flow rate, pressure gradients, and shear stress distributions. Among these, shear stress magnitude and exposure time play a dominant role in hemolysis development, whereas pressure and flow rate reflect global cardiac performance and energetic efficiency [7]. Consequently, time-resolved numerical simulations are essential for capturing transient hemodynamic conditions relevant to blood damage assessment.

Accurate representation of blood behavior requires appropriate rheological modeling. While Newtonian assumptions offer computational efficiency, they fail to capture the shear-dependent viscosity intrinsic to blood flow. For this reason, non-Newtonian formulations, such as the Carreau–Yasuda and Walburn–Schneck models, are widely employed in cardiovascular simulations [8,9]. Yield-stress-based models, including the Casson and Quemada formulations, are particularly relevant in regions characterized by low shear rates or flow recirculation [10]. Importantly, select models explicitly incorporate hematocrit as a governing parameter, rendering them well suited for studies focused on hematological variability [11].

Although blood density is commonly assumed constant in computational analyses, shear rate and hematocrit remain the primary determinants of blood rheology [12]. These factors might influence erythrocyte aggregation and deformability, suggesting a potential governance over non-Newtonian flow behavior [13]. Incorporating hematocrit-dependent viscosity enhances the realism and predictive capacity of numerical models, particularly in disease-specific and confined flow environments [14].

Previous investigations of mitral PVL have demonstrated that simplified ventricular geometries preserve hemolysis-related shear metrics while significantly reducing computational cost [15]. Detailed numerical studies further indicate that both the spatial distribution and temporal persistence of elevated shear stresses within PVL channels critically influence hemolytic potential [7,16]. However, most PVL-oriented simulations assume fixed blood rheology, whereas hematocrit-dependent viscosity models have predominantly been applied in non-valvular cardiovascular settings, including vascular stenoses and aneurysms [9,17–19].

Consequently, the systematic integration of hematocrit-dependent blood rheology into intracardiac computational fluid dynamics models of mitral PVL remains insufficiently explored. Addressing this gap may improve the mechanistic understanding of PVL-related hemolysis and support more patient-specific risk assessment strategies.

2. Materials and methods

The left ventricle generates systemic blood flow by ejecting blood through the aortic valve during systole, while diastolic filling occurs via the mitral valve. In prosthetic valve recipients, PVL may develop due to incomplete sealing between the prosthesis and native tissue, leading to high velocity regurgitant jets. The flow characteristics within PVL are primarily governed by leak geometry and fluid properties and are associated with elevated shear stresses that may induce erythrocyte damage and hemolysis [3,5]. This configuration constitutes a well-controlled framework for analyzing pathological intracardiac flow using computational fluid dynamics [6].

Hematocrit, defined as the volumetric fraction of erythrocytes in the total blood volume, is a key determinant of blood rheology. In healthy adults, it typically ranges between ∼0.3 and 0.55 (Figure 1) [20], depending on sex and physiological conditions. Even non-pathological variations, such as those induced by dehydration, can significantly alter blood viscosity [1]. Changes in hematocrit affect flow behavior particularly in confined or high shear regions, modifying velocity gradients and stress distributions [21], and directly influence erythrocyte aggregation and deformability, which govern non-Newtonian effects [13].

Distribution of hematocrit levels in the general population.

Accurate hemodynamic modeling therefore requires non-Newtonian constitutive laws. While commonly used models such as Carreau–Yasuda capture shear thinning, they do not explicitly account for hematocrit [8,10]. In contrast, models such as Walburn–Schneck and Quemada incorporate hematocrit as a governing parameter [11]. The Quemada model is particularly suitable for representing suspension behavior of blood, as it directly relates viscosity to shear rate and red blood cell concentration, making it relevant for simulations addressing hematocrit-dependent changes in PVL hemodynamics [14,21], as follows:

(1)

η = η_{F} {1 - \frac{1}{2} \frac{k_{0} + k_{\infty} \sqrt{\hat{\dot{γ}}}}{1 + \sqrt{\hat{\dot{γ}}}} H_{c})}^{- 2}

(2)

\hat{\dot{γ}} = \frac{{\dot{γ}}_{eff}}{{\dot{γ}}_{c}}

where

\dot{γ}

is the shear rate tensor;

{\dot{γ}}_{eff}

is the effective shear rate;

{\dot{γ}}_{c}

is the critical shear rate, above which elasticity can be neglected;

η_{F}

= 1.2 mPa∙s is the plasma viscosity.

The values of the parameters ${\dot{γ}}_{c}$ , $k_{0}$ , and $k_{\infty}$ are calculated as follows:

(3)

k_{0} = 25.608 {H_{c}}^{2} - 35.091 H_{c} + 14.981

(4)

k_{\infty} = 1.3125 {H_{c}}^{2} - 2.3001 H_{c} + 2.6884

(5)

{\dot{γ}}_{c} = 30.443 {H_{c}}^{2} - 13.528 H_{c} + 1.49

The hematocrit range (0.30–0.55) investigated in this study was selected to represent physiologically relevant values commonly observed in the adult population and associated with normal inter-individual variability as well as hydration status. Although severe hemolysis may lead to substantially reduced hematocrit values, reported cases indicate that hematocrit may decrease below 0.30 and, in extreme conditions, even approach 0.10–0.20 [22–25]. However, such reductions are generally considered a consequence of hemolysis rather than a primary factor governing its initiation. Therefore, the present study focused on the influence of physiologically realistic hematocrit variations on hemodynamic conditions associated with hemolysis risk. Extending the investigated range toward severe anemia would require coupling the computational fluid dynamics (CFD) model with erythrocyte destruction and hematocrit transport models, which is beyond the scope of the present work.

A three-dimensional computational model of intracardiac flow with a mitral PVL was developed using the finite volume method with time step equal to 0.00025 s, solving the unsteady Navier–Stokes equations.

(6)

\nabla \cdot u = 0

(7)

ρ (\frac{\partial u}{\partial t} + u \cdot \nabla u) = - \nabla p + \nabla \cdot (μ_{eff} \nabla u)

where u is the velocity vector,

p

is the pressure,

ρ

is the fluid density,

μ_{eff}

is the effective dynamic viscosity (molecular viscosity + turbulent viscosity), and

t

– time

.

The simulations were performed using a pressure-based transient solver implemented in ANSYS Fluent. Pressure–velocity coupling was handled with the Pressure-Implicit with Splitting of Operators algorithm. Second-order spatial discretization schemes were applied to all transport equations, while first-order temporal discretization was used to maintain numerical stability during dynamic mesh deformation. The simulations were conducted with a time step of 2.5 × 10⁻⁴ s and a maximum of 100 iterations per time step. The selected time step ensured that the Courant number remained below 40 throughout the simulation (based on [26]). Dynamic mesh calculations were performed using smoothing, layering, and remeshing techniques to preserve mesh quality during ventricular motion.

The computational domain was constructed as a simplified three-dimensional computer-aided design (CAD) model of the left heart, including the left atrium, left ventricle, and a PVL channel (Figure 2a). The ventricular and atrial volumes were approximately 114 and 41 ml, respectively, with a stroke volume of 62 ml. The diameter of the aorta was set to 25 mm [15] and diameters of pulmonary veins were set as: 17.4 mm, 16.5 mm, 16.95 mm, and 14.45 mm [15]. The PVL orifice (7.9 mm²) was modeled as a localized gap along the mitral annulus, enabling formation of a high velocity regurgitant jet.

Computational model: (a) CAD geometry of the left heart with PVL; (b) computational mesh with local refinement in the PVL region and boundary layers; (c) ventricular motion model based on the centroid method; and (d) mesh convergence study.

The mesh (Figure 2b) consisted of an unstructured grid with approximately 687,537 nodes and 164,845 control volumes. Local refinement was applied in the PVL region and near walls, and three boundary layers were introduced to capture near-wall gradients. Mesh quality metrics ensured stability, with skewness below 0.79 and orthogonal quality above 0.21.

A mesh convergence study (Figure 2d) was performed using 110 successively refined meshes on Helios HPC supercomputer with 75264 AMD Zen4 computing cores and 200 TB of DDR5 RAM. Differences in key parameters, including maximum velocity and shear stress, were below 1%, confirming mesh convergence. The mesh was selected as a trade-off between accuracy and computational cost. To find the most stable mesh configuration, the standard deviation was used as the convergence criterion for each mesh size and number of layers. The mesh features an average global cell volume of 19 mm³, which is refined to approximately 0.03 mm³ inside the leak. The leak cross-section contains roughly 290 cells (4,640 cells in whole PVL), fully satisfying the requirement for accurate velocity profile representation.

Blood density was modeled using the Tait equation [26,27].

(8)

{\frac{ρ}{ρ_{0}})}^{n} = \frac{K}{K_{0}}, K = K_{0} + n Δ p

where

ρ

is the fluid density,

ρ_{0} = 1,060 {kg ∙ m}^{- 3}

is the reference density,

K

is the bulk modulus at pressure

p

K_{0} = 2.2 \cdot 10^{9} Pa

is the reference bulk modulus,

n

is the exponent describing compressibility,

Δ p = p - p_{0}

is the pressure difference relative to reference pressure, and

p_{0}

is the reference pressure.

Turbulence was modeled using the k–ω shear stress transport model with parameter based on the study by Varghese and Frankel [28]. A constant pressure of 30 mmHg was imposed at the inlet, while a time-dependent pressure was applied at the outlet to reproduce physiological pulsatile flow. The cardiac cycle was assumed to last 1 s (60 bpm), consisting of a systolic phase (∼0.3 s) followed by diastolic filling with E-wave and A-wave components. The prescribed flow ensured consistency between pressure, flow rate, and ventricular volume.

Ventricular contraction was modeled using a centroid-based kinematic approach (Figure 2c), in which node displacement is defined as follows:

(9)

P_{i} = P_{start} + f \cdot (S - P_{start}) \cdot \frac{d - d_{\min}}{d_{\max} - d_{\min}}

where P _i is the current position of mesh node,

P_{start}

is the initial position of mesh node,

S

is the centroid position of the ventricle,

f

is the contraction scaling factor,

d

is the distance between node and characteristic damping point,

d_{\min}, d_{\max}

is the minimum and maximum distance from the characteristic damping point

.

Characteristic damping point was defined for controlling the unchangeability of aorta and PVL.

Based on the theoretical dependence between $f$ and $V$ , the model was created as

(10)

V = a \cdot f^{β} + V_{start}

where

V

is the ventricular volume,

V_{start} = 185.4378629

is the starting ventricular volume,

f

is the contraction parameter, fitted empirically parameter

a = - 64.1844

, and fitted empirically parameter

β = 0.7092

Additionally, the expected volume change resulting from flow was calculated analytically from the study by Seo et al. [29] as

(11)

\int_{t_{\min}}^{t} Q (t) d t = \begin{matrix} \frac{Q_{s}}{2} t - t_{peak} - \frac{t_{peak}}{π} Sin \frac{π t}{t_{peak}})], t < t_{peak} \\ \frac{Q_{s}}{2} t - t_{peak} + \frac{t_{\max} - t_{peak}}{π} Sin \frac{π (t - t_{peak})}{t_{\max} - t_{peak}})], t \geq t_{peak} \end{matrix}

where

Q (t)

is the instantaneous mass flow rate,

Q_{s} is the

characteristic systolic mass flow amplitude,

t

is the time

,

t_{\min}

is the initial time of the cardiac cycle,

t_{peak}

is the time of pressure peak in aorta,

t_{\max}

is the total duration of the cardiac cycle, and

\sin (\cdot)

is the sinusoidal function describing pulsatile flow variation.

Equation (11) describes a physiologically motivated pulsatile flow waveform, with smooth transitions between cardiac phases ensured by sinusoidal functions.

Equation (12) is related to mass conservation

(12)

\int_{t_{\min}}^{t} Q (t) dt = Δ V (t) \cdot ρ

where

Q (t)

is the mass flow rate,

t

is the time,

Δ V (t)

is the change in ventricular volume, and

ρ

is the fluid density.

This formulation ensures consistency between kinematic wall motion and mass flow variation during the cardiac cycle.

Shear stress was calculated as

(13)

τ = μ_{eff} \dot{γ}

where

τ

is the shear stress,

μ_{eff}

is the effective viscosity (molecular viscosity + turbulent viscosity), and

\dot{γ}

is the shear rate and compared to a critical threshold of 300 Pa to assess hemolysis risk within the PVL region.

The model assumes a simplified geometry and prescribed ventricular motion without fluid–structure interaction; therefore, the applied boundary conditions and material models represent approximations of physiological conditions. Despite these limitations, the model captures the key mechanisms governing PVL flow and shear stress generation.

3. Results

The developed computational model enabled detailed analysis of flow dynamics within the left heart and the PVL region under pulsatile conditions. All reported shear stress values correspond to the maximum shear stress within the PVL region, where the highest velocity gradients occur.

The time-dependent pressure distribution within the left heart is shown in Figure 3. The results demonstrate a clear evolution of pressure during the cardiac cycle. At the initial phase (t = 0 s), relatively uniform low-pressure regions are observed. During systole, a significant increase in pressure occurs within the left ventricle, creating a strong pressure gradient between the ventricle and the atrium. This pressure difference drives blood flow both through the aortic valve and through the PVL channel. The highest values of pressure are localized within the ventricle, while the atrium remains at lower pressure levels throughout the cycle.

Time-dependent pressure distribution in the left ventricle during the cardiac cycle at selected time instances: t = 0 s, t = 0.2 s, and t = 0.25 s.

Figure 4a illustrates the relationship between hematocrit and maximum shear stress within the PVL region. A strong correlation (R ² ≈ 0.99) is observed, with increasing hematocrit leading to higher shear stress values. This behavior results from increased effective blood viscosity, which enhances velocity gradients within the regurgitant jet. At higher hematocrit values, shear stress exceeds the critical threshold of 300 Pa, indicating an elevated risk of hemolysis. This behavior results from increased effective blood viscosity, which elevates shear stress levels within the regurgitant jet even under similar velocity gradients. While the critical threshold of 300 Pa is exceeded across all analyzed hematocrit values, higher hematocrit levels further amplify these stresses, potentially exacerbating the risk of hemolysis.

(a) Relationship between hematocrit and maximum shear stress in the PVL region. (b) Time required to reach the critical shear stress threshold (300 Pa) in the PVL region as a function of hematocrit.

The influence of hematocrit on the temporal behavior of shear stress is presented in Figure 4b. The results show that increasing hematocrit leads to an earlier onset of critical shear stress conditions. Specifically, the time required to reach the hemolysis threshold (300 Pa) decreases as hematocrit increases. This indicates that higher hematocrit not only increases the magnitude of shear stress but also accelerates the occurrence of potentially hemolytic flow conditions.

The spatial distribution of shear stress for different hematocrit values is shown in Figure 5. In all cases, the highest shear stress values are localized at the inlet of the PVL channel and along the channel walls, where steep velocity gradients are present. A concentrated high-velocity jet is observed, with clearly defined shear layers separating the jet from the surrounding fluid.

Distribution of shear stress in the PVL region for different hematocrit values (φ = 0.3–0.55).

As hematocrit increases, the shear stress distribution becomes more localized, with higher peak values and sharper gradients. Lower hematocrit values result in more diffuse stress fields and reduce peak stresses. These spatial patterns are consistent with the quantitative trends observed in Figures 4a and b.

4. Discussion

The results of the present study demonstrate that flow through a mitral PVL is dominated by the formation of a high velocity regurgitant jet, which governs the spatial distribution of shear stress. This observation is consistent with previous CFD studies, where PVL-induced jets were identified as the primary source of localized high shear stress and hemolysis risk [6,30]. The confined geometry of the leak promotes steep velocity gradients at the jet boundary, leading to localized regions of elevated shear, in agreement with earlier findings reported for patient-specific PVL simulations [5].

A key outcome of this study is the quantitative relationship between hematocrit and shear stress within the PVL region. The results show a strong linear increase in maximum shear stress with increasing hematocrit, which is consistent with the role of viscosity as a dominant factor in blood flow resistance and shear generation [11,12]. Similar trends have been reported in vascular flow studies, where higher hematocrit values led to increased shear stress in confined geometries such as stenoses and aneurysms [2,18]. The present findings extend these observations to intracardiac flow, confirming the importance of hematocrit-dependent rheology in PVL simulations.

The temporal analysis revealed that increasing hematocrit leads to an earlier onset of critical shear stress conditions. In other words, the time required to reach the hemolysis threshold (300 Pa) decreases with increasing hematocrit. This indicates that hematocrit affects not only the magnitude of shear stress but also its temporal development, which is particularly relevant for hemolysis risk assessment. This finding complements previous studies emphasizing the role of both shear stress magnitude and exposure time in blood damage mechanisms [7].

Despite similar global flow patterns across all cases, significant differences were observed in local flow structures. Higher hematocrit values produced more concentrated shear layers and stronger localization of stress within the PVL channel, as seen in the spatial distributions. In contrast, lower hematocrit values resulted in more diffused flow fields with reduced peak stresses. This response is consistent with the non-Newtonian nature of blood and the influence of erythrocyte interactions on flow structure [13]. Together, the quantitative trends and spatial distributions confirm that even moderate physiological variations in hematocrit can significantly affect local hemodynamics.

The applied ventricular motion model enabled reproduction of physiologically realistic pulsatile flow conditions while maintaining computational efficiency. This approach is consistent with simplified kinematic models used in previous intracardiac simulations [15]. Although the model does not include fully coupled fluid–structure interaction, it captures the dominant flow features relevant to PVL response, particularly those governing shear stress generation.

Hematocrit significantly affects both the magnitude and onset of elevated shear stress. Future studies should incorporate patient-specific geometries, exposure-time-based hemolysis models, and fully coupled fluid–structure interaction.

5. Conclusion

The developed computational model demonstrated that flow in a mitral PVL is dominated by a high‑velocity jet generating localized regions of elevated shear stress, particularly during systole.
Hematocrit was identified as a key parameter influencing both the magnitude and onset of shear stress, with higher values leading to increased peak stresses and earlier attainment of the hemolysis threshold.
These findings highlight the importance of incorporating hematocrit‑dependent blood rheology in CFD analyses and suggest that patient‑specific hematological conditions may significantly affect the assessment of hemolysis risk in PVL.

Acknowledgements

We gratefully acknowledge the Polish high-performance computing infrastructure PLGrid (HPC Center: ACK Cyfronet AGH) for providing computational resources and support within computational grant no. PLG/2026/019728.

Funding information

This research was conducted within the framework of statutory activities no. 16.16.110.663 at AGH University of Krakow.

Author contributions

Igor Szumny: investigation, numerical simulations, literature review, and responses to reviewers. Magdalena Kopernik: conceptualization, methodology, writing – original draft, and writing – review & editing. Both authors contributed equally to this work. Both authors reviewed and approved the final version of the manuscript.

Conflict of interest statement

Authors state no conflict of interest.

Influence of hematocrit-dependent blood viscosity on hemodynamics of mitral paravalvular leak

Full Article

1. Introduction

2. Materials and methods

Figure 1

(1)

(2)

(3)

(4)

(5)

(6)

(7)

Figure 2

(8)

(9)

(10)

(11)

(12)

(13)

3. Results

Figure 3

Figure 4

Figure 5

4. Discussion

5. Conclusion

Acknowledgements

Funding information

Author contributions

Conflict of interest statement

Paradigm

My account