Comparison of Non-Newtonian Models of One-Dimensional Hemodynamics

Abstract

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.

1. 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 in arteries are analysed. It is demonstrated that the Newtonian model is valid for mass transport at low Reynolds numbers. At high values, the non-Newtonian models provide more accurate results.

Some non-Newtonian models are dependent on the hematocrit—the typical examples are: the Quemada model, used in [16,20,27,41,42], and the Modified Yeleswarapu model [21,43].

In most works devoted to 1D models of blood flow, the non-Newtonian properties are ignored. As it is mentioned above, the complete 1D model, where the time-dependent non-Newtonian properties are considered, is proposed by Ghigo et al. [13]. Some attempts to construct the non-Newtonian models of blood flow in large vascular networks have been made in [44,45].

In the present paper, the non-Newtonian 1D models of blood flow, obtained by the averaging of Navier–Stokes equations on the vessel cross-section, are presented. It is demonstrated that for the large arteries, the small parameters take place in the models. It provides the possibility to use the perturbation method for obtaining analytical solutions of the problems for a 1D system of hemodynamical equations. The simplified but nonlinear problem for the case of a semi-infinite vessel is solved analytically. Thus, in the paper, the theoretical problem for the comparison of different 1D models is considered. The solutions obtained for the non-Newtonian models are compared with each other. The effects of the flattening of the velocity profile and hematocrit are analysed.

The paper has the following structure. In Section 2, the 1D models of blood flow are considered. The nonlinear expressions for the frictional term are presented. Section 3 is devoted to the application of the perturbation method to obtain the analytical solution of the problem for the semi-infinite interval. The solutions for the different models are compared by the values of the non-Newtonian factor, which estimates the deviations of the solutions for Newtonian and non-Newtonian models. Some concluding remarks are made in Section 4.

2. Non-Newtonian Models of One-Dimensional Hemodynamics

2.1. Non-Newtonian Models of Blood Flow

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:

$$\nabla \times \mathbf{V}=0,\rho \left(\frac{\partial \mathbf{V}}{\partial t}+\mathbf{V}\times \nabla \mathbf{V}\right)=-\nabla p+\nabla \times \mathbf{T},$$

(1)

where $\rho$ is a constant density, t is time, $\mathbf{V}$ is the velocity, p is the pressure, $\mathbf{T}$ is the viscous stress tensor.

The time-independent non-Newtonian model (so-called generalized Newtonian fluid) is based on the following relation [15]:

$$\mathbf{T}=2\mu \left(I_2\right)\mathbf{D},$$

where $\mathbf{D}$ is the strain rate tensor, $I_2$ is its second invariant and $\mu$ is a dynamic viscosity. The case of constant value of $\mu$ corresponds to the Newtonian fluid. The expressions for $\mu \left(I_2\right)$ for the models, widely-used for the blood, are presented in

Table 1. Rheological models of blood as generalized Newtonian fluid.

Model	$\mathit{\mu }\left({\mathit{I}}_2\right)$
Newtonian	$\mu \left(I_2\right) = \mu = const$
Power Law	$\mu \left(I_2\right) = k{I}_2^{\frac{n - 1}{2}}$
Carreau	$\mu \left(I_2\right) = {\mu }_{\infty } + ({\mu }_0 - {\mu }_{\infty }){\left(1 + {\left(\lambda \sqrt{I_2}\right)}^2\right)}^{\frac{n - 1}{2}}$
Carreau–Yasuda	$\mu \left(I_2\right) = {\mu }_{\infty } + ({\mu }_0 - {\mu }_{\infty }){\left(1 + {\left(\lambda \sqrt{I_2}\right)}^a\right)}^{\frac{n - 1}{a}}$
Cross	$\mu \left(I_2\right) = {\mu }_{\infty } + \frac{{\mu }_0-{\mu }_{\infty }}{1 + {\left(\lambda \sqrt{I_2}\right)}^m}$
Simplified Cross	$\mu \left(I_2\right) = {\mu }_{\infty } + \frac{{\mu }_0 - {\mu }_{\infty }}{1 + \lambda \sqrt{I_2}}$
Modified Cross	$\mu \left(I_2\right)={\mu }_{\infty } + \frac{{\mu }_0-{\mu }_{\infty }}{{\left(1 + {\left(\lambda \sqrt{I_2}\right)}^m\right)}^a}$
Yeleswarapu	$\mu \left(I_2\right)={\mu }_{\infty } + ({\mu }_0-{\mu }_{\infty })\frac{1 + \ln \left(1 + \lambda \sqrt{I_2}\right)}{1 + \lambda \sqrt{I_2}}$
Modified Yeleswarapu
Quemada	$\mu \left(I_2\right) = \mu {\left(1-\frac{1}{2}{k}_{Q}H\right)}^{ - 2}$, $k_Q = \frac{k_0 + k_{\infty }\sqrt{\frac{\sqrt{I_2}}{{\dot{\gamma }}_c}}}{1 + \sqrt{\frac{\sqrt{I_2}}{{\dot{\gamma }}_c}}}$

