Comparison of Non-Newtonian Models of One-Dimensional Hemodynamics

The paper is devoted to the comparison of different one-dimensional models of blood flow. In such models, the non-Newtonian property of blood is considered. It is demonstrated that for the large arteries, the small parameter is observed in the models, and the perturbation method can be used for the analytical solution. In the paper, the simplified nonlinear problem for the semi-infinite vessel with constant properties is solved analytically, and the solutions for different models are compared. The effects of the flattening of the velocity profile and hematocrit value on the deviation from the Newtonian model are investigated.


Introduction
The one-dimensional (1D) models are used for the simulation of blood flow in large vascular systems, where the application of 2D and 3D models leads to computational difficulties [1][2][3][4]. In clinical practice, such models are used as a predictive tool in vascular surgery [5]. The 1D models are constructed by the averaging of the equations of hydrodynamics on the vessel cross-section [6,7]. The results obtained by 1D models are compatible with the averaged results of the 3D simulations [8].
As it is obtained in experiments, blood demonstrates the non-Newtonian behaviour [9], with the shear-thinning [10], thixotropic [11] and viscoelastic [12] properties. The complete 1D model of blood, with time-dependent characteristics, where all these properties are taken into account, is presented in [13]. As it is mentioned in [14], only the shear-thinning property is important for the blood flow models because the other mentioned properties do not significantly change the velocity field. The shear-thinning property of blood can be successively described by the generalized Newtonian models, where the dynamic viscosity is dependent only on the shear rate tensor [15]. For example, the models of such type are applied to the 2D and 3D blood flow simulations in [16][17][18][19][20][21].
The typical example of the generalized Newtonian fluid is a so-called Power Law model. This model is used for the 2D and 3D modelling of blood flow dynamics in a coronary artery [22], aorta [17], aorta-iliac bifurcation [23], thoracic aorta [24], carotid artery [25] and in various types of stenosed vessels [18,20,26]. In [22,27], it is demonstrated that at high shear rates, the solutions are close to the solutions for the Newtonian model, but at low shear rates, serious deviations are observed.
In many works, models with two asymptotic values of dynamic viscosity, corresponding to infinitely small and infinitely large shear rates, are used. The examples of such models are: Carreau model [9,[16][17][18][19][20]22,24,[27][28][29][30][31]; Carreau-Yasuda model [9,17,18,22,23,25,27,[31][32][33]; Cross model [9,16,17,20,27,34,35]; Simplified Cross model [9,16,36]; Modified Cross model [9,16,22,27,36]; and Yeleswarapu model [37][38][39][40]. In [16], models of such types are applied to the simulation of blood flow in atherosclerotic coronary arteries, and results are compared with clinical measurements. The lowest average deviations are obtained for Cross and Simplified Cross models. The largest deviations are observed for the Modified Cross model. In [23], the non-Newtonian effects on the low-density lipoprotein transport The model of blood flow as a flow of viscous incompressible fluid is described by the incompressibility condition and momentum equation, when the body force action is ignored: where ρ is a constant density, t is time, V is the velocity, p is the pressure, T is the viscous stress tensor. The time-independent non-Newtonian model (so-called generalized Newtonian fluid) is based on the following relation [15]: where D is the strain rate tensor, I 2 is its second invariant and µ is a dynamic viscosity. The case of constant value of µ corresponds to the Newtonian fluid. The expressions for µ(I 2 ) for the models, widely-used for the blood, are presented in Table 1. The value of parameter µ for the Newtonian model can be taken from [16,18,21,28]; values of n and k for the Power Law model can be taken from [23,46,47]; values of µ 0 , µ ∞ , λ and n for the Carreau model can be taken from [9,[16][17][18][19][20][21][22][23][24]48]; values of µ 0 , µ ∞ , λ, n and a for the Carreau-Yasuda model can be taken from [9,16,23,33]; values of µ 0 , µ ∞ , λ and m for the Cross model can be taken from [9,16,17,20,35]; values of µ 0 , µ ∞ and λ for the Simplified Cross model can be taken from [9]; values of µ 0 , µ ∞ , λ, m and a for the Modified Cross model can be taken from [9]; values of µ 0 , µ ∞ and λ for the Yeleswarapu model can be taken from [16,38]; the dependence of these parameters on hematocrit H for the Modified Yeleswarapu model are presented in [21] and values of µ, k 0 , k ∞ andγ c for the Quemada model are presented in [16,20,49]. Table 1. Rheological models of blood as generalized Newtonian fluid.

Model µ(I 2 )
Newtonian The system (1) under some assumptions (e.g., see [2,6]) can be reduced to the following form: where (r, z) are the cylindrical coordinates, P = P(t, z) is the pressure, corresponding to the vessel cross-section at coordinate z, V r and V z are the components of V. System (2) can be averaged on the vessel cross-section with the following representation of V z [6]: where U(t, z) is the mean velocity, R(t, z) is the vessel radius and s is the dimensionless velocity profile.
As it is demonstrated in [6,7], system (2) after the averaging is reduced to the following form: ∂A ∂t where A = A(t, z) is the cross-sectional area, Q = Q(t, z) is the flow rate (Q = AU), f is the frictional term, defined as [13]: and α is the momentum correction coefficient, obtained as: where σ is the vessel cross-section. System (4) must be closed by the equation-of-state P = P(A). For the case of flow in arteries, the following equation is widely used [8]: where P ext is the external pressure, A d and P d are the diastolic cross-sectional area and pressure, β = 4 3 √ πEh, where E is the Young's modulus and h is the vessel wall thickness. The case of constant values of these parameters is considered in the present paper.
After the substitution of (5) into (4), the second equation is rewritten as:

∂Q ∂t
The system (4) can be rewritten using the new dimensionless variables: where L C is considered as a vessel length and U C is estimated by the wave speed at mean pressure, T C = L C /U C . The tilde sign will be ignored in the text below.
In the dimensionless variables, system (4) is rewritten as:

∂A ∂t
where The expressions of f (A, Q) for the models from Table 1 in the dimensionless case are presented in Table 2. In these expressions, the values of s (1) are used. The values of α are computed from the values of s(y), where y = r/R. Thus, the system (6) will be closed if the dependence of s on y will be introduced. This expression can be obtained from the solution of the problem for the steady flow in a circular tube. For the Power Law model, this problem is solved exactly: and α is computed as: In the cases of other models, s(y) can be obtained only numerically. It can be dependent on the flow conditions (e.g., see [30,31]). According to the existence of the cellular part, the flattened velocity profile is typical for the blood [50]. Such profiles are obtained in 2D and 3D simulations of blood flow, based on the non-Newtonian models in [9,18,31,32]. For the introduction of the flattened profile by (3), into the 1D Newtonian model, in [1,4,6,7] the following expression is used: where d is the dimensionless parameter. Plots of s(y) at various values of d are presented in Figure 1. As it can be seen, the flattening of the profile can be regulated by the proper value of d.
A similar approach is used in [13] for the time-dependent non-Newtonian 1D model, but s is dependent on t and y. The expression for the velocity profile is represented by the Womersley solution for the pulsatile flow of the Newtonian fluid in a circular tube. Table 2. Expressions for the frictional term.

Model f (A, Q)
Newtonian Being motivated by the works mentioned above, in the presented paper we decided to use the expression (7) for the velocity profile. It must be noted that the use of the prescribed expression for the profile is a lack of the considered 1D models. However, as it is mentioned above, it is used in many other works. Moreover, as it is demonstrated in [8], the effects of a more realistic representation of the profile, where effects such as the existence of boundary layers, recirculation zones and others that can be taken into account, do not seriously change the values of the averaged characteristics. As a result, the model representation (7) is used as an approximation or simplification in this study, but parameter d will be varied to estimate the effect of profile flattening on the solution.
The value of α, corresponding to (7), is computed as:

Dimensionless Parameters
For the estimation of parameters for dimensionless Equation (6), the values of L C , U C and A d are taken from [8] for the large arteries, such as a carotid artery, aorta and iliac artery. The value of ρ is taken as 1.05 g/cm 3 . The value of χ is estimated as 0.9783 for the carotid artery, 0.9501 for the aorta and 0.9739 for the iliac artery. Parameters ε, ξ and ζ are dependent on the parameters of rheological models and vessel properties. In the presented study, it is proposed that the following relationship between ε and ξ takes place: For the Newtonian model, ε is approximately equal to 0.0137 for the carotid artery and ≈0.0020 for the aorta and iliac artery. For the Power Law model and parameters from [46] ε ≈ 0.0063 for the carotid artery, 0.0010 for aorta and 9.9 × 10 −4 for the iliac artery. For the parameters from [23,47], ε ≈ 0.0096 for carotid artery, 0.0018 for aorta and 0.0015 for iliac artery.
For the Carreau-Yasuda model and parameters from [9,16,23,33], C is approximately equal to 15. The parameter ε has the same values as for the Carreau model. Parameter ζ has the following values: 3.4552 × 10 4 for the carotid artery; 4.1662 × 10 3 for aorta and 1.6325 × 10 4 for iliac artery. For the Cross model and parameters from [9], the value of C is estimated as 15. Parameter ε has the same values as for Carreau and Carreau-Yasuda model. For ζ, the following values take place: 2.8088 × 10 3 for the carotid artery, 4.9312 × 10 2 for aorta and 0.0084, 1.5161 × 10 3 for the iliac artery. For the parameters from [16,17,20,35] C ≈ 9.55. Parameter ε has the same values. Parameter ζ is equal to 1.7801 × 10 4 for the carotid artery, 1.5301 × 10 3 for aorta and 7.4602 × 10 3 for the iliac artery.
For the case of the Modified Cross model C is estimated as 25 and ε has the same values as for the Cross model. For the carotid artery, ζ ≈ 2.76116896 × 10 9 ; 4.70660529 × 10 7 for aorta and 6.52156883 × 10 8 for iliac artery.
As can be seen, the small parameter ε is observed in the system (6), so the perturbation method can be applied to obtain analytical solutions. For most of the considered models, the large value of ζ takes place. The first term −εQ/A in the expression for f (A, Q) prevails, and if the values of ε are close, it can lead to the close values of the solutions for different models.

Solution of the Problem for the Semi-Infinite Interval
In this section, for the comparison of the models, the problem for the system (6) in the semi-infinite spatial interval is solved analytically by the perturbation method. The problem is considered as a model of the flow in a semi-infinite vessel with constant mechanical properties, which is situated at the interval z ∈ [0, + ∞). At the left side of this interval, the flow is induced by the periodic time functions, which are presented as the value of the flow rate. The perturbations induced by this function propagate at z → + ∞. Thus, we try to compare models on the problem, which can be solved analytically. For the solution, the assumption Q > 0 is used. The values of dimensionless parameters are taken for the case of the carotid artery.
It must be noted that the considered problem is analytically solved for the following purposes: (1) To compare the results for different non-Newtonian models and estimate the deviation from the Newtonian solution; (2) To estimate the effect of the velocity profile; (3) To estimate the effect of hematocrit H (for the Quemada and modified Yeleswarapu models); (4) To provide a tool for testing the programs that implement the numerical methods algorithms.

Perturbation Method
Let the semi-infinite interval z ∈ [0, + ∞) be considered. The initial conditions are presented as: where ϕ i (z),ψ i (z) are bounded at z → ±∞, A 0 ,Q 0 are the constants. The solutions obtained by this method correspond to the situation when the steady state, defined by the constant values A 0 and Q 0 , is perturbed by the small additive terms, as shown in (8) and (9). The boundary condition for (6) is presented as: where ω i (t) are discussed below. This condition induces the flow at the left boundary z = 0. The solution of problem (6), (8)-(10) is presented in the following form: After the substitution of (11), (12) into (6), the problems for A i (t, z), Q i (t, z) are formulated. In the presented study, the terms only up to and including the second order of ε are considered. It is easy to obtain that A 0 (t, z) For the first-order terms A 1 and Q 1 , the following problem takes place: where b 1 , b 2 and f 0 are defined as: For the second-order terms A 2 , Q 2 , the following problem is considered: where and constants δ 1 , δ 2 for every model have its own expressions. For the Power Law model (including Newtonian at n = 1), they are written as: For the Carreau-Yasuda model (including Carreau at a = 2), the following expressions take place: for the Modified Cross model (a = 1 corresponds to the Cross model, a = 1, m = 1-to the Simplified Cross model): For the Yeleswarapu model and its modified form, δ 1 and δ 2 are presented as: For the Quemada model, the following expressions are obtained: Let the initial functions in (14) and (17) be presented as: where i 2 = −1 and A, Q, A, Q are the real constants.
For the simulation of the oscillatory behaviour of flow rate at z = 0, it is proposed that function ω 1 (t) in (15) is presented as: where q 0 (t) and q 1 (t) are the solutions of the following initial problems: It is easy to obtain that: where γ 1,2 = i 2 −b 2 ± 4b 1 + b 2 2 and at α ≥ 1. Constants C 1 and C 2 are computed as: The solution of the problem (13)-(15) is written as: Let the function ω 2 (t) be presented as: where q 0 , q 1 , q 2 are the solutions of the following initial problems: where g 1 (t) = δ 1 a 1 (t) + δ 2 q 1 (t) + iβ 1 q 0 (t)a 1 (t) + 2iβ 2 q 0 (t)q 1 (t), where g 2 (t) = 2iβ 1 a 1 (t)q 1 (t) + 2iβ 2 q 2 1 (t) + 2iβ 3 a 2 1 (t). The initial problems (19) and (20) can be rewritten as: where P is a constant complex matrix. The solution of (21) is obtained as: where Y(t) is the fundamental matrix of systemẏ = Py. The expression for q 0 (t) is written as: q 0 (t) = δ 2 f 0 t 2 2 , expressions for a 1 (t), q 1 (t), a 2 (t), q 2 (t) are not presented according to its cumbersomeness, but they can be easily obtained using the systems of symbolic computations.

Comparison of Solutions
For the comparison of the solutions, obtained for the non-Newtonian and Newtonian models, the following criterion, which is called the non-Newtonian factor, is used: where L 2 norm is used, Q Newt is the solution for the Newtonian model. The values of Re(Q) are considered in NNF because the damping of the solution, induced by the viscosity, is most evident for this function. For various models, the following factors affect the obtained solutions: the value of ε and the expression for the nonlinear frictional term f (A, Q), which are defined by the characteristics of rheological models. These factors can lead to the difference between the solutions. For the comparison of the solutions, the following parameter values are used: 15]. The plots of the real part of ω(t) and corresponding Re(A(t, 0)) for the Power Law model with parameters from [47] are presented in Figure 2. Figure 3 shows typical plots of the real parts of Q and A.
The value of NNF for the Power Law model is equal to 2.5826 % for the parameters from [47] and to 4.4378 % for parameters from [46]. The plots of the Re(Q) at the midpoint of the considered space interval are presented in Figure 4. As can be seen, the strongest damping is demonstrated for the Newtonian model.
The plots of NNF for the models without the hematocrit at different values of d are shown in Figure 5. The maximum values take place for Simplified Cross and Yeleswarapu models (≈14%) at d = 6. The maximum values for other models are realized also at this value of d and are approximately equal to 8 %. For the parabolic profile (d = 2), typical for the Newtonian model, the NNF is equal to 3 % for the Simplified Cross and Yeleswarapu model and can be considered as negligible for the other models.
In Figure 6, the plots of NNF for the models, where the value of H is included, are presented. At most H values, the value of NNF for the parabolic profile is greater than for the Modified Yeleswarapu model, but at other d values, the situation is reversed. The greatest deviations from the solution for the Newtonian model are observed at d = 6 (approximately 11% for the Quemada model and 17% for the modified Yeleswarapu model). The plots of Re(Q) at selected values of d and H are presented in Figure 7 in comparison with the solution for the Newtonian case.

Conclusions
In the presented paper, the following problems are considered: (1) The non-Newtonian models of 1D hemodynamics are presented. Such models are constructed by the averaging of the equation of incompressible flow of viscous fluid over a vessel's cross-section. The models are characterized by the nonlinear frictional terms, obtained from the rheological relations on the stress tensor.
(2) It is demonstrated that in the case of large arteries, small parameters are observed in the models, so the perturbation method can be used for the analytical solution of model problems.
(3) The problem for the semi-infinite space interval is solved. The expressions for the first-and second-order terms in the expansion on the small parameter are obtained.
(4) For the comparison of the solutions, the non-Newtonian factor NNF is introduced. For the Power Law model, it is demonstrated that the flattening of the velocity profile does not lead to more significant damping of the solution than for the Newtonian case. For the other models, the deviation from the Newtonian case starts to increase with the increase in the flattening of the profile. For the models, where the hematocrit is taken into account, the deviation from the Newtonian solution increases with the increase of hematocrit value at all considered velocity profiles.
It must be noted that the obtained analytical solutions can be used for the testing of the programs, which implement the algorithms for the simulation of blood flow in large vascular systems.
The comparison of models is based on a simplified problem which can be solved analytically. For a whole comparison, the models will be compared in the cases of flows in large vascular systems, which will be considered in future works.

Acknowledgments:
The author wishes to thank the anonymous referees for careful checking of the article and helpful comments.

Conflicts of Interest:
The author declares no conflict of interest.