Blood rheology modeling effects in aortic flow simulations
Alexander Fuchs1,2,*, Niclas Berg2, Lisa Prahl Wittberg2
1. Department of Radiology, Karolinska University Hospital, Stockholm and Department of Medical and Health Sciences, Linköping University, Linköping, Sweden.
2. FLOW & BioMEx, Department of Engineering Mechanics, KTH, 100 44 Stockholm, Sweden.
*Corresponding Author
Email for correspondence: alex@mech.kth.se
Abstract
Purpose: The purpose of the study is to assess the importance of non-Newtonian rheological models on blood flow in the human thoracic aorta.
Methods: The pulsatile flow in the aorta is simulated using the models of Casson, Quemada and Walburn- Schneck in addition to a case of fixed (Newtonian) viscosity. The impact of the four rheological models was assessed with respect to the following quantities: i. Magnitude of the viscosity relative to a reference value (the Newtonian case) and the relative mean deviation from that value. ii. Mechanical kinetic energy,
vorticity, viscous dissipation rate. iii. WSS and its time derivative. iv. WSS-related indicators; OSI, TAWSS and RRT.
Results: The flow in the thoracic aorta is characterized by shear-rates leading to an increase in viscosity by a factor of up to six. The different models had negligible impact on the kinetic energy and viscous dissipation rate. The effect on WSS related parameters was quantified and was found to be modest. Largest effect was observed for low shear-rates (below 100 s-2).
Conclusions: The choice of a non-Newtonian model is important whenever the flow is viscosity dominated.
Blood flow in larger arteries is weakly dependent on viscosity and can be handled by a model with weak dependence on shear-rate (e.g. Quemada or Newtonian). Blood flows with regions with low shear-rate and strong temporal variation requires rheological models that better account for low shear and explicitly includes temporal variation effects.
Keywords: Whole blood viscosity, non-Newtonian fluid, Thoracic aorta, Wall shear-stress, viscous losses.
Introduction
Blood is a fluid that is mainly composed of water but also a wide range of cells, micelles and molecules of widely different sizes. The multiple functionality of blood and its components includes transporting
substances needed by the different organs, playing a major role in maintaining optimal environment (e.g. pH, temperature) in the body as well as defending the body against microorganisms and stopping bleeding.
Hence, blood composition may change in terms of numbers and types of cells, micelles and molecules as response to needs. This adaptive behavior is rather quick, implying that for each individual, the blood rheological properties may vary during the day. Also, pathological variations of the blood constituents affect blood viscosity. For example, RBC size distribution (or red blood cell distribution width (RDW)) has been found to be a predictor of morbidity and mortality (cf Lippi et al (2017) and Danese et al (2015)). Non- uniform distribution of RBC within the circulatory system contributes to the difficulty in determining whole blood viscosity (WBV). Presence of macromolecules in a Newtonian fluid implies that the mixture
characteristics becomes non-Newtonian due to the macromolecules’ ability to energy take-up and release during a finite time (so called relaxation time). In contrast, Newtonian fluids respond instantly to the forces acting on it. Hence, the constitutive relation between stress and strain is linear and contains no response time.
Blood rheology has been the subject of a large number of studies found in literature. Often, WBV measurements have been carried out using standard (shear) rheometers (cf Cowan et al (2012)), implying that the flow is laminar and that the shear-rate is well-defined and depending solely on the rotation rate (cf Agarwal et al (2019)). Yamamoto et al (2010) used a compact-sized falling needle rheometer on fresh blood samples measuring the relationship between the shear stress (τ) and shear-rate (γ). The study identified three
typical “regions”: The “Casson” region for a low shear-rate range (below 140 s−1), the transition region (up to160 s−1), and the Newtonian fluid region for higher shear-rates (above 160 s−1). The range of human blood viscosity was found to be in the range of 5.5 to 6.4 mPa s, and 4.5 to 5.3 mPa s for males and females, respectively. Blood is a non-Newtonian fluid and exhibits a shear-thinning behavior, i.e. reduction in viscosity with increasing shear-rate. As indicated earlier, WBV differs among individuals and in the same individual at different points in time. Wang et al (2003) measured WBV for a group of healthy individuals showing that both inter- and intra-individual variations were higher in the morning than later in the day for all shear-rates. One factor influencing WBV is the types and number of lipoproteins (Lowe (1982)).
The flow of blood in large arteries differs significantly from the corresponding situation in small arteriole. In large arteries, the blood is often considered as a homogenous mixture of RBC and solvents in the plasma.
The blood can be assumed to behave as a continuum although the blood components are not homogenously distributed, given that the length scales of blood cells and mean distance among the cells are much smaller as compared to the flow scales under consideration. Obviously, this approximation is not appropriate for vessels of the size of the RBCs themselves (or slightly larger). In smaller arteries, a clear RBC free layer forms that may lead to separation of RBC:s from plasma at vessel bifurcations. In even smaller vessels, RBC move individually and the continuity assumptions, used in deriving the common Navier-Stokes equations, are not valid for the RBC “phase”. As consequence, several modeling frameworks for blood flow have been proposed. The most common formulation uses a mixture model without or with a transport equation for the RBC concentration (hematocrit), where the latter requires a constitutive model for the mixture and a diffusion model for the RBC transport. This approach is a simplification of a two-fluid model in which the plasma (which is a Newtonian fluid in shear flow but viscoelastic in pure extensional flow (Brust et al (2013)) and the RBC phase are assumed to be handled as two separate continuums interacting with each other. This approach requires a constitutive model for the viscosity properties of each of the phases accounting for diffusion within each phase in addition to the effects of interphase forces.
Over past decades, different types of whole blood viscosity models based on rheological data have been proposed where power-law based models were suggested due to the resemblance between blood and other shear-thinning fluids. A recent review of most of existing models of whole blood viscosity was presented by Hund et al (2017). All models include several parameters determined by fitting the models to measurements.
For example, Marcinkowska-Gapinska et al (2007) measured the viscosity of 100 whole blood and plasma samples over a range of shear-rates (between 0.01 to 100 s-1). Applying the measured data to the models of Casson, Ree-Eyring and Quemada, the Quemada model was found to provide the best fit. Several studies have assessed and compared different rheological models, resulting in discrepancy among the conclusions.
Gallagher et al (2019) considered models of Bird, Carreau, Cross and Yasuda. The aim of the study was to address the problem of inferring model parameters by fitting to experiments. By refitting published data, families of parameter sets capturing the data equally well were identified, implying that such parameters cannot be used to draw conclusions about fluid physical properties. In order to assess the sensitivity of the models, random perturbations were added to the measured data from which new model parameters were derived, showing that the problem was inherent to the models considered. The effects of the Casson and Carreau-Yasuda models on a steady and oscillatory 2D-flow in a straight and curved pipe geometry was studied by Boyd et al (2007) using the Lattice Boltzmann method. Differences in velocity and shear was found at low Reynolds and Womersley numbers, although rather small in terms of velocity profiles. Akhera et al (2017) simulated blood flow in arteriovenous (AV) fistula using reconstructed 3D geometries and the models of Quemada and Casson. The results displayed no major differences in the flow field and the flow characteristics. Instead, the shape of the geometry was found to be far more important for the WSS
distribution as compared to the effect due to the rheological model. Investigating the flow in a patient-specific aorta model, Karimi et al (2014) applied nine rheological models (three Casson model variants, Carreau, Carreau-Yasuda, Cross, Power-law, Modified Power-law, and Generalized Power-law Parameters), focusing on WSS and the deviation of the computed viscosity as compared to a reference viscosity (3.45 10-3 Pa s).
The largest differences in WSS were located near the branches and found more pronounced at low flow rate (diastole) where the viscosity was considerably larger than the value used in the Newtonian computations.
Karimi et al (2014) concluded that various rheological models may yield equally good results with the exception of the Cross model. A similar conclusion was presented by Johnston et al (2004), numerically investigating the effects of six rheological models (including a Newtonian) on WSS during the cardiac cycle.
The results were comparable for all models under steady flow conditions and at unsteady mid-range flow
velocities (around 0.2 m/s). Moreover, Newtonian model was reported as a good approximation in regions of mid-range to high shear, whereas for low shear-rates, the power-law model was found more appropriate.
However, in the numerical study of a pipe flow by Jahangiri et al (2018), all rheological models applied (Carreau, Carreau-Yasuda, modified Casson, Power law, generalized power law, and Walburn-Schneck) except the generalized power law model resulted in graphically the same axial velocity profile and WSS Returning to patient specific related geometries, Liepsch et al (2018) studied the impact of non-Newtonian viscosity models on the hemodynamics in a cerebral aneurysm. The flow in the aneurysm was computed using Newtonian, Power Law, Bird-Carreau, Casson and Local viscosity model. Although similar flow patterns were observed both for Newtonian and non-Newtonian models, a quantitative comparison done on a group of monitoring points, revealed an average difference between the Newtonian and non-Newtonian models of about 12%. Furthermore, in low speed regions the differences were even larger (20% to 63%).
Skiadopoulos et al (2017) considered Newtonian, Quemada and Casson blood viscosity models for simulating pulsatile flow in patient specific geometry of the iliac bifurcation. The effect of the rheological models was monitored through the WSS distribution, magnitude and oscillations, and viscosity behavior as function of the shear-rate. In addition to the commonly used WSS-related indicators (OSI and TAWSS), the studied
considered the (wall) area averaged WSS and the corresponding area averaged shear-rate. The magnitude of WSS and its oscillations were found to depend on the shear-rate and the rheological model. The main conclusion of Skiadopoulos et al (2017) was that the Newtonian approximation is mostly applicable for high shear and flow rates. The Newtonian model was found to overestimate the possibility of formation of
atherosclerotic lesions in regions with oscillatory WSS. In a recent paper, Mendieta et al (2020) compared the effect of rheological models on the WSS-related parameters in stenotic carotid flows. Four rheological models were considered (Carreau, Cross, Quemada and Power-law). Largest differences between the Newtonian and non-Newtonian models were noted for OSI of about 12% (in terms of maximum and mean values). Regarding TAWSS, the difference was less than 6%, except for the Quemada model for which the difference was as much as 26%. The authors conclude that the assumption of a Newtonian model can be reasonable, however, non-Newtonian models were found necessary in low TAWSS regions.
Most of the above cases considered laminar and/or transitional flow regimes. In contrast, Molla and Paul (2012) studied a turbulent flow in a channel with a constriction, computed by Large Eddy Simulations. Five rheological models (Power-law, Carreau, Quemada, Cross and a modified-Casson) were compared in terms of peak shear-rate, mean shear-stress and pressure, re-circulation zones and turbulent kinetic energy. The main finding was that the non-Newtonian viscosity models extended the post-stenotic re-circulation region and reduced the turbulent kinetic energy downstream of the stenosis. In the simulations, the shear-rate was limited to lower than 100 s-1, i.e. in the range where the viscosity differs strongly from the high shear-rate range that in turn implies higher levels of viscosity and lower level of turbulent kinetic energy (TKE). In terms of TKE, the differences between the rheological models were rather modest.
Rheological blood models are calibrated to in-vitro measured data. However, when applied to patient- specific simulations, the main goal with the models is to be able to account for important fluid physics and capture relevant clinical observations. An equally important question is the sensitivity of the models and their impact on parameters of clinical significance. The aim of this study is to address the issues related to the impact of rheological models on hemodynamically related parameters (i.e. mechanical kinetic energy, vorticity, viscous dissipation rate and WSS with time derivatives) when applied to blood flow in the thoracic aorta.
Numerical method and Case-setup
The geometry of the thoracic aorta was derived from a computed tomography angiography (CTA) study of a healthy patient. The computational domain consisted of the ascending aorta, the aortic arch, the descending aorta and the tree main branching arteries; the Brachiocephalic Artery (BCA), left Common Carotid Artery (LCCA) and left Subclavian Artery (LSCA)). The blood was assumed to be incompressible with non- constant bulk density depending on the concentration (hematocrit) of Red Blood Cells (RBC). The effects of other blood components were neglected. Thus, the blood was modelled as a non-homogenous mixture satisfying conservation of mass (eq. 1a) and momentum (eq. 1b).
( i j)
i i
j i j j
u u u p u
t x x x x
(1a)
i 0
i
u
t x
(1b)
Where are the density and viscosity of the mixture respectively, p is the pressure and ui is the Cartesian velocity component in the i:th direction. The bulk viscosity depends on the local blood composition and the local flow conditions, accounted for by applying several different mixture viscosity models; Newtonian fluid, and the models of Walburn and Schneck (1976), Casson (1959) and Quemada (1977 and 1978).
Mixture viscosity models
The simplest non-Newtonian model, used in multiple publications and available in several common software, is the power law model of Ostwald de Waele. In its simplest form, the model directly relates effective
viscosity to the shear-rate:
eff k
n1, with is the magnitude of the shear-rate tensor 𝛾 = + . The power of the shear-rate is negative for shear-thinning fluids (n < 1).The Walburn-Schneck model for the effective viscosity (eff) of the blood (mixture), is an extension of the Power law model. In the current work with =0.45, effective has the following form:
) ,
0034 . 0
max( 0.0225
eff (2)Casson (1959) suggested a model based on a calibrated power law concept. The following form and parameters were used herein.
( ) ( ) 2
min , C y
eff max
k
(3)
The mode parameters were derived for a hematocrit of =0.45 (Cokelet et al. 1963); Perktold et al.
1999).
The Quemada model (Quemada (1977, 1978)), may include two variables, namely the hematocrit () and the shear-rate (). Here we used a constant =0.45 leading to a simpler formulation:
10.225 ( )
2
eff p k(4)
with plasma viscosity p =1.32 ∙ 10-3 Pa s. k is an expression of (details may be found in Fuchs et al (2019).
Boundary conditions
The aortic domain was discretized by about 5 million computational cells, found to yield adequately accurate results. Regarding boundary conditions, no-slip conditions were set on the (rigid) walls of the thoracic aorta.
At the inlet, a time-dependent flow-rate profile derived from measured human cardiac profile was imposed (Fig 1a), representing 90 heart beats per minute (BPM) and 9 liters per minute (LPM). The flow rate in the main branches, BCA, LCCA and LSCA was set to 15%, 7.5% and 7.5%, respectively (Benim et al (2011)).
The flow out from the descending aorta was set to 70%. Four viscosity modeling approaches were employed:
a Newtonian fluid with kinematic viscosity of 3.310-6 m2/s, Casson-, Walburn-Schneck- and Quemada- models, using the parameter values given above. The hematocrit was kept at a constant value of 0.45. The variation of viscosity as function of shear-rate for the models given by equations (2-4) for hematocrit of 0.45 is depicted in Fig 1b.
The governing equation were discretized using a formally second order finite-volume scheme. The discrete equations were advanced in time using an implicit solver (with the OpenFoam 5.0). The results were processed using MATLAB, Paraview and VTK.
Viscosity model and WSS sensitivity indicators Monitoring effects of Non-Newtonian viscosity
The different non-Newtonian models lead to blood mixture viscosity coefficients varying in space and time.
To quantify the variation of the viscosity, the following indices were used (Johnston et al (2004)); the relative viscosity, IL, and the non-Newtonian importance factors, Ig-space and Ig-time capturing the space and time effects, respectively.
ref
I
L
;
2 / 1
1
)2
) , 1 (
) (
ref N
i
ref i
space g
t x t N
I
;
2 / 1
1
)2
) , 1 (
) (
ref M
j j ref
time g
t x x M
I
(5)
The impact of the local viscosity on the flow within the thoracic aorta was quantified through the spatial and temporal distribution of the mechanical kinetic energy (MKE) per unit mass and the specific viscous
dissipation rate (Eps), defined by;
1 2 i i
MKE u u and Eps ij ij (6)
where is the local kinematic viscosity.
WSS related indictors
The WSS plays a major role for the processes in the arterial wall and depends directly on the near wall viscosity. To assess the impact of the rheological models, different WSS indicators were used.
a. WSS(x,t) and TD_WSS(x,t)
ij ij
ij j
i
x t n
WSS ( , ) | |
_ ( , ) | |
t t WSS x WSS
TD
(7)
b. Time average based expression of WSS: TAWSS, OSI and RRT.
Time Averaged Wall Shear Stress (TAWSS) cf. Suo et al (2008) and Chen et al (2016), a local time averaged WSS:
0
1T
TAWSS WSS dti
T
(8)nj is the wall j-th component of the wall normal vector. Spatial variation of WSS provides an indication of the WSS level and its spatial non-uniformity (i.e. WSS-gradient), but do not contain any information about WSS temporal variation.
Oscillatory Shear Index (OSI) (cf He et al (1996), Chen et al (2016)), on the other hand, is a measure for temporal sign-change of WSS vector. OSI is defined as
0
0
| |
1 1 2
T i T
i
WSS dt OSI
WSS dt
(9)
Hence, OSI varies between 0 and 0.5. When WSSi has a constant sign (possibly oscillatory but maintaining the same sign of WSS vector), OSI = 0. When WSSi changes signs such that the integral of the positive and negative sequences are equal, OSI gets the value 0.5. Values of 0 < OSI < 0.5 indicate a sign oscillatory WSS field.
Relative Residence Time (RRT) has been formulated (cf Rikhtegar et al (2012), Gallo et al (2016) as follows):
0
1 (1 2 ) 1
T i
RRT WSS
OSI dt
T
(10)
Results
To elucidate the differences and similarities of the flow and shear stress characteristics due to the different rheological models used, the results are presented in terms of different parameters within the lumen or near the aortic wall:
i. The relative size of the viscosity coefficient, IL, and the so called non-Newtonian importance factor, Ig.
ii. The impact of the viscosity model on the mechanical kinetic energy (MKE) and the vorticity, related to the viscous dissipation rate (Eps).
iii. The impact of blood rheological model on WSS and its time derivative (TD_WSS).
iv. The impact of the rheological model on the WSS-related indicators, namely; OSI, TAWSS and RRT.
The results are presented in form of spatial and/or temporal distribution of the different parameters.
a b
c d
Figure 1: a) Inlet flow-rate profile vs time. b) Viscosity vs shear-rate for the three non-Newtonian models (equations 2-4). c) The space averaged relative viscosity IL and d) the importance factor (Ig-space(t)), plotted over the cardiac cycle.
a. Viscosity coefficient
First, the variations in the local viscosity due to the different models are considered. Space- averaging the relative viscosity (IL) and the non-Newtonian importance factor (Ig-space), enables assessment of the temporal variation. Figure 1c and d depicts these variables for the three non-Newtonian models. Obviously, for the Newtonian model IL=1 and Ig=0 and not shown in Fig 1. As noted, the Walburn-Schneck model yields the
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Time [s]
0 1 2 3 4 5 6 7
8 10-4 Flow-rate vs time: 90 BPM & 9 LPM
100 101 102 103 104 105
Shear-rate [1/s]) 0
10 20 30 40 50 60 70 80 90
100 Viscosity vs shear-rate
Quemada Casson Walburn-Schenk Newtonian
largest values as compared to the other two models. The Casson model presents a mean increase in viscosity only by a factor less than 1.4 whereas the corresponding value for the Quemada model is about 1.8. All three models have largest viscosity in early systole, decreasing during systole and increasing at late diastole. This behavior reflects the shear-thinning property. In systole and early-diastole, the shear-rates are largest, leading to lower viscosity whereas with lower shear-rate, the viscosity increases. The importance factor (Ig) also shows that the deviation from the Newtonian reference value is smallest for the Casson model and largest for Walburn-Schneck. The local minimum is about 0.1 s and about 0.45 s correspond to the instantaneous peaks of flow rate and largest shear. The local peaks at about 0.2 s is linked to the strongest of flow-rate
deceleration inducing retrograde flow and temporarily, decreasing shear-rates in parts of the domain. A similar effect is seen in the late diastolic phase.
The probability distribution of IL and Ig-space for the rheological models of Quemada, Casson and Walburn- Schneck is depicted in Fig 2. It is clearly shown that the values of the viscosities as computed by the
different models differ significantly. The largest number of the relative time-averaged viscosity is about 1.6, 1.25 and 1.7 for the Quemada, Casson and Walburn-Schneck models, respectively. The corresponding highest probability for Ig-space is at about 0.5, 0.2 and 0.5, respectively. IL values close to unity and correspondingly Ig-space close to 0 imply values close to the reference (Newtonian) viscosity. Peaks of the Quemada model indicate that during the cardiac cycle, there are regions where the viscosity due to the Quemada model is lower than that of the reference value. At low shear-rates the models yield considerably larger viscosity values as compared to the reference viscosity. This leads to IL > 1 and larger Ig-space values.
The largest viscosity values, although associated with low probability, are noted for the WS and Quemada models.
a b
c d
e f
Figure 2: The probability distribution of the relative viscosity (IL) and importance factor (Ig-time) computed by the Quemada model (a and b); Casson (c and d) and Walburn-Schneck (e and f).
A quantitative but less detailed comparison of the two viscosity parameters IL and Ig, is given in Table 1. The table shows the mean, standard deviation (STD), peak and minimal values of these parameters. The largest peak viscosity value was found with the Walburn-Schneck model. Yet, this model has smallest root mean square (RMS) deviation from the reference viscosity value. The Casson model shows the smallest mean values but relatively large peak value.
Model Quemada IL Ig
Casson
IL Ig
Walburn-Schneck IL Ig
Mean 1.536 0.538 1.215 0.215 1.690 0.631 std 0.288 0.266 0.166 0.158 0.290 0.535 max 5.755 4.745 6.002 4.991 4.480 14.352 min 0.817 0.011 0.906 0.017 0.876 0.009
Table 1: Time and space statistical values for IL and Ig, with the three non-Newtonian models.
b. Non-wall impact
Next, the impact of the models on the kinetic energy and vorticity or dissipation rate of kinetic energy are addressed. Fig 3 depict the temporal behavior of the space averaged MKE and Eps during the cardiac cycle.
Both MKE and Eps are similar for the four cases, although the Walburn-Schneck model has a lower peak dissipation rate value.
a b
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
cardiac cycle [s]
0 0.1 0.2 0.3 0.4 0.5
0.6 MKE - Rheological models: 90BPM 9LPM Casson Quemada Walburn-Scneck Newtonian
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
cardiac cycle [s]
0 0.5 1 1.5 2 2.5 3 3.5 4
Eps [m2/s3]
10-4Disipation-rate: Rheological models: 90BPM 9LPM Casson Quemada Walburn-Schneck Newtonian
Figure 3: a) Space averaged MKE and b) Eps for the four viscosity models investigated over a cardiac cycle.
The behavior of vorticity exposes the dynamics of the flow. Therefore, it is instructive to evaluate the probability distribution of this variable along with the MKE, Fig 4. Low values of MKE are found during diastole and corresponds therefore to higher probability values. The peak in MKE (at about 0.15 m2/s2) and larger MKE values are associated to systole. The flow in the aorta reaches its largest instantaneous value (of about 0.7 ∙10-3 m3/s) at peak systole. The peak MKE values for the Newtonian and the Quemada models are slightly larger than those of the Casson and Walburn-Schneck models. In terms of vorticity, the Walburn- Schneck model is characterized by a larger peak at lower vorticity levels, resulting in a lower probability for large vorticity values. As the vorticity and viscous dissipation are directly related (cf. Koh (1994)), the lower vorticity values may explain the finding of lower viscous dissipation during systole. For the same reason, similarity between the vorticity probabilities of the three other viscosity models is expected. On the other hand, the low probability at higher levels of vorticity suggests that the effect of viscosity is not prominent in the computed flow rate of 9 LPM.
a b
c d
e f
g h
Figure 4: Time-averaged probability distribution of MKE and vorticity for the four models: Newtonian (a and b); Casson (c and d); Quemada (e and f) and Walburn-Schneck (g and h). Note the different probability scales in the different frames as the integral of the distribution equals always unity.
Model WSS [Pa] WSS
t
[Pa/s]
Newtonian 2.122 ±3.294 36.349 ± 115.49 Quemada 2.044 ± 2.678 32.251 ± 90.47
Casson 2.045 ± 2.994 35.058 ± 112.17
Walburn-Schneck 2.049 ± 3.145 46.360 ± 171.82 Quemada 10.4LPM 2.145 ± 2.532 50.942 ± 99.61
Table 2: The mean and RMS of spatial and temporal statistics of WSS and its time-derivative. For
comparison, the corresponding values computed for a higher flow rate (10.4 LPM), with the Quemada model is given in the bottom row.
c. Model impact on WSS-related quantities
Viscous effects are most important near the walls and the stresses at the wall are known to be significant for the formation of wall pathologies. The time-averaged probability distributions of WSS and its temporal derivative (equation 7) are depicted in Fig 5, where the close shape similarities of the distributions of the different models can be observed. However, differences in terms of the peaks are noted in the WSS and its
time-derivative. Such peaks found at almost 2 Pa and about 20 Pa/s, respectively. Largest WSS peaks are found for the Walburn-Schneck model followed by the Newtonian model. Lowest WSS peak is noted for the Quemada model, implying the presence of larger WSS values, although with low level of probability. A common observation is that most of the aortic wall is subjected to low WSS and only a smaller portion of the wall is subjected to large and oscillatory WSS, reflected by the long and low-level tails. The statistics of WSS and its time derivative are given in Table 2. The mean values of WSS, about 2 Pa, and spatial WSS fluctuation (noted in terms of root mean square (RMS)) are similar comparing the models. The time- derivative of WSS shows larger differences, where the Walburn-Schneck model displays the largest mean and RMS values. The near wall behavior of the model is a possible explanation of the larger dissipation generated by this model during systole. The peak of the time derivative of WSS is smallest for the Casson and the Quemada models, also reflected in Table 2.
a b
c d
e f
g h
Figure 5: The probability of distribution of time-averaged WSS and TD_WSS using different rheological models: Newtonian (a and b), Casson (c and d), Quemada (e and f) and Walburn-Schneck (g and h).
Fig 6 depicts the spatial distribution of TAWSS (eq 8), OSI (eq 9) and RRT (eq 10), on the walls of the thoracic aorta. All four models show strong similarity to each other. OSI has mostly values below about 0.1 with the exception of some well localized regions. Larger OSI values (close to 0.5, indicating sign
oscillations of the WSS vector) are seen in the inner wall of the ascending aorta, in the artic sinus and near the branches from the aortic arch. RRT indicate analogous behavior with largest values at the same locations as the OSI. TAWSS displays stronger values at the junction between the aortic sinus and the ascending aorta and proximal parts of the brachiocephalic artery. In terms of the averaged parameters the differences between the models are rather small and require close examination of the non-averaged data for making a qualitative assessment. Visualizations of wall indicators have commonly been adopted in the literature. Here, a different way of assessing the differences between the models is proposed, by considering the probability distribution of the OSI, RRT and TAWSS. These parameters are time-averaged and hence the probability of distribution, as shown in Fig 7, resembles the procedure used to present the WSS and its time-derivative.
a) Newtonian model
b) Casson
c) Quemada
d) Walburn-Schneck
Figure 6: WSS indicators: OSI, RRT and TAWSS projected on the surface of the thoracic aorta for a) Newtonian; b) Casson; c) Quemada and d) Walburn-Schneck models.
OSI probability is largest for low OSI, where the distribution of probability clearly depends on the
rheological model. Lower OSI values indicate the change of sign in WSS. The Newtonian and the Walburn- Schneck models have a larger probability for OSI < 0.1 as compared to the other two models. Thus, the Newtonian and Walburn-Schneck models have larger rate of sign change in the computed WSS. For larger OSI (i.e. > 0.3), the probability level is similar for the different models, an effect occurring in surface regions without retrograde flow. The TAWSS probabilities are insensitive to the viscosity model used, both in terms of values and distribution, with the exception of the irregularity at 1 Pa. The reason for this behavior was not identified. The RRT distribution for the Casson model is narrower around the peak value (about 0.0012 s), resulting in a larger probability for greater RRT values.
Table 3 provides quantitative statistical data for the three WSS-indicators. Despite the probability distribution variations in Fig 7, the mean and RMS values of OSI are close to each other for the different models. This observation demonstrates the averaging effect of the statistical data. As the probability distributions of TAWSS and RRT do not differ strongly, the mean and RMS values in Table 3 are close to each other.
Model OSI [-] TAWSS [Pa] RRT [1/s]
Newtonian 0.112 ± 0.100 3.146 ± 4.889 0.0027 ± 0.0049
Quemada 0.105 ± 0.100 3.023 ± 3.975 0.0020 ± 0.0033
Casson 0.113 ± 0.100 3.031 ± 4.443 0.0026 ± 0.0049
Walburn-Schneck 0.112 ± 0.100 3.032 ± 4.667 0.0026 ± 0.0046 Quemada 10.4LPM 0.125 ± 0.106 3.170 ± 3.751 0.0033 ± 0.0060 Table 3: The space-time averages of the WSS-indicators
a) Newtonian
b) Casson
c) Quemada
d) Walburn-Schneck
Figures 7: Probability distribution of OSI (left), TAWSS (middle) and RRT (right) for a) Newtonian; b) Casson; c) Quemada and d) Walburn-Schneck models.
Probability
Discussion
Studies related to the importance and significance of different rheological models for blood flow simulations have been discussed extensively in the literature: The conclusions have not always been in agreement in terms of the importance of using non-Newtonian models versus using a constant viscosity value. Several authors concluded that the flow was significantly affected by certain rheological models and raised caution in using those (Jahangiri et al (2018)). Other studies reported that the computed flow field general
characteristics exhibit no major differences when using Newtonian and non-Newtonian models (Akherat et al (2017)), and additional publications find the Newtonian approximation reasonable for high shear flows (Skiadopoulos et al (2017)). The significance and importance of the models were in other studies determined by visual judgment of graphical data (e.g. surface values of OSI, TAWSS and RRT, or velocity distribution at certain cross-sections). In addition to graphical data, the study presented herein also suggests
quantification of the viscosity variation in space and time and supplemental (WSS-related) variables relevant to the formation of atherosclerosis.
The direct comparison of the relative viscosity and the RMS deviation from a reference (Newtonian case) value shows that the Walburn-Schneck model yields largest viscosity values throughout the cardiac cycle, and in particular during flow acceleration with lowest values during flow deceleration. The latter effect is related directly to the formation of retrograde flow in a substantial part of the lumen. The importance factor (Ig) reaches lowest values at local peaks in flow rate (i.e. peak systole and a local peak at early diastole, Figs 1c and d. The probability distribution of the time-averaged IL and Ig show clear differences between the non- Newtonian viscosity models. It also shows that the Casson model has least increase in viscosity (near the peak with smaller range) followed the Quemada model. The largest increase in viscosity as compared to the reference value is noted for the Walburn-Schneck model. The same effect is observed in terms of the RMS of the distributions. The impact of the models on the temporal development of the space averaged kinetic energy and the corresponding viscous dissipation rate show small differences among the models. Yet, a closer look at details such as the distribution of the kinetic energy, vorticity in the aorta and wall related parameters, WSS and its temporal derivative, OSI, TAWSS and RRT better expose the behavior of the models. Regarding WSS and its temporal derivative or the WSS-related indicators (Tables 2 and 3), the differences are modest. Tables 2 and 3 also show that increasing the flow rate to 10.4 LPM has larger impact on most parameters, as compared to the differences due to viscosity model used.
It should be emphasized that the aortic flow is characterized by large regions with relatively high shear-rates, mostly greater than order of 10 s-1. For large shear-rates (Fig 1b), the differences among the models is considerably smaller than for shear-rates below O(100 s-1). The sensitivity of the models (i.e. the derivative of viscosity with respect to shear-rate), is less than 0.01 Pa s/Pa for the Quemada model for shear-rate of 100 s-1 to 1000 s-1while it is an order of magnitude greater for the Walburn-Schneck model with a shear-rate of about 100 s-1 (not shown here). The increase in viscosity due to the non-Newtonian models is of order one.
Yet, vorticity generation within the aorta is dominated by vortex stretching and with only marginal contribution of viscosity. WSS, on the other hand is more closely dependent on the local value of the viscosity. Hence, the viscosity model could be expected to be most important for medium and smaller size arteries where the flow rate and thereby shear-rate is relatively high. This effect is clearly seen for arteries often subject to atherosclerosis; such as the coronary, carotid and renal arteries as shown by van Wyk et al (2013) for a generic aortic bifurcation. For smaller arterioles, the shear-rate is low and flow losses are viscosity dominated and hence the impact of the rheological models could be more important than for large arteries. Other situations where viscosity plays a central role is in computing viscous forces (drag and lift) acting on particles transported by blood (i.e. lipoproteins and cells). Viscous forces have the potential affecting the path of the substances, making predicting for example drug delivery a challenging problem, requiring also advance rheological models.
The results indicate that the choice of rheological model approach should depend on the problem itself (inertia or viscosity dominated flow with high or low shear-rate). Another criterion for choosing a rheological model will depend on the parameters of interest and their dependence on viscosity. In case of blood flow in larger arteries (flows with high shear-rate (mean > 100 s-1)) and with small regions of
stagnation, the Quemada model is appropriate for the range of shear-rates under consideration. The reason is due to the characteristics of the model capturing the variability in viscosity characteristics also at the low end of shear rates. For predominantly high shear-rate flows, the Newtonian approach is a natural choice. In case
of a flow containing larger regions of low shear, non-Newtonian models calibrated for appropriate ranges of shear-rates are mandatory. It should be noted though that a major shortcoming of existing models is due to being calibrated against stationary/steady flow situations. For flows with strong temporal variation and in particular in combination low shear-rate, existing models may be inadequate. This fact is important as atherosclerosis is believed to be dependent on low level but fluctuating WSS.
The limitation of the results presented here is mainly due to the fact that only one geometry with rigid walls was considered. Furthermore, the models used do not include the effects of varying hematocrit. However, the sensitivity of the three non-Newtonian models to hematocrit variations was found to be smaller than the sensitivity to shear-rate (results not shown here). Modeling of RBC transport requires further models for RBC diffusivity (for example the models of (Zydney et al. 1991), Leighton and Acrivos (1987), Phillips et al.
(1992)) implying additional modeling parameters with own uncertainties.
References:
Agarwal R, Sarkar A, Paul S, and Chakraborty S - A portable rotating disc as blood rheometer.
Biomicrofluidics 13, 064120 (2019); doi: 10.1063/1.5128937
Akherat ASMJM, Cassel K, Boghosian M, Dhar P, and Hammes M, - Are Non-Newtonian Effects Important in Hemodynamic Simulations of Patients With Autogenous Fistula? ASME Journal of Biomechanical Engineering 2017, 139 / 044504-1. DOI: 10.1115/1.4035915
Benim AC, A. Nahavandi, A. Assmann, D. Schubert, P. Feindt, S.H. Suh - Simulation of blood flow in human aorta with emphasis on outlet boundary conditions. Applied Mathematical Modelling 35 (2011) 3175–3188
Boyd J, Buick JM and Green S - Analysis of the Casson and Carreau-Yasuda non-Newtonian blood models in steady and oscillatory flows using the lattice Boltzmann method. Phys Fluids 19, 093103 (2007).
Brust, M. & Schaefer, C. & Pan, L. & Garcia, Mike & Arratia, Paulo & Wagner, C.. (2013). Rheology of Human Blood Plasma: Viscoelastic Versus Newtonian Behavior. Physical Review Letters. 110. 078305.
10.1103/PhysRevLett.110.078305.
Carreau, P.J., 1972. Rheological equations from molecular network theories. Transactions of the Society of Rheology, 16(1), 99-127.
Casson, N., 1959. Rheology of disperse systems. Pergamon Press, London, 84.
Charm, S. E. and Kurland, G. S., 1972. Blood rheology. In Cardiooascular Fluid Dynamics, Vol. 2 (Edited by Bergel)
Chen X, Gao Y, Lu B, Jia X, Zhong L, Kassab GS, et al. Hemodynamics in Coronary Arterial Tree of Serial Stenoses. PLoS ONE, 2016, 11(9): e0163715. doi:10.1371/journal. pone.0163715
Cokelet, G., Merrill, E., Gilliland, E., Shin, H., Britten, A. & Wells Jr, R. 1963 The rheology of human blood{measurement near and at zero shear rate. Journal of Rheology 7, 303{307
Cole JS, Watterson JK and O’Reilly MJG, Numerical investigation of the hemodynamics at a patched arterial bypass anastomosis, Medical Engineering & Physics 24 (2002), 393–401. Doi: 10.1016/S1350- 4533(02)00038-3
Cowan, Aimee, Cho, Daniel, Rosenson, Robert. Importance of blood rheology in the pathophysiology of atherothrombosis. Cardiovascular drugs and therapy, 26(4), 339-348. Cardiovascular Drugs and Therapy, 2012, Vol.26(4), pp.339-348
Danese E, Lippi G, Montagnana M. Red blood cell distribution width and cardiovascular diseases. J Thorac Dis. 2015;7(10):E402-11. doi: 10.3978/j.issn.2072-1439.2015.10.04
Fuchs A, Berg N, and Prahl Wittberg, L. (2019). Stenosis Indicators Applied to Patient-Specific Renal Arteries without and with Stenosis. Fluids. 4.26. DOI:10.3390/fluids4010026.
Gallagher MT, Wain RAJ, Dari S, Whitty JP and Smith DJ - Non-identifiability of parameters for a class of shear-thinning rheological models, with implications for haematological fluid dynamics. Journal of
Biomechanics 85 (2019) 230–238. Doi: 10.1016/j.jbiomech.2019.01.036
Gallo D, Steinman DA and Morbiducci U - Insights into the co-localization of magnitude-based versus direction-based indicators of disturbed shear at the carotid bifurcation. J Biomech 49 (2016) 2413–2419 Hund SJ, Kameneva MV and Antaki JF - A Quasi-Mechanistic Mathematical Representation for Blood Viscosity. Fluids, 2, 10; 2017. doi:10.3390/fluids2010010
Jahangiri, M, Haghani A, Ghaderi R, and Harat SMS, - Effect of Non-Newtonian Models on Blood Flow in Artery with Different Consecutive Stenosis. Int J Adv Design and Manufacturing Technology (ADMT). 11.
89-96. 2018
Johnston BM, Johnston PR, Corney S, Kilpatrick D - Non-Newtonian blood flow in human right coronary arteries: steady state simulations, J. Biomech. 37 (2004) 709–720
Karimi S, Dabagh M, Vasava P, Dadvar M, Dabir B, and Jalali P. - Effect of rheological models on the hemodynamics within human aorta: CFD study on CT image-based geometry. J. Nonnewton. Fluid Mech. 207:42–52, 2014.
Koh Y-M – Vorticity and Viscous Dissipation in an Incompressible Flow. KSME J., 8(1); pp.35-42,1994 LaDisa JF Jr, Dholakia RJ, Figueroa CA, Vignon-Clementel IE, Chan FP, Samyn MM, et al. Computational simulations demonstrate altered wall shear stress in aortic coarctation patients treated by resection with end- to-end anastomosis. Congenit Heart Dis. 2011, 6(5),432-443. doi: 10.1111/j.1747-0803.2011.00553.x Leighton D and Acrivos A - The shear-induced migration of particles in concentrated suspensions. J. Fluid Mech. (1987), vol. 181, p~ 415-439
Lippi G, Bovo C, Buonocore R, Mitaritonno M, Cervellin G. Red blood cell distribution width in patients with limb, chest and head trauma. Arch Med Sci. 2017;13(3):606-611.
Lowe GDO, McArdle BM, Stromberg P, Lorimer AR, Forbes CD, Prentice CRM: Increased blood viscosity and fibrinolytic inhibitor in type II hyperlipoproteinaemia. Lancet 1982;1:472-475
Mendieta JB, Fontanarosa D, Wang J, et al. The importance of blood rheology in patient-specific
computational fluid dynamics simulation of stenotic carotid arteries [published online ahead of print, 2020 Jan 2]. Biomech Model Mechanobiol. 2020;10.1007/s10237-019-01282-7.
Molla MM, and Paul MC. – LES of non-Newtonian physiological blood flow in a model of arterial stenosis.
Med Eng and Physics 34 (2012) 1079–1087
Perktold K, Karner G, Leuprecht A, Hofer M, Influence of non-Newtonian flow behavior on local hemodynamics, Z Angew Math Mech 79(S1):187–190, 1999, doi: 10.1002/zamm.19990791349.
Phillips RJ, Armstrong RC, Brown RA, Graham AL, and Abbott JR - A constitutive equation for concentrated suspensions that accounts for shear induced particle migration. Physics of Fluids A: Fluid Dynamics 4, 30 (1992); doi: 10.1063/1.858498
Quemada D. 1977 Rheology of concentrated disperse systems and minimum energy dissipation principle.
Rheologica Acta 16 (1), 82{94.
Quemada D. 1978 Rheology of concentrated disperse systems ii. a model for non-newtonian shear viscosity in steady ows. Rheologica Acta 17 (6), 632-642.
Rikhtegar F, Knight JA, Olgac U, Saur SC, Poulikakos D, Marshall W, Cattin PC, Alkadhi H, Kurtcuoglu V - Choosing the optimal wall shear parameter for the prediction of plaque location–a patient-specific
computational study in human left coronary arteries, Atherosclerosis 221 (2012) 432–437. doi:
10.1016/j.atherosclerosis.2012.01.018
Skiadopoulos A, Neofytou P, and Housiadas Ch. - Comparison of blood rheological models in patient specific cardiovascular system simulations. J hydrodynamics, 2017,29(2):293-304 DOI: 10.1016/S1001- 6058(16)60739-4
Suo J, Oshinski JN, Giddens DP. Blood flow patterns in the proximal human coronary arteries: relationship to atherosclerotic plaque occurrence. Mol Cell Biomech. 2008;5:9–18
Suzuki T et al., Variability of hemodynamic parameters using the common viscosity assumption in a computational fluid dynamics analysis of intracranial aneurysms, Technol Health Care 15(1):37–47, 2017.
van Wyk S, Prahl Wittberg L, Fuchs L, - Wall shear stress variations and unsteadiness of pulsatile blood-like flows in 90-degree bifurcations. Computers in Biology and Medicine, 43(8), 2013, 1025-1036,
https://doi.org/10.1016/j.compbiomed.2013.05.008.
Walburn, F. & Schneck, D. 1976 A constitutive equation for whole human blood. Biorheology, 1976; 13(3), 201-210.
Yamamoto H, Kawamura K, Omura K, Tokudome S. Development of a compact-sized falling needle rheometer for measurement of flow properties of fresh human blood. Int J Thermophys. 2010; 31(11): 2361- 79, DOI 10.1007/s10765-010-0842-5.
Zarins CK, Giddens DP, Bharadvaj BK, Sottiurai VS, Mabon RF, Glagov S. Carotid bifurcation
atherosclerosis. Quantitative correlation of plaque localization with flow velocity profiles and wall shear stress. Circ Res. 1983, 53(4), 502-514. doi: 10.1161/01.RES.53.4.502
Zydney AL, Oliver JD, Colton CK. - A constitutive equation for the viscosity of stored red cell suspensions:
Effect of hematocrit, shear rate, and suspending phase. J Rheology, 1991, 35(8), pp.1639-1680