. The value of parameter $\mu$ 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 ${\mu }_0$, ${\mu }_{\infty }$, $\lambda$ and n for the Carreau model can be taken from [9,16,17,18,19,20,21,22,23,24,48]; values of ${\mu }_0$, ${\mu }_{\infty }$, $\lambda$, n and a for the Carreau–Yasuda model can be taken from [9,16,23,33]; values of ${\mu }_0$, ${\mu }_{\infty }$, $\lambda$ and m for the Cross model can be taken from [9,16,17,20,35]; values of ${\mu }_0$, ${\mu }_{\infty }$ and $\lambda$ for the Simplified Cross model can be taken from [9]; values of ${\mu }_0$, ${\mu }_{\infty }$, $\lambda$, m and a for the Modified Cross model can be taken from [9]; values of ${\mu }_0$, ${\mu }_{\infty }$ and $\lambda$ 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 $\mu$, $k_0$, $k_{\infty }$ and ${\dot{\gamma }}_c$ for the Quemada model are presented in [16,20,49].

The system (1) under some assumptions (e.g., see [2,6]) can be reduced to the following form:

$$\nabla \times \mathbf{V}=0,\frac{\partial V_z}{\partial t}+V_r\frac{\partial V_z}{\partial r}+V_z\frac{\partial V_z}{\partial z}=-\frac{1}{\rho }\frac{\partial P}{\partial z}+\frac{1}{\rho r}\frac{\partial }{\partial r}\left(r{T}_{rz}\right),$$

(2)

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 $\mathbf{V}$. System (2) can be averaged on the vessel cross-section with the following representation of $V_z$ [6]:

$$V_z(t,z)=U(t,z)s\left(\frac{r}{R(t,z)}\right),$$

(3)

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:

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial z}=0,\frac{\partial Q}{\partial t}+\frac{\partial }{\partial z}\left(\alpha \frac{Q^2}{A}\right)+\frac{A}{\rho }\frac{\partial P}{\partial z}=f(A,Q),$$

(4)

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]:

$$f(A,Q)=\frac{2\pi }{R}{T}_{rz}{|}_{r=R},$$

and $\alpha$ is the momentum correction coefficient, obtained as:

$$\alpha =\frac{\underset{\sigma }{\int }{s}^{2}dS}{A},$$

where $\sigma$ is the vessel cross-section. System (4) must be closed by the equation-of-state $P=P\left(A\right)$. For the case of flow in arteries, the following equation is widely used [8]:

$$P-P_{ext}=P_d+\frac{\beta }{A_d}(\sqrt{A}-\sqrt{A_d}),$$

(5)

where $P_{ext}$ is the external pressure, $A_d$ and $P_d$ are the diastolic cross-sectional area and pressure, $\beta =\frac{4}{3}\sqrt{\pi }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:

$$\frac{\partial Q}{\partial t}+\frac{\partial }{\partial z}\left(\alpha \frac{Q^2}{A}\right)+\gamma \sqrt{A}\frac{\partial A}{\partial z}=f(A,Q),$$

where $\gamma =\frac{\beta }{2A_d\rho }$.

The system (4) can be rewritten using the new dimensionless variables:

$$\tilde{z}=\frac{z}{L_C},\tilde{t}=\frac{t}{T_C},\tilde{A}=\frac{A}{A_d},\tilde{p}=\frac{p}{\rho U_C^2},\tilde{Q}=\frac{Q}{A_d{U}_C}.$$

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:

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial z}=0,\frac{\partial Q}{\partial t}+\frac{\partial }{\partial z}\left(\alpha \frac{Q^2}{A}\right) + \chi \sqrt{A}\frac{\partial A}{\partial z}=f(A,Q),$$

(6)

where $\chi =\gamma \sqrt{A_C}/U_C^2$.

The expressions of $f(A,Q)$ for the models from Table 1 in the dimensionless case are presented in

Table 2. Expressions for the frictional term.

Model	$\mathit{f}(\mathit{A},\mathit{Q})$
Newtonian	$f(A,Q)=s^{\prime }\left(1\right)\epsilon \frac{Q}{A}$, $\epsilon =-\frac{K{T}_C}{A_C}$, $K=-\frac{2\pi \mu }{\rho }$,
Power Law	$f(A,Q)=s^{\prime }\left(1\right){\left|s^{\prime }\left(1\right)\right|}^{n-1}\epsilon \frac{{Q\left|Q\right|}^{n-1}}{A^{\frac{3n-1}{2}}}$, $\epsilon =-\frac{K{U}_C^{n-1}{T}_C}{A_C^{\frac{n+1}{2}}}$,
	$K=-\frac{2k{\pi }^{\frac{n+1}{2}}}{\rho }$
Carreau	$f(A,Q)=s^{\prime }\left(1\right)\epsilon \frac{Q}{A}+s^{\prime }\left(1\right)\xi {\left(1+{\left(s^{\prime }\left(1\right)\right)}^2\zeta \frac{Q^2}{A^3}\right)}^{\frac{n-1}{2}}\frac{Q}{A}$,
	$\epsilon =-\frac{K_1{T}_C}{A_C}$, $\xi =-\frac{K_2{T}_C}{A_C}$, $\zeta =\frac{K_3{U}_C^2}{A_C}$,
	$K_1=-\frac{2\pi {\mu }_{\infty }}{\rho }$, $K_2=-\frac{2\pi ({\mu }_0-{\mu }_{\infty })}{\rho }$, $K_3={\lambda }^2\pi$
Carreau–Yasuda	$f(A,Q)=s^{\prime }\left(1\right)\epsilon \frac{Q}{A}+s^{\prime }\left(1\right)\xi {\left(1+|s^{\prime }{\left(1\right)|}^a\zeta \frac{{\left|Q\right|}^a}{A^{\frac{3a}{2}}}\right)}^{\frac{n-1}{a}}\frac{Q}{A}$,
	$\epsilon =-\frac{K_1{T}_C}{A_C}$, $\xi =-\frac{K_2{T}_C}{A_C}$, $\zeta =\frac{K_3{U}_C^a}{A_C^{\frac{a}{2}}}$,
	$K_1=-\frac{2\pi {\mu }_{\infty }}{\rho }$, $K_2=-\frac{2\pi ({\mu }_0-{\mu }_{\infty })}{\rho }$, $K_3={\lambda }^a{\pi }^{\frac{a}{2}}$
Cross	$f(A,Q)=s^{\prime }\left(1\right)\epsilon \frac{Q}{A}+s^{\prime }\left(1\right)\xi \frac{1}{1+|s^{\prime }{\left(1\right)|}^m\zeta \frac{{\left|Q\right|}^m}{A^{\frac{3m}{2}}}}\frac{Q}{A}$,
	$\epsilon =-\frac{K_1{T}_C}{A_C}$, $\xi =-\frac{K_2{T}_C}{A_C}$, $\zeta =\frac{K_3{U}_C^m}{A_C^{\frac{m}{2}}}$,
	$K_1=-\frac{2\pi {\mu }_{\infty }}{\rho }$, $K_2=-\frac{2\pi ({\mu }_0-{\mu }_{\infty })}{\rho }$, $K_3={\lambda }^m{\pi }^{\frac{m}{2}}$
Simplified Cross	$f(A,Q)=s^{\prime }\left(1\right)\epsilon \frac{Q}{A}+s^{\prime }\left(1\right)\xi \frac{1}{1+|s^{\prime }\left(1\right)|\zeta \frac{\left|Q\right|}{A^{\frac{3}{2}}}}\frac{Q}{A}$,
	$\epsilon =-\frac{K_1{T}_C}{A_C}$, $\xi =-\frac{K_2{T}_C}{A_C}$, $\zeta =\frac{K_3{U}_C}{\sqrt{A_C}}$,
	$K_1=-\frac{2\pi {\mu }_{\infty }}{\rho }$, $K_2=-\frac{2\pi ({\mu }_0-{\mu }_{\infty })}{\rho }$, $K_3=\lambda {\pi }^{\frac{1}{2}}$
Modified Cross	$f(A,Q)=s^{\prime }\left(1\right)\epsilon \frac{Q}{A}+s^{\prime }\left(1\right)\xi \frac{1}{{\left(1+|s^{\prime }{\left(1\right)|}^m\zeta \frac{{\left|Q\right|}^m}{A^{\frac{3m}{2}}}\right)}^a}\frac{Q}{A}$,
	$\epsilon =-\frac{K_1{T}_C}{A_C}$, $\xi =-\frac{K_2{T}_C}{A_C}$, $\zeta =\frac{K_3{U}_C^m}{A_C^{\frac{m}{2}}}$,
	$K_1=-\frac{2\pi {\mu }_{\infty }}{\rho }$, $K_2=-\frac{2\pi ({\mu }_0-{\mu }_{\infty })}{\rho }$, $K_3={\lambda }^m{\pi }^{\frac{m}{2}}$
Yeleswarapu	$f(A,Q)=s^{\prime }\left(1\right)\epsilon \frac{Q}{A}+s^{\prime }\left(1\right)\xi \frac{1+\ln \left(1+|s^{\prime }\left(1\right)|\zeta \frac{\left|Q\right|}{A^{\frac{3}{2}}}\right)}{1+|s^{\prime }\left(1\right)|\zeta \frac{\left|Q\right|}{A^{\frac{3}{2}}}}\frac{Q}{A}$,
	$\epsilon =-\frac{K_1{T}_C}{A_C}$, $\xi =-\frac{K_2{T}_C}{A_C}$, $\zeta =\frac{K_3{U}_C}{\sqrt{A_C}}$,
	$K_1=-\frac{2\pi {\mu }_{\infty }}{\rho }$, $K_2=-\frac{2\pi ({\mu }_0-{\mu }_{\infty })}{\rho }$, $K_3=\lambda \sqrt{\pi }$
Quemada	$f(A,Q)=s^{\prime }\left(1\right)\epsilon {\left(1-\frac{H}{2}\frac{k_0+k_{\infty }\sqrt{|s^{\prime }\left(1\right)|}\zeta \frac{\sqrt{\left|Q\right|}}{A^{\frac{3}{4}}}}{1+\sqrt{|s^{\prime }\left(1\right)|}\zeta \frac{\sqrt{\left|Q\right|}}{A^{\frac{3}{4}}}}\right)}^{-2}\frac{Q}{A}$
	$\epsilon =-\frac{K_1{T}_C}{A_C}$, $\zeta =\frac{K_2\sqrt{U_C}}{A_C^{\frac{1}{4}}}$
	$K_1=-\frac{2\pi \mu }{\rho }$, $K_2=\frac{{\pi }^{\frac{1}{4}}}{\sqrt{{\dot{\gamma }}_c}}$

. In these expressions, the values of $s^{\prime }\left(1\right)$ are used. The values of $\alpha$ are computed from the values of $s\left(y\right)$, 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:

$$s\left(y\right)=\frac{3n+1}{n+1}\left(1-y^{1+\frac{1}{n}}\right),$$

and $\alpha$ is computed as:

$$\alpha =\frac{3n+1}{2n+1}.$$

In the cases of other models, $s\left(y\right)$ 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:

$$s\left(y\right)=\frac{d+2}{d}\left(1-y^d\right).$$

(7)

where d is the dimensionless parameter. Plots of $s\left(y\right)$ 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.

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 $\alpha$, corresponding to (7), is computed as:

$$\alpha =\frac{d+2}{d+1}.$$

2.2. 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 $\rho$ is taken as 1.05 $\mathrm{g}/\mathrm{c}{\mathrm{m}}^3$. The value of $\chi$ is estimated as 0.9783 for the carotid artery, 0.9501 for the aorta and 0.9739 for the iliac artery. Parameters $\epsilon$, $\xi$ and $\zeta$ are dependent on the parameters of rheological models and vessel properties. In the presented study, it is proposed that the following relationship between $\epsilon$ and $\xi$ takes place: $\xi =\mathcal{C}\epsilon$, where $\mathcal{C}$ is a constant.

For the Newtonian model, $\epsilon$ 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] $\epsilon \approx 0.0063$ for the carotid artery, $0.0010$ for aorta and $9.9\times {10}^{-4}$ for the iliac artery. For the parameters from [23,47], $\epsilon \approx 0.0096$ for carotid artery, $0.0018$ for aorta and $0.0015$ for iliac artery.

For the Carreau model and parameters from [9,16,17,18,20,21,22,23,24,48], $\mathcal{C}$ is approximately equal to 15. For the dimensionless parameters, the following values take place: $\epsilon \approx 0.0135$, $\zeta \approx 5.54098186\times {10}^7$ for carotid artery; $0.0019$ and $1.87761843\times {10}^6$ for aorta; $0.0020$, $1.66957179\times {10}^7$ for the iliac artery.

For the Carreau–Yasuda model and parameters from [9,16,23,33], $\mathcal{C}$ is approximately equal to 15. The parameter $\epsilon$ has the same values as for the Carreau model. Parameter $\zeta$ has the following values: $3.4552\times {10}^4$ for the carotid artery; $4.1662\times {10}^3$ for aorta and $1.6325\times {10}^4$ for iliac artery.

For the Cross model and parameters from [9], the value of $\mathcal{C}$ is estimated as 15. Parameter $\epsilon$ has the same values as for Carreau and Carreau–Yasuda model. For $\zeta$, the following values take place: $2.8088\times {10}^3$ for the carotid artery, $4.9312\times {10}^2$ for aorta and 0.0084, $1.5161\times {10}^3$ for the iliac artery. For the parameters from [16,17,20,35] $\mathcal{C}\approx 9.55$. Parameter $\epsilon$ has the same values. Parameter $\zeta$ is equal to $1.7801\times {10}^4$ for the carotid artery, $1.5301\times {10}^3$ for aorta and $7.4602\times {10}^3$ for the iliac artery.

For the Simplified Cross model, $\mathcal{C}$ is estimated as: $\mathcal{C}\approx 25$. Parameters have the following values: $\epsilon \approx 0.01968$, $\zeta \approx 1.7974\times {10}^4$ for the carotid artery; 0.0028, $3.3088\times {10}^3$ for aorta and 0.0029, $9.8666\times {10}^3$ for iliac artery.

For the case of the Modified Cross model $\mathcal{C}$ is estimated as 25 and $\epsilon$ has the same values as for the Cross model. For the carotid artery, $\zeta \approx 2.76116896\times {10}^9$; $4.70660529\times {10}^7$ for aorta and $6.52156883\times {10}^8$ for iliac artery.

For the Yeleswarapu model and parameters from [16,38], $\mathcal{C}\approx 13.7$ and $\epsilon \approx 0.01968$ as for the Simplified Cross model, $\zeta \approx 3.3275\times {10}^4$ for the carotid artery, 0.0028, $6.1254\times {10}^3$ for aorta and 0.0029, $1.8265\times {10}^4$ for the iliac artery.

For the Modified Yeleswarapu model, the parameters are dependent on the hematocrit H [21], so the dimensionless parameters $\epsilon$, $\xi$ and $\zeta$ are also dependent on its values. Let us obtain its values for two values of H: 0.4 and 0.5. At $H=0.4$, the value of $\mathcal{C}$ is approximately equal to 9.3876 and $\epsilon \approx 0.0175$, $\zeta \approx 1.9299\times {10}^4$ for the carotid artery, 0.0025, $3.5527\times {10}^3$ for aorta and 0.0026, $1.0594\times {10}^4$ for iliac artery. At $H=0.5$, the value of $\mathcal{C}$ is estimated as 13.3797, $\epsilon \approx 0.0228$, $\zeta \approx 2.0163\times {10}^4$ for the carotid artery, 0.0032, $3.7117\times {10}^3$ for aorta and $0.0033$, $\zeta \approx 1.1068\times {10}^4$ for iliac artery.

For the Quemada model, two parameters—$\epsilon$ and $\zeta$—are considered. For the parameters from [16,20,49], the following values are obtained: $\epsilon \approx 0.0047$, $\zeta \approx 34.57$ for the carotid artery, $6.81\times {10}^{-4}$, $14.83$ for aorta and $6.98\times {10}^{-4}$, 25.61 for the iliac artery.

As can be seen, the small parameter $\epsilon$ 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 $\zeta$ takes place. The first term $-\epsilon Q/A$ in the expression for $f(A,Q)$ prevails, and if the values of $\epsilon$ are close, it can lead to the close values of the solutions for different models.

3. 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\in [0,+\infty )$. 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\to +\infty$. 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.

3.1. Perturbation Method

Let the semi-infinite interval $z\in [0,+\infty )$ be considered. The initial conditions are presented as:

$$A(0,z)=A^0+\epsilon {\phi }_1\left(z\right)+{\epsilon }^2{\phi }_2\left(z\right)+\dots ,$$

(8)

$$Q(0,z)=Q^0+\epsilon {\psi }_1\left(z\right)+{\epsilon }^2{\psi }_2\left(z\right)+\dots ,$$

(9)

where ${\phi }_i\left(z\right)$,${\psi }_i\left(z\right)$ are bounded at $z\to \pm \infty$, $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:

$$Q(t,0)=\omega \left(t\right)=Q^0+\epsilon {\omega }_1\left(t\right)+{\epsilon }^2{\omega }_2\left(t\right)+\dots ,$$

(10)

where ${\omega }_i\left(t\right)$ 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:

$$A(t,z)=A_0(t,z)+\epsilon A_1(t,z)+{\epsilon }^2{A}_2(t,z)+\dots ,$$

(11)

$$Q(t,z)=Q_0(t,z)+\epsilon Q_1(t,z)+{\epsilon }^2{Q}_2(t,z)+\dots .$$

(12)

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 $\epsilon$ are considered. It is easy to obtain that $A_0(t,z)=A^0$, $Q_0(t,z)=Q^0$.

For the first-order terms $A_1$ and $Q_1$, the following problem takes place:

$$\frac{\partial A_1}{\partial t}+\frac{\partial Q_1}{\partial z}=0,\frac{\partial Q_1}{\partial t}+b_1\frac{\partial A_1}{\partial z}+b_2\frac{\partial Q_1}{\partial z}=f_0,$$

(13)

$$A_1(0,z)={\phi }_1\left(z\right),Q_1(0,z)={\psi }_1\left(z\right),$$

(14)

$$Q_1(t,0)={\omega }_1\left(t\right),$$

(15)

where $b_1$, $b_2$ and $f_0$ are defined as:

$$b_1=\chi \sqrt{A_0}-\alpha \frac{Q_0^2}{A_0^2},b_2=\frac{2\alpha Q_0}{A_0},f_0=\frac{f(A_0,Q_0)}{\epsilon }.$$

For the second-order terms $A_2$, $Q_2$, the following problem is considered:

$$\frac{\partial A_2}{\partial t}+\frac{\partial Q_2}{\partial z}=0,\frac{\partial Q_2}{\partial t}+b_1\frac{\partial A_2}{\partial z}+b_2\frac{\partial Q_2}{\partial z}=F(t,z),$$

(16)

$$A_2(0,z)={\phi }_2\left(z\right),Q_2(0,z)={\psi }_2\left(z\right),$$

(17)

$$Q_2(t,0)={\omega }_2\left(t\right),$$

(18)

where

$$F(t,z)={\beta }_1\frac{\partial \left(A_1{Q}_1\right)}{\partial z}+{\beta }_2\frac{\partial Q_1^2}{\partial z}+{\beta }_3\frac{\partial A_1^2}{\partial z}+{\delta }_1{A}_1+{\delta }_2{Q}_1,$$

where

$${\beta }_1=\frac{2\alpha Q_0}{A_0^2},{\beta }_2=-\frac{\alpha }{A_0},{\beta }_3=-\left(\alpha \frac{Q_0^2}{A_0^3}+\frac{\chi }{4\sqrt{A_0}}\right),$$

and constants ${\delta }_1$, ${\delta }_2$ for every model have its own expressions. For the Power Law model (including Newtonian at $n=1$), they are written as:

$${\delta }_1=-s^{\prime }\left(1\right)|s^{\prime }{\left(1\right)|}^{n-1}\frac{(3n-1)Q_0^n}{2A_0^{\frac{3n+1}{2}}},{\delta }_2=s^{\prime }\left(1\right){\left|s^{\prime }\left(1\right)\right|}^{n-1}\frac{n{Q}_0^{n-1}}{A_0^{\frac{3n-1}{2}}}.$$

For the Carreau–Yasuda model (including Carreau at $a=2$), the following expressions take place:

$${\delta }_1=-\frac{s^{\prime }\left(1\right)}{2A_0^2\left(A_0^{\frac{3a}{2}}+{\left|s^{\prime }\left(1\right)\right|}^a\zeta Q_0^a\right)}(2Q_0{A}_0^{\frac{3a}{2}}+2{\left|s^{\prime }\left(1\right)\right|}^a\zeta Q_0^{a+1}+$$

$$+\mathcal{C}{\left(1+|s^{\prime }{\left(1\right)|}^a\zeta Q_0^a{A}_0^{-\frac{3a}{2}}\right)}^{\frac{n-1}{a}}\left(2Q_0{A}_0^{\frac{3a}{2}}-|s^{\prime }{\left(1\right)|}^a\zeta Q_0^{a+1}+3{\left|s^{\prime }\left(1\right)\right|}^a\zeta Q_0^{a+1}n\right)),$$

$${\delta }_2=s^{\prime }\left(1\right)\frac{A_0^{\frac{3a}{2}}+{\left|s^{\prime }\left(1\right)\right|}^a\zeta Q_0^a+\mathcal{C}{\left(1+|s^{\prime }{\left(1\right)|}^a\zeta Q_0^a{A}_0^{-\frac{3a}{2}}\right)}^{\frac{n-1}{a}}\left(A_0^{\frac{3a}{2}}+{\left|s^{\prime }\left(1\right)\right|}^a\zeta n{Q}_0^a\right)}{A_0\left(A_0^{\frac{3a}{2}}+{\left|s^{\prime }\left(1\right)\right|}^a\zeta Q_0^a\right)},$$

for the Modified Cross model ($a=1$ corresponds to the Cross model, $a=1$, $m=1$—to the Simplified Cross model):

$${\delta }_1=-\frac{s^{\prime }\left(1\right)}{2A_0^2\left(A_0^{\frac{3m}{2}}+{\left|s^{\prime }\left(1\right)\right|}^m\zeta Q_0^m\right)}(2Q_0{A}_0^{\frac{3m}{2}}+2Q_0^{1+m}{\left|s^{\prime }\left(1\right)\right|}^m\zeta +$$

$$+\mathcal{C}\left(-3Q_0^{1+m}a|s^{\prime }{\left(1\right)|}^m\zeta m+2Q_0{A}_0^{\frac{3m}{2}}+2Q_0^{1+m}{\left|s^{\prime }\left(1\right)\right|}^m\zeta \right){\left(1+|s^{\prime }{\left(1\right)|}^m\zeta Q_0^m{A}_0^{-\frac{3m}{2}}\right)}^{-a}),$$

$${\delta }_2=s^{\prime }\left(1\right)\frac{1}{A_0\left(A_0^{\frac{3m}{2}}+{\left|s^{\prime }\left(1\right)\right|}^m\zeta Q_0^m\right)}(A_0^{\frac{3m}{2}}+{\left|s^{\prime }\left(1\right)\right|}^m\zeta Q_0^m+$$

$$+\mathcal{C}\left(-a|s^{\prime }{\left(1\right)|}^m\zeta Q_0^{m}m+A_0^{\frac{3m}{2}}+{\left|s^{\prime }\left(1\right)\right|}^m\zeta Q_0^m\right){\left(1+|s^{\prime }{\left(1\right)|}^m\zeta Q_0^m{A}_0^{-\frac{3m}{2}}\right)}^{-a}).$$

For the Yeleswarapu model and its modified form, ${\delta }_1$ and ${\delta }_2$ are presented as:

$${\delta }_1=-\frac{s^{\prime }\left(1\right)}{2A_0^{\frac{5}{2}}{\left(A_0^{\frac{3}{2}}+\left|s^{\prime }\left(1\right)\right|\zeta Q_0\right)}^2}(2Q_0\sqrt{A_0}(A_0^{\frac{3}{2}}+|s^{\prime }\left(1\right){\left|\zeta Q_0\right)}^2+$$

$$+\mathcal{C}{A}_0^2{Q}_0\left(2A_0^{\frac{3}{2}}+2\left|s^{\prime }\left(1\right)\right|\zeta Q_0+\left(2A_0^{\frac{3}{2}}-\left|s^{\prime }\left(1\right)\right|\zeta Q_0\right)ln\left(\frac{A_0^{\frac{3}{2}}+\left|s^{\prime }\left(1\right)\right|\zeta Q_0}{A_0^{\frac{3}{2}}}\right)\right)),$$

$${\delta }_2=s^{\prime }\left(1\right)\frac{A_0{\left(A_0^{\frac{3}{2}}+\left|s^{\prime }\left(1\right)\right|\zeta Q_0\right)}^2+\mathcal{C}{A}_0^{\frac{3}{2}}\left(A_0^{\frac{3}{2}}+\left|s^{\prime }\left(1\right)\right|\zeta Q_0+A_0^{\frac{3}{2}}ln\left(\frac{A_0^{\frac{3}{2}}+\left|s^{\prime }\left(1\right)\right|\zeta Q_0}{A_0^{\frac{3}{2}}}\right)\right)}{{\left(A_0^{\frac{3}{2}}+\left|s^{\prime }\left(1\right)\right|\zeta Q_0\right)}^2},$$

For the Quemada model, the following expressions are obtained:

$${\delta }_1=-\frac{s^{\prime }\left(1\right)}{A_0^{\frac{9}{4}}{\left(-2A_0^{\frac{3}{4}}-2\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q_0}+H{k}_0{A}_0^{\frac{3}{4}}+H{k}_{\infty }\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q_0}\right)}^3}\times$$

$$\left(2\left(A_0^{\frac{3}{4}}+\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q_0}\right)Q_0\right(-4A_0^{\frac{7}{4}}-8A_0\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q_0}+2A_0^{\frac{7}{4}}H{k}_0-$$

$$-A_{0}H{k}_{\infty }\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q0}-4A_0^{\frac{1}{4}}|s^{\prime }\left(1\right)|{\zeta }^2{Q}_0+5\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q_0}H{k}_0{A}_0+2A_0^{\frac{1}{4}}\left|s^{\prime }\left(1\right)\right|{\zeta }^2{Q}_{0}H{k}_{\infty }),$$

$${\delta }_2=\frac{s^{\prime }\left(1\right)}{A_0\sqrt{Q_0}{\left(-2A_0^{\frac{3}{4}}-2\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q_0}+H{k}_0{A}_0^{\frac{3}{4}}+H{k}_{\infty }\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q_0}\right)}^3}\times$$

$$\left(4\left(A_0^{\frac{3}{4}}+\sqrt{|s^{\prime }\left(1\right)|}\zeta \sqrt{Q_0}\right)\right(-2\sqrt{Q_0}{A}_0^{\frac{3}{2}}-4A_0^{\frac{3}{4}}\sqrt{|s^{\prime }\left(1\right)|}\zeta Q_0+\sqrt{Q_0}{A}_0^{\frac{3}{2}}H{k}_0-$$

$$-2|s^{\prime }\left(1\right)|{\zeta }^2{Q}_0^{\frac{3}{2}}+2\sqrt{|s^{\prime }\left(1\right)|}\zeta Q_{0}H{k}_0{A}_0^{\frac{3}{4}}+\left|s^{\prime }\left(1\right)\right|{\zeta }^2{Q}_0^{\frac{3}{2}}H{k}_{\infty }).$$

Let the initial functions in (14) and (17) be presented as:

$${\phi }_1\left(z\right)=\overline{A}{e}^{iz},{\psi }_1\left(z\right)=\overline{Q}{e}^{iz},{\phi }_2\left(z\right)=\overline{\overline{A}}{e}^{2iz},{\psi }_2\left(z\right)=\overline{\overline{Q}}{e}^{2iz},$$

where $i^2=-1$ and $\overline{A}$, $\overline{Q}$, $\overline{\overline{A}}$, $\overline{\overline{Q}}$ are the real constants.

For the simulation of the oscillatory behaviour of flow rate at $z=0$, it is proposed that function ${\omega }_1\left(t\right)$ in (15) is presented as:

$${\omega }_1\left(t\right)={\overline{q}}_0\left(t\right)+{\overline{q}}_1\left(t\right),$$

where ${\overline{q}}_0\left(t\right)$ and ${\overline{q}}_1\left(t\right)$ are the solutions of the following initial problems:

$${\dot{\overline{q}}}_0=f_0,{\overline{q}}_0\left(0\right)=0,$$

$${\dot{\overline{a}}}_1=-i{\overline{q}}_1,{\dot{\overline{q}}}_1=-i{b}_1{\overline{a}}_1-i{b}_2{\overline{q}}_1,{\overline{a}}_1\left(0\right)=\overline{A},{\overline{q}}_1\left(0\right)=\overline{Q}.$$

It is easy to obtain that:

$${\overline{q}}_0\left(t\right)=f_{0}t,{\overline{a}}_1\left(t\right)=C_1{e}^{{\gamma }_{1}t}+C_2{e}^{{\gamma }_{2}t},{\overline{q}}_1\left(t\right)=C_{1}i{\gamma }_1{e}^{{\gamma }_{1}t}+C_{2}i{\gamma }_2{e}^{{\gamma }_{2}t},$$

where ${\gamma }_{1,2}=\frac{i}{2}\left(-b_2\pm \sqrt{4b_1+b_2^2}\right)$ and

$$4b_1+b_2^2=4\left(\chi \sqrt{A_0}+({\alpha }^2-\alpha ){\left(\frac{Q_0}{A_0}\right)}^2\right)>0,$$

at $\alpha \ge 1$. Constants $C_1$ and $C_2$ are computed as:

$$C_1=-\frac{\overline{A}{\gamma }_2+i\overline{Q}}{{\gamma }_1-{\gamma }_2},C_2=\frac{\overline{A}{\gamma }_1+i\overline{Q}}{{\gamma }_1-{\gamma }_2}.$$

The solution of the problem (13)–(15) is written as:

$$A_1(t,z)={\overline{a}}_1\left(t\right)e^{iz},Q_1(t,z)={\overline{q}}_0\left(t\right)+{\overline{q}}_1\left(t\right)e^{iz}.$$

Let the function ${\omega }_2\left(t\right)$ be presented as:

$${\omega }_2\left(t\right)={\overline{\overline{q}}}_0\left(t\right)+{\overline{\overline{q}}}_1\left(t\right)+{\overline{\overline{q}}}_2\left(t\right),$$

where ${\overline{\overline{q}}}_0$, ${\overline{\overline{q}}}_1$, ${\overline{\overline{q}}}_2$ are the solutions of the following initial problems:

$${\dot{\overline{\overline{q}}}}_0={\delta }_2{f}_{0}t,{\overline{\overline{q}}}_0\left(0\right)=0,$$

$${\dot{\overline{\overline{a}}}}_1=-i{\overline{\overline{q}}}_1,{\dot{\overline{\overline{q}}}}_1=-i{b}_1{\overline{\overline{a}}}_1-i{b}_2{\overline{\overline{q}}}_1+g_1\left(t\right),{\overline{\overline{a}}}_1\left(0\right)={\overline{\overline{q}}}_1\left(0\right)=0,$$

(19)

where $g_1\left(t\right)={\delta }_1{\overline{a}}_1\left(t\right)+{\delta }_2{\overline{q}}_1\left(t\right)+i{\beta }_1{\overline{q}}_0\left(t\right){\overline{a}}_1\left(t\right)+2i{\beta }_2{\overline{q}}_0\left(t\right){\overline{q}}_1\left(t\right)$,

$${\dot{\overline{\overline{a}}}}_2=-2i{\overline{\overline{q}}}_2,{\dot{\overline{\overline{q}}}}_2=-2i{b}_1{\overline{\overline{a}}}_2-2i{b}_2{\overline{\overline{q}}}_2+g_2\left(t\right),{\overline{\overline{a}}}_2\left(0\right)=\overline{\overline{A}},{\overline{\overline{q}}}_2\left(0\right)=\overline{\overline{Q}},$$

(20)

where $g_2\left(t\right)=2i{\beta }_1{\overline{a}}_1\left(t\right){\overline{q}}_1\left(t\right)+2i{\beta }_2{\overline{q}}_1^2\left(t\right)+2i{\beta }_3{\overline{a}}_1^2\left(t\right)$.

The initial problems (19) and (20) can be rewritten as:

$$\dot{\mathbf{y}}=\mathbf{P}\mathbf{y}+\mathbf{g}\left(t\right),\mathbf{y}\left(0\right)={\mathbf{y}}_0,$$

(21)

where $\mathbf{P}$ is a constant complex matrix.

The solution of (21) is obtained as:

$$\mathbf{y}\left(t\right)=\mathbf{Y}\left(t\right)\left({\mathbf{Y}}^{-1}\left(0\right){\mathbf{y}}_0+\underset{0}{\overset{t}{\int }}{\mathbf{Y}}^{-1}\left(\tau \right)\mathbf{f}\left(\tau \right)d\tau \right),$$

where $\mathbf{Y}\left(t\right)$ is the fundamental matrix of system $\dot{\mathbf{y}}=\mathbf{P}\mathbf{y}$. The expression for ${\overline{\overline{q}}}_0\left(t\right)$ is written as: ${\overline{\overline{q}}}_0\left(t\right)=\frac{{\delta }_2{f}_0{t}^2}{2}$, expressions for ${\overline{\overline{a}}}_1\left(t\right)$, ${\overline{\overline{q}}}_1\left(t\right)$, ${\overline{\overline{a}}}_2\left(t\right)$, ${\overline{\overline{q}}}_2\left(t\right)$ are not presented according to its cumbersomeness, but they can be easily obtained using the systems of symbolic computations.

The solution of problem (16)–(18) is presented as:

$$A_2(t,z)={\overline{\overline{a}}}_1\left(t\right)e^{iz}+{\overline{\overline{a}}}_2\left(t\right)e^{2iz},Q_2(t,z)={\overline{\overline{q}}}_0\left(t\right)+{\overline{\overline{q}}}_1\left(t\right)e^{iz}+{\overline{\overline{q}}}_2\left(t\right)e^{2iz}.$$

3.2. 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:

$$\mathrm{NNF}=\frac{\left|\right|Re\left(Q(t,z)\right)-Re\left(Q_{Newt}(t,z)\right)\left|\right|}{\left|\right|Re\left(Q_{Newt}(t,z)\right)\left|\right|}\times 100\%,$$

where $L_2$ norm is used, $Q_{Newt}$ is the solution for the Newtonian model. The values of $Re\left(Q\right)$ 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 $\epsilon$ 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: $A^0=\pi$, $Q^0=1$, $\overline{A}=\overline{Q}=\overline{\overline{A}}=\overline{\overline{Q}}=1$, t is considered on the time interval $[0,2\pi ]$, $z\in [0,15]$. The plots of the real part of $\omega \left(t\right)$ and corresponding $Re\left(A\right(t,0\left)\right)$ 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\left(Q\right)$ 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\left(Q\right)$ at selected values of d and H are presented in Figure 7 in comparison with the solution for the Newtonian case.

4. 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.