Clinical Hemorheology and Microcirculation 62 (2016) 109–121 DOI 10.3233/CH-151955 IOS Press
109
Numerical simulation of blood flow through a capillary using a non-linear viscoelastic model Amin Shariatkhaha , Mahmood Norouzia,∗ and Mohammad Reza Heyrani Nobarib a b
Department of Mechanical Engineering, Shahrood University of Technology, Shahrood, Iran Department of Mechanical Engineering, Amirkabir University of Technology, Tehran, Iran
Abstract. In this article, a periodic developing blood flow in a capillary is simulated using a non-linear viscoelastic model for the first time. Here, the Giesekus model is used as the constitutive equation, and based on the experimental data, the best value for the mobility factor and zero shear rate viscosity are derived. The numerical solution of the problem is obtained using the finite volume method. The algorithm of the solution is pressure implicit with splitting of operators (PISO). The simulation carried out using the Giesekus, Oldroyd-B and Newtonian models and the results indicate that the Giesekus model presents a more accurate solution for the stress and velocity fields than the Newtonian and Oldroyd-B models. The previous studies on this problem were restricted to the linear and quasi-linear viscoelastic models. It is shown that only non-linear viscoelastic models can accurately describe the experimental data of unsteady blood flow in capillaries. Keywords: Giesekus model, oscillating flow, blood viscosity, mobility factor, capillary
1. Introduction Clinical science strongly depends on the experiments, which are usually very difficult to set up; therefore, the accurate simulations of human physiological processes can benefit the clinical researchers. However, the numerical results will be reliable if they provide a suitable approximation for the available experimental data. Blood is a combination of red blood cells (RBC), white blood cells (WBC) and plasma [12]. Although, 55% of blood is plasma, which is mostly water (92% by volume) [12] and appears to be a Newtonian fluid, the existence of RBC and WBC make blood very complicated. Arising from the complex characteristics of the blood, it cannot be simulated without considering the effects of those cells on the properties. One of the most important characteristics is related to the red blood cells aggregation which cause the non-linearity and non-Newtonian properties of blood [9, 20, 23, 30, 40, 41]. At large shear rates [1, 8], blood can be assumed as a Newtonian fluid [14, 16–18]. When heart pumps blood, the aorta delivers blood into the vessel network, where the size of vessels decrease but the number of them increases. The total cross sectional area of these small vessels is 1000 times greater than the aorta [12]. It is clear that since the blood volume is constant, by increasing the cross sectional area, the velocity becomes slower. Hence, the low velocity makes shear rate of flow small and the blood inside capillaries can be no longer investigated accurately via Newtonian assumption. ∗ Corresponding author: Mahmood Norouzi, Department of Mechanical Engineering, Shahrood University of Technology, Shahrood, Iran. Tel.: +98 9123726933; Fax: +98 2733335445; E-mails: [email protected] and [email protected].
1386-0291/16/$35.00 © 2016 – IOS Press and the authors. All rights reserved
110
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
Deformation and aggregation of RBCs and also the large diameter ratio of RBCs to capillaries are the reason of some non-Newtonian characteristics of blood flow [30]. Previous studies showed that blood has viscoelastic properties [30, 40]. RBCs in blood are similar to polymeric additives. They make blood a viscoelastic fluid in which temporal effects (relaxation and retardation time) and other mechanical properties (dependency of stress on the deformation in addition to the deformation rate) [5] should be considered. Although many studies considered blood as a generalized Newtonian fluid (GNF) and employed models such as Power-law, Herschel Bulky, Cross and Casson [4, 6, 31, 36–37], studying blood becomes crucial especially in complex flows like stenosis and branching of capillaries, where simple models are unable to describe the physics of flow accurately. Upper Convected Maxwell (UCM) and Oldroyd-B are two quasi-linear viscoelastic constitutive equations that have been used in former studies [32]. Viscoelastic properties of fluids are like a discrete spectrum and various elements should be used to cover the relaxation spectra [5]. Thurston [41] used a generalized Maxwell model which is combination of six Maxwell elements and obtained the constants of this model based on the oscillation test of blood. Hern´an et al. [35] used a generalized Maxwell model similar to the work of Thurston [41]. Owens [30], Gonzalez and Owens [21] and Gonzalez et al. [22] introduced new viscoelastic model based on the temporal network theory which was used previously for modeling the polymeric flows. According to the literature, the previous studies about blood flow are restricted to quasi-linear viscoelastic and generalized Newtonian models. In the present study, the Giesekus constitutive equation (a non-linear viscoelastic model) is employed to simulate the unsteady blood flow in an arteriole for the first time. The empirical investigations indicated that the blood flow shows the non-linear viscoelastic behavior in capillaries. It is important to remember that the non-linear viscoelastic properties of blood flow cause some phenomena: 1. The non-linear viscosity which was previously studied by simple models such as power law. This property causes the blunting of velocity profiles. 2. The first and second normal stress differences. This property could not be modeled via pure viscose and quasi-linear viscoelastic constitutive equations. 3. A phase difference between velocity and stress fields in periodic viscoelastic flows. This effect is happen due to the relaxation spectra of viscoelastic materials. In this paper, the non-linear Giesekus constitutive equation is used to model all of above mentioned effects for unsteady blood flow in a capillary. This model is also able to spot the effect of the third invariant of shear rate in stress tensor and elongational effect. The constants of Giesekus model are obtained by implementing an optimization process on experimental data of rheological tests. Here, the governing equations are solved using finite volume method based on the pressure implicit with splitting of operators (PISO) algorithm. The numerical simulations are also performed via Newtonian and Oldroyd-B model and the results are compared with Giesekus constitutive equation. 2. Methods 2.1. Viscosity of blood Viscosity of blood depends on different parameters such as hematocrit, shear rate and diameter of vessel [11, 20, 24, 34]. Hematocrit of normal human is 45% [12] and its high shear viscosity for vessels larger than 20 m is about 0.0032 pa.s [11, 34]. It is possible to consider the blood as a dilute viscoelastic solution in which the Plasma is the Newtonian solvent and other components play the role 0.001 of viscoelastic additive. Based on this assumption, the viscosity ratio of blood is ηηps = 0.0022 [21, 22,
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
111
Fig. 1. Diagram of blood viscosity versus the shear rate for three rheological models and experimental data [46].
30], where ηs is the plasma viscosity (as the solvent) and ηp is the viscosity of other components (as the viscoelastic additives). In this paper, the constants of Giesekus constitutive equation are calculated optimally based on the non-linear viscosity of blood (see Fig. 1). Here, these constants are obtained using the least square of error between the experimental data and viscometric response of Giesekus model in steady shear test [46]. The results of optimization process indicated that α ≈ 0.5 and η0 = 0.0063pa.s are the best approximations. These are close to the results of similar study done by Hakim et al. [13] to find the constants of Giesekus model. In order to obtain the whole diagram of viscometric function of a material, specific kinds of rheometer should be used at any range of shear rates [7]. However, the former researchers reported only the viscosity of blood using one kind of viscometer [15, 20, 46]. The experimental data which was used in this study is obtained using Couette viscometer [46]. This viscometer is not suitable for very small shear rate [7], so noticeable deviation is expected in this region. According to Fig. 1, it is clear that the Newtonian constitutive equation is not a suitable model for simulating the blood flow at small shear rates where the changes of viscosity are significant, but it has an acceptable approximation at large enough shear rates. The quasi-linear viscoelastic models (such as UCM, second-order and Oldroyd-B models) predict the elastic material modules of blood flow but they considered a constant viscosity for all shear rates in steady shear flows [5]. Therefore, it would be expected that the non-linear constitutive equations simulate the blood flow in capillaries better than the Newtonian or quasi-linear models.
2.2. Governing equations The governing equations consist of the continuity and momentum equations are presented as: ∇ · v = 0
ρ
∂v + v · ∇v = ∇ · T ∂t
(1a) (1b)
where v = (v r , v θ , v z ) is the velocity vector, T , t and ρ are the Cauchy stress tensor, time and density of fluid, respectively. Assuming the body force to be zero and the flow to be incompressible, T is the sum of hydrostatic pressure and extra stress:
112
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
T = −pI + τ
(2)
where τ, p and I are extra stress, hydrostatic pressure and identity matrix. Cylindrical coordinate system (r, θ, z) is used to analyze the flow in a straight circular pipe, where the flow is axisymmetric (∂ ∂θ = 0 , v θ = 0) [47]. Using the axisymmetric condition, the governing equations are simplified as follows: 1 ∂ ∂v z (rv r ) + =0 r ∂r ∂z ∂v r 1 ∂ ∂v r ∂v r ∂ = + vr + vz (rTrr ) + Tzr ρ ∂t r ∂r ∂r ∂z ∂z
∂v z ∂v z ∂v z ρ + vr + vz ∂t ∂r ∂z
=
1 ∂(rTrz ) ∂ + Tzz r ∂r ∂z
(3a) (3b)
(3c)
2.3. Boundary conditions In this study, the inlet velocity is uniform and changes with time. Unlike the previous works, instead of using a cosine function for temporal changes of inlet velocity, the real data of Tsukada [42] is used. There are many experimental data which report velocity profile at some specific time [8, 10, 36] or velocity at some specific point for an interval time [42], but there is no information regarding the changes of mean velocity with time, which is needed in this study for inlet velocity. Therefore, an approximate relation was derived between mean and maximum of velocity profiles [37]. After that in order to have mean velocity versus time, this relation was applied to temporal data of Tsukada [42] which reports changes of maximum velocity in arteriole. At the outlet, the flow is assumed to be fully developed. Also zero gradient condition is used for the velocity and stress components in flow direction. Furthermore, no-slip condition is used for the velocity at walls. On the axis of symmetry, vr and radial gradient of vz are zero. Pressure has zero gradients across all boundaries. The stress field is known at the inlet to satisfy the prescribe inlet velocity. The zero value for τrz and τrθ , and zero radial gradient for other stress components are applied at the symmetry boundary. Also, the zero gradients are used for stress field at walls. 2.4. Initial condition The Giesekus constitutive equation has high order non-linear terms and the convergence of numerical simulation is highly depended on the initial condition. To overcome this problem, the initial condition for the velocity is set as the inlet velocity, while for pressure and stress are set as zero [37, 45]. Due to the periodic behavior of flow field, the effect of initial condition on the final solution is completely waned after some cycles. 2.5. Constitutive equation There are some families of quasi-linear and non-linear viscoelastic models which are obtained based on the different approaches in rheology [2]. The Oldroyd-B constitutive equation is a well-known quasi-linear model. The model is derived by replacing the time derivative of stress by upper convected derivative in Jeffreys model. The Oldroyd-B model is suitable for dilute viscoelastic solutions and its extra stress is usually derived by considering the Newtonian behavior for solvent and UCM model for viscoelastic additives:
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
113
τ = τs + τp
(4a)
τs = ηs γ˙
(4b)
τp + λ1 τp(1) = η0 γ˙
(4c)
where τ is the extra stress tensor, λ1 is the relaxation time, η0 is the zero shear rate viscosity, γ˙ is the shear rate tensor, ‘s’ is the solvent contribution (Newtonian) and ‘p’ is the polymeric contribution (viscoelastic). Here, τ(1) is the upper convected time derivative of stress [5] and can be expressed as τ(1) =
D τ − {(∇v)† . τ + τ. (∇v)} Dt
(5)
D is the substantive derivative, ∇v is the tensor of derivative for the fluid and superscript “†” where Dt means transpose of the tensor. In a steady shear flow test, this model presents the following viscometric functions:
η = η0
(6a)
1 = 2η0 (λ1 − λ2 )
(6b)
2 = 0
(6c)
where 1 and 2 are the coefficients of the first and second normal stress differences and λ2 is the retardation time. The model is linear because the viscosity and first normal stress difference coefficient are revealed to be constant. Also the second normal stress coefficient is constant and zero. In order to have a non-linear model, a non-linear term of stress should be added to the constitutive equation. This modification is implemented in Giesekus model [5]: τp + λ1 τp(1) + α
λ1 {τp . τp } = ηp γ˙ ηp
(7)
The Giesekus model contains four parameters: a relaxation time (λ1 ), the solvent and polymeric contribution of zero shear rate viscosity (ηs and ηp ) and mobility factor (α). The origin of the term involving α can be associated to the anisotropic Brownian motion and/or anisotropic hydrodynamic drag on the constituent polymer molecules [5]. In a steady shear flow test, Giesekus model for a simplified case (λ2 = 0, α = / 0, 1 and γ˙ → ∞) present the non-linear viscometric functions: η∼
(1 − α)/α
√
η0 λ1 γ˙
2 η0 λ1 (α(1 − α))1/4 ˙ 3/2 α (λ1 γ) η0 λ1
2 ∼ − ˙ 2 (λ1 γ)
1 ∼
(8a) (8b) (8c)
Unlike the Oldroyd-B model, the viscometric functions of Giesekus model shows a fractional dependency on shear rate which can be fitted appropriately on experimental data. It should be noted that this formulation is defined for polymeric solutions and extra stress is the sum of polymeric and Newtonian contributions. The polymeric part of solution causes the viscoelastic properties. This is the same as blood in which the RBCs plays the viscoelastic role while Plasma (and other components with small molecular weight) acts as the Newtonian counterpart. The Giesekus constitutive equation can be used in various modes to model the relaxation spectra. Normally, viscoelastic
114
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
fluids have a spectrum of relaxation time instead of only one mode [5]. Some dominant modes can be used if more accurate simulation is required. 2.6. Numerical method In this study, finite volume method is used to solve the governing equations. Algorithm of solution is chosen as Pressure Implicit with Splitting of Operators (PISO). PISO algorithm is a pressure-velocity calculation algorithm, which has one predictor step and two corrector steps and may be assumed as an extended version of Semi-Implicit Method for Pressure Linked Equations (SIMPLE) [43]. It’s possible to use more than two corrector steps for a rapid convergence. All of flow parameters are discretized using central differences, except for the convection terms which are approximated by the linear-upwind differencing scheme (LUDS). This is a generalization of the well-known upwind differencing scheme (UDS), where the value of a convected variable at a cell face location is given by its value at the first upstream cell center. In the linear-upwind differencing scheme, the value of that convected variable at the same cell face is given by a linear extrapolation based on the values of the variable at the two upstream cells. In general, LUDS is the second-order accurate, as compared with first-order accuracy of UDS, and thus, the problem of numerical diffusion would be reduced by using it [29]. The time discretization is carried out using the second order Crank Nicholson scheme. It should be noted that all of the governing equations can be written in the form of the general transport equation as follows: ∂ ∂t
∂ Aφ) + (AUj φ ∂xj
∂ = ∂xj
∂φ ∂xj
+ Sφ
(9)
where the working variable φ is a component of a vector or tensor and even an unknown constant, U is velocity and t is time. The coefficients A and have different meanings for different working variables and Sφ is a source term which includes all the terms that cannot be accommodated in the convective and diffusion terms. This has different contents for different equations. The flow domain is divided into a set of non–overlapping control volumes ∇V with bounding surface area [44]. By Integrating Equation (9) over any control volume and a time increment δt and using the divergence theory, we have: A
V (φp − φpn ) + δt
(Auj φ −
∂φ )nj dA = S¯ φ ∂xj
(10)
where the subscript p refers to the grid point, superscript n denotes the value evaluated at time–level n and S¯ φ is the volume integral of source term Sφ . This can be linearized as: S¯ φ =
Sφ dV = S¯ c + S¯ p φp
(11)
where S¯ c does not explicitly depend on φ and S¯ p is the coefficient of φp , which is made negative to enhance the numerical stability of the discretized equation system. Using a proper spatial variation approximation scheme, the final discretized equation, which relates φp to its neighboring grid point values, can be symbolically expressed in a general form for each control volume: ap φp =
nb
anb φnb + S¯ c + ap0 apn
(12)
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
115
where ‘nb’ means neighborhood and ap0 = A
V , ap = anb + ap0 − S¯ p δt nb
(13)
It is then solved by PISO iterative algorithm. The system of linear equations for pressure is solved by Multi-Grid (MG) algorithm. The convergence criterion of numerical solution is reaching to the residual less than 5 × 10−8 . The system of linear equations for the velocity and stress is solved using Cholesky algorithm with the residual of 5 × 10−7 . The maximum time interval (δt) is set as 5 × 10−8 s for the Giesekus model, 10−6 s for the Oldroyd-B model, and 10−5 s for the Newtonian model. The shorter value of δt for Giesekus model is related to the higher order non-linearity of this constitutive equation. The vessel is simulated as a straight pipe with circular cross section. Due to the axisymmetry condition, a 2D rectangular mesh was used. The ratio of duct length to its radius is considered equal to 40 and its radius is 25 m. The reference of cylindrical coordinate system is located at the left bottom corner of rectangular. The length of pipe is considered long enough to satisfy the fully developed conditions. A structured uniform mesh with rectangular cells is used, where 15 nodes in the r direction and 150 nodes in the z direction are considered. The density of blood is 1053 kg.m–3 [16, 30], its viscosity at large shear rates (η∞ ) is set to 0.0032 pa.s [11, 12, 30, 34], and the relaxation time of blood is assumed to be 0.008 s [20]. 3. Results and discussion 3.1. Grid study and validation In order to check the grid dependency of numerical solution, we used four types of meshes with different grid numbers. Characteristics of these grids are: M1 (5∗ 10 nodes), M2 (10∗ 100 nodes), M3 (15∗ 150 nodes) and M4 (20∗ 200 nodes). In order to have independent results from mesh size, these grids are tested for a known flow and finally one of them is selected to be used for other numerical simulations. The selected grid has an acceptable accuracy while it does not have excessive runtime. The flow of Oldroyd-B model with constant inlet velocity is chosen for the test case. Here, the maximum value of main flow velocity in steady fully developed condition is used to check the accuracy of different grids. The exact analytical solution of this problem indicates that the ratio of maximum value of axial velocity to inlet velocity is 2.0 for steady fully developed Oldroyd-B fluid flow in a straight pipe [5]. The solution of Oldroyd-B flow is obtained by considering the zero value for α parameter of Giesekus constitutive equation in our CFD code. The error of maximum fully developed axial velocity according to the analytical solution and runtime of simulation, in order are: M1 (26.65%, 69 min), M2 (3.9%, 175 min), M3 (0.75%, 238 min) and M4 (0.35%, 491 min). It’s clear that M1 and M2 have considerable errors while the accuracy of M3 and M4 is acceptable (less than 1%). Because the runtime of M3 is less than the half of M4 so M3 is selected as the reference mesh grid structure. It is also possible to check the order of CFD simulation for non-Newtonian flows based on the grid numbers [25–28]. Applying the Richardson extrapolation on data of these four types of mesh, the order of numerical solution is obtained around 2. 3.2. Numerical simulation In this section, the effect of non-linearity of viscoelastic model on periodic blood flow is studied and its results are compared by the results of Newtonian and quasi-linear models which are used in previous studies. The numerical simulation indicated that the Newtonian and Oldroyd-B cases run about 18.4 and
116
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
Fig. 2. Velocity profile in two specific times: (a) Diastole (minimum velocity) and (b) Systole (maximum velocity).
6.2 times faster than the Giesekus case, respectively. In other words, the CFD simulation of unsteady blood flow by Giesekus model is slower and more unstable which is related to the non-linearity of this model. Figure 2 shows the systolic and diastolic velocity profiles (maximum and minimum value of velocity in a period) for different rheological models. According to the figure, the velocity profile of Giesekus model is flatter than the Newtonian and Oldroyd-B fluids. This type of velocity distribution is attributed to the shear thinning behavior of Giesekus model that provides larger values of viscosity in the vicinity of central region of vessel and consequently the small value of shear rate. The Maximum to mean ratio of velocity (a parameter that shows flatness of velocity profile) is approximately 1.5 for blood in experimental studies [10, 36]. This ratio is exactly 2 for Newtonian fluid flow [5]. According to the Fig. 2, this parameter is exactly 2 for Newtonian model (maximum velocity is twice of mean (inlet) velocity) and for the Giesekus model is closer to experimental value so the velocity profile of Giesekus model is shown to fit better with the experimental observations [10, 36] of blood flow in vessels in comparison to other models. According to the study of Popel and Johnson [33], the velocity profile becomes flatter in vessels smaller than 32 m, where the size of RBCs becomes more significant fraction of vessel diameter. Mean difference between results of Giesekus and Oldroyd-B model is 2.1% in Fig. 2a and 1.09% in Fig. 2b. Mean difference between results of Giesekus and Newtonian model is 4.12% in Fig. 2a and 2.24% in Fig. 2b. Figure 3a and 3b shows systolic and diastolic velocity contours for Giesekus model. Here it’s shown again that velocity become fully developed so quickly. According to these figures, it seems that selection of 40 as the value of geometry aspect ratio looks overestimated. Developing region of shear stress of Giesekus model is shown in Fig. 3c and 3d. It’s clear that the developing length is longer for stress and suitable length is selected for the pipe. The stress has a maximum in the upper corner of the entrance. In this region, the inlet velocity faces the wall and has to become zero; therefore high shear rate exists in this region. Figure 4 compares developing length for velocity and stress in centerline for systolic flow. Here, it is shown that for having fully developed flow in outlet, the shear stress is the determinant
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
117
Fig. 3. Contours of velocity for blood flow using Giesekus model: (a) Diastole and (b) Systole. Contours of shear stress for blood flow using Giesekus model: (c) Diastole and (d) Systole.
Fig. 4. Comparing developing length of velocity with shear stress of Giesekus model in systole (maximum velocity).
parameter and not the velocity. Therefore, the blood flow in short vessels has not enough length to get to the fully developed situation and this point should be taken into account in CFD simulations. The history of the Giesekus, Oldroyd-B and Newtonian velocity at the center of pipe in fully developed conditions (at downstream) and also inlet velocity (at upstream) are shown in Fig. 5. The fully
118
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
Fig. 5. History of inlet velocity and velocity in center of vessel for three models.
Fig. 6. History of wall shear stress in center of vessel for three models.
Fig. 7. Profile of Newtonian part, viscoelastic part and total shear stress of Giesekus model in systole.
developed velocities are shown to be larger than the inlet velocity because the growth of boundary layers is completed in this region. According to this figure, the maximum value of Newtonian velocity is exactly twice the size of inlet velocity and is larger than the velocity derived by Giesekus or Oldroyd-B constitutive equations. Algorithm of changes and time of minimum and maximum values are the same for three models. Mean difference between results of Giesekus and Oldroyd-B model is 2.96% and between Giesekus and Newtonian model is 6.93%. One of the most important parameters in studying the cardiovascular systems is the stress of walls of vessel, which is significant in many diseases [8]. Shear stresses on the walls are shown in Fig. 6. The velocity gradients near the walls for the three models are very close, but it’s seen that the wall shear stress is smaller for the Giesekus case. Mean difference between results of Giesekus and Oldroyd-B model is 8.89% and between Giesekus and Newtonian model is 19.32% and seems significant.
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
119
Fig. 8. Diagrams of shear rate and velocity profile versus the radius of vessel for systole.
This event interestingly match with the characteristics of blood flow in small vessels: the presence of RBCs in center of vessel and existence of a layer of plasma near to the wall [3]. In this scenario, blood is divided into two independent parts: RBCs (as viscoelastic additives) and plasma (as a Newtonian solvent). The concentration of RBCs is very low in vicinity of walls so the Newtonian part of stress is the dominant term in this region. However, the Newtonian part of Giesekus wall shear stress (dominant part) is smaller than the pure Newtonian case (wall shear rate is almost equal but viscosity of Newtonian part of Giesekus case is smaller than pure Newtonian case), and the viscoelastic part of stress has a more concentrated distribution in the vicinity of center line of vessel due to presence of RBCs (see Fig. 7). In other words, the viscoelastic portion of Giesekus shear stress obtains a small value around the wall because the shear rate obtains its maximum value and the viscometric function of this constitutive equation is a fractional shear-thinning model. Therefore, it is expected that the total wall shear stress to be smaller in the Giesekus fluids. Figure 8 illustrates the diagrams of shear rate in a cross section versus the radius of vessel. This graph is presented for systole in fully developed region. As can be seen, the shear rate is in order of 103 (s–1 ) at the region near the wall so the Giesekus model can be more reliable while the linear and quasi-linear models can not be suitable choices in this range of shear rate. 4. Conclusions In this paper, the suitability of linear and non-linear constitutive equations in simulation of the unsteady blood flow in capillaries is investigated. It is shown that the Newtonian model is more suitable in the large shear rate cases. On the other hand, at small shear rates, the Giesekus model can more accurately estimate the unsteady blood flow behavior. The Giesekus model gives a flatter profile of velocity, which has revealed a better agreement with experimental observations [10, 33, 36]. The CFD simulations indicate that the results of Oldroyd-B model stay between the Newtonian and Giesekus models. The wall shear stress in the Giesekus model coincides with the other studies concerning with the algorithm of RBCs movement in small vessels. As it is shown the Newtonian model cannot describe stress in small vessels correctly so it recommends utilizing the complicated viscoelastic models. References [1] M. Anand and K.R. Rajagopal, A shear-thinning viscoelastic fluid model for describing the flow of blood, Int J Cardiovasc Med Sci 4 (2004), 59–68. [2] H.A. Barnes, J.F. Hutton and K. Walters, Rheology series - volume 3: An introduction to rheology. 1st ed. Elsevier, Amsterdam (1989), 152–156.
120
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
[3] O. K. Baskurt, H. J. Meiselman, Blood rheology and hemodynamics, In Seminars in thrombosis and hemostasis (Vol 29). New York, stratton intercontinental medical book corporation, c1974-(2003), 435–450. [4] O.K. Baskurt, M. Hardeman, M.W. Rampling and H.J. Meiselman, Handbook of hemorheology and hemodynamicsvolume 69 biomedical and health research. IOS Press, Amsterdam (2007). [5] R.B. Bird, R.C. Armstrong and O. Hassager, Dynamics of polymeric liquids- Volume 1: fluid mechanics. Wiley (1987). [6] J. Chen and L.U. Xi-Yun, Numerical investigation of the non-Newtonian pulsatile blood flow in a bifurcation model with a non-planar branch, J Biomech 39 (2006), 818–832. [7] R.P. Chhabra, J.F. Richardson, Non-newtonian flow and applied rheology – Engineering application. 2nd ed. Elsevier (2008), 5–8. [8] S. Chien, R.G. King, R. Skalak, S. Usami and A.L. Copley, Viscoelastic properties of human blood and red cell suspensions, Biorheology 12 (1975), 341–346. [9] S. Cho, B. Namgung, H.S. Kim, H.L. Leo and S. Kim, Effect of erythrocyte aggregation at pathological levels on NO/O2 transport in small arterioles. Clin Hemorheol Microcirc 59 (2015), 163–175. [10] S. Einav, H.J. Berman, R.L. Fuhro, P.R. DiGiovanni, S. Fine, J. D. Fridman, Measurement of velocity profiles of red blood cells in the microcirculation by laser doppler anemometry (LDA). Biorheology 12 (1975), 207-210. [11] R. F˚ahræus and T. Lindqvist, The viscosity of the blood in narrow capillary tubes, Am J Physiol 96 (1931), 562–568. [12] A.C. Guyton and J.E. Hall, Textbook of Medical Physiology. 11th ed. Elsevier, Saunders (2006), 162–169. [13] S. Hakim, J. Morshedian, M. R. G. Narenji, P. K. Nia, Rheological modelling of Caspian pony blood. Iran Polym J 10(5) (2001), 293–303. [14] Z. Keshavarz-Motamed and L., Kadem, 3D pulsatile flow in a curved tube with coexisting model of aortic stenosis and coarctation of the aorta. Med Eng Phys 33(3) (2011), 315–324. [15] S. Kim, A study of non-Newtonian viscosity and yield stress of blood in a scanning capillary-tube rheometer. PhD thesis, Drexel University (2002). [16] R.E. Klabunde, Cardiovascular physiology concepts. Lippincott, Williams & Wilkins (2011). [17] D.N. Ku, Blood flow in arteries, Annu Rev Fluid Mech 29 (1997), 399–434. [18] M.X. Li, J.J. Beech-Brandt, L.R. John, P.R. Hoskins and W.J. Easson, Numerical analysis of pulsatile blood flow and vessel wall mechanics in different degrees of stenosis, J Biomech 40 (2007), 3715–3724. [19] D.S. Long, M.L. Smith, A.R. Pries, K. Ley and E.R. Damiano, Microviscometry reveals reduced blood viscosity and altered shear rate and shear stress profiles in microvessels after hemodilution, Proc Natl Acad Sci 101 (2004), 10060–10065. ¨ [20] J.A. Long, A. Undar, K.B. Manning and S. Deutsch, Viscoelasticity of pediatric blood and its implications for the testing of a pulsatile pediatric blood pump, ASAIO J 51 (2005), 263–566. [21] M.A. Moyers-Gonzalez and R.G. Owens, A non-homogeneous constitutive model for human blood Part II. Asymptotic solution for large P´eclet numbers, J Non-Newtonian Fluid Mech 155 (2008), 146–160. [22] M.A. Moyers-Gonzalez, R.G. Owens and J. Fang, A non-homogeneous constitutive model for human blood Part III. Oscillatory flow, J Non-Newtonian Fluid Mech 155 (2008), 161–173. [23] B. Namgung, H. Sakai and S. Kim, Influence of erythrocyte aggregation at pathological levels on cell-free marginal layer in a narrow circular tube, Clin Hemorheol Microcirc (2014), doi: 10.3233/CH-141909 [24] J. Newman, Physics of the Life Sciences. 1st ed. Springer, (2008), 231–247. [25] M. Norouzi, M.H. Kayhanm, M.R.H. Nobari and M.K. Demneh, Convective heat transfer of viscoelastic flow in a curved duct, World Acad Sci Eng Technol 56 (2009), 327–333. [26] M. Norouzi, M.H. Kayhani, C. Shu and M.R.H. Nobari, Flow of second-order fluid in a curved duct with square cross-section, J Non-Newtonian Fluid Mech 165 (2010), 323–329. [27] M. Norouzi, M.R.H. Nobari, M.H. Kayhani and F. Talebi, Instability investigation of creeping viscoelastic flow in a curved duct with rectangular cross-section, Int. J Non-Linear Mech 47 (2012), 14–25. [28] M. Norouzi, M.M. Shahmardan and A. Shabani Zahiri, Bifurcation phenomenon of inertial viscoelastic flow through gradual expansions, Rheol Acta 54(5) (2015), 423–435. [29] P.J. Oliveira, F.T. Pinho and G.A. Pinto, Numerical simulation of non-linear elastic flows with a general collocated finite-volume method, J Non-Newtonian Fluid Mech 79 (1998), 1–43. [30] R.G. Owens, A new microstructure-based constitutive model for human blood, J Non-Newtonian Fluid Mech 140 (2006), 57–70. [31] R. Ponalagusamy, R. Tamil Selvi and A.K. Banerjee, Mathematical model of pulsatile flow of non-Newtonian fluid in tubes of varying cross-sections and its implications to blood flow, J Franklin Inst 349 (2012), 1681–1698. [32] G. Pontrelli, Blood flow through a circular pipe with an impulsive pressure gradient, Math Models Methods Appl Sci 10 (2000), 187–202. [33] A.S. Popel and P.C. Johnson, Microcirculation and hemorheology, Annu Rev Fluid Mech 37 (2005), 43–69.
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
121
[34] A.R. Pries and T.W. Secomb, Microvascular blood viscosity in vivo and the endothelial surface layer, Am J Physiol Heart Circ Physiol 289 (2005), H2657–H2664. [35] H.A.G. Rojas, Numerical implementation of viscoelastic blood flow in a simplified arterial geometry, Med Eng Phys 29 (2007), 491–496. [36] D.S. Sankar and K. Hemalatha, A non-Newtonian fluid flow model for blood flow through a catheterized artery—Steady flow, Appl Math Model 31 (2007), 1847–1864. [37] D.S. Sankar and U. Lee, Mathematical modeling of pulsatile flow of non-Newtonian fluid in stenosed arteries, Commun Nonlinear Sci Numer Simul 14 (2009), 2971–2981. [38] A. Shariatkhah, Numerical Investigation of Unsteady Viscoelastic Flows inside the Closed Channels. MSc thesis, Shahrood University of Technology (2013). [39] Y. Sugii, S. Nishio and K. Okamoto, Measurement of a velocity field in microvessels using a high resolution PIV technique. Ann NY Acad Sci 972 (2002), 331–336. [40] G.B. Thurston, Viscoelasticity of human blood, Biophys J 12 (1972), 1205–1217. [41] G.B. Thurston, Rheological parameters for the viscosity viscoelasticity and thixotropy of blood, Biorheology 16 (1979), 149–162. [42] K. Tsukada, H. Minamitani, E. Sekizuka and C. Oshio, Image correlation method for measuring blood flow velocity in microcirculation: Correlation ‘window’ simulation and in vivo image analysis, Physiol Meas 21 (2000), 459–471. [43] H.K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics- THE FINITE VOLUME METHOD. 2nd ed. Pearson (2007), 193-197. [44] S.C. Xue, N. Phan-Thein and R.I. Tanner, Numerical study of secondary flows of viscoelastic fluid in straight pipes by an implicit finite volume method, J Non-Newtonian Fluid Mech 59 (1995), 191–213. [45] S.C. Xue, R.I. Tanner and N. Phan-Thien, Numerical modeling of transient viscoelastic flows, J Non-Newtonian Fluid Mech 123 (2004), 33–58. [46] Y.L. Yeow, S.R. Wickramasinghe, Y.K. Leong and B. Han, Model-independent relationships between hematocrit, blood viscosity, and yield stress derived from Couette viscometry data, Biotechnol Progr 18 (2002), 1068–1075. [47] P. Yue and J.J. Feng, A general criterion for viscoelastic secondary flow in pipes of noncircular cross section, J Rheol 52 (2008), 315–332.
Appendix The components of stress for the Giesekus constitutive equation in axisymmetric cylindrical coordinate system ( ∂θ∂ = 0, vθ = 0) are given as follows:
∂τrr ∂τrr ∂τrr ∂v r ∂v r + vr + vz −2 τrr + τzr τrr + λ1 ∂t ∂r ∂z ∂r ∂z
λ1 2 ∂v r 2 +α τrr + τrz = 2ηp ηp ∂r
∂τzz ∂τzz ∂τzz ∂v z ∂v z τzz + λ1 + vr + vz −2 τrz + τzz ∂t ∂r ∂z ∂r ∂z
∂v z λ1 2 2 +α τrz + τzz = 2ηp ηp ∂z
∂τrz ∂τrz ∂τrz ∂v r ∂v r τrz + λ1 + vr + vz − τrz + τzz ∂t ∂r ∂z ∂r ∂z λ1 ∂v z ∂v r +α τrr τrz + τrz τzz = ηp + ηp ∂z ∂r τzr = τrz
(A1)
(A2)
∂v z ∂v z − τrr + τrz ∂r ∂z
(A3) (A4)
109
Numerical simulation of blood flow through a capillary using a non-linear viscoelastic model Amin Shariatkhaha , Mahmood Norouzia,∗ and Mohammad Reza Heyrani Nobarib a b
Department of Mechanical Engineering, Shahrood University of Technology, Shahrood, Iran Department of Mechanical Engineering, Amirkabir University of Technology, Tehran, Iran
Abstract. In this article, a periodic developing blood flow in a capillary is simulated using a non-linear viscoelastic model for the first time. Here, the Giesekus model is used as the constitutive equation, and based on the experimental data, the best value for the mobility factor and zero shear rate viscosity are derived. The numerical solution of the problem is obtained using the finite volume method. The algorithm of the solution is pressure implicit with splitting of operators (PISO). The simulation carried out using the Giesekus, Oldroyd-B and Newtonian models and the results indicate that the Giesekus model presents a more accurate solution for the stress and velocity fields than the Newtonian and Oldroyd-B models. The previous studies on this problem were restricted to the linear and quasi-linear viscoelastic models. It is shown that only non-linear viscoelastic models can accurately describe the experimental data of unsteady blood flow in capillaries. Keywords: Giesekus model, oscillating flow, blood viscosity, mobility factor, capillary
1. Introduction Clinical science strongly depends on the experiments, which are usually very difficult to set up; therefore, the accurate simulations of human physiological processes can benefit the clinical researchers. However, the numerical results will be reliable if they provide a suitable approximation for the available experimental data. Blood is a combination of red blood cells (RBC), white blood cells (WBC) and plasma [12]. Although, 55% of blood is plasma, which is mostly water (92% by volume) [12] and appears to be a Newtonian fluid, the existence of RBC and WBC make blood very complicated. Arising from the complex characteristics of the blood, it cannot be simulated without considering the effects of those cells on the properties. One of the most important characteristics is related to the red blood cells aggregation which cause the non-linearity and non-Newtonian properties of blood [9, 20, 23, 30, 40, 41]. At large shear rates [1, 8], blood can be assumed as a Newtonian fluid [14, 16–18]. When heart pumps blood, the aorta delivers blood into the vessel network, where the size of vessels decrease but the number of them increases. The total cross sectional area of these small vessels is 1000 times greater than the aorta [12]. It is clear that since the blood volume is constant, by increasing the cross sectional area, the velocity becomes slower. Hence, the low velocity makes shear rate of flow small and the blood inside capillaries can be no longer investigated accurately via Newtonian assumption. ∗ Corresponding author: Mahmood Norouzi, Department of Mechanical Engineering, Shahrood University of Technology, Shahrood, Iran. Tel.: +98 9123726933; Fax: +98 2733335445; E-mails: [email protected] and [email protected].
1386-0291/16/$35.00 © 2016 – IOS Press and the authors. All rights reserved
110
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
Deformation and aggregation of RBCs and also the large diameter ratio of RBCs to capillaries are the reason of some non-Newtonian characteristics of blood flow [30]. Previous studies showed that blood has viscoelastic properties [30, 40]. RBCs in blood are similar to polymeric additives. They make blood a viscoelastic fluid in which temporal effects (relaxation and retardation time) and other mechanical properties (dependency of stress on the deformation in addition to the deformation rate) [5] should be considered. Although many studies considered blood as a generalized Newtonian fluid (GNF) and employed models such as Power-law, Herschel Bulky, Cross and Casson [4, 6, 31, 36–37], studying blood becomes crucial especially in complex flows like stenosis and branching of capillaries, where simple models are unable to describe the physics of flow accurately. Upper Convected Maxwell (UCM) and Oldroyd-B are two quasi-linear viscoelastic constitutive equations that have been used in former studies [32]. Viscoelastic properties of fluids are like a discrete spectrum and various elements should be used to cover the relaxation spectra [5]. Thurston [41] used a generalized Maxwell model which is combination of six Maxwell elements and obtained the constants of this model based on the oscillation test of blood. Hern´an et al. [35] used a generalized Maxwell model similar to the work of Thurston [41]. Owens [30], Gonzalez and Owens [21] and Gonzalez et al. [22] introduced new viscoelastic model based on the temporal network theory which was used previously for modeling the polymeric flows. According to the literature, the previous studies about blood flow are restricted to quasi-linear viscoelastic and generalized Newtonian models. In the present study, the Giesekus constitutive equation (a non-linear viscoelastic model) is employed to simulate the unsteady blood flow in an arteriole for the first time. The empirical investigations indicated that the blood flow shows the non-linear viscoelastic behavior in capillaries. It is important to remember that the non-linear viscoelastic properties of blood flow cause some phenomena: 1. The non-linear viscosity which was previously studied by simple models such as power law. This property causes the blunting of velocity profiles. 2. The first and second normal stress differences. This property could not be modeled via pure viscose and quasi-linear viscoelastic constitutive equations. 3. A phase difference between velocity and stress fields in periodic viscoelastic flows. This effect is happen due to the relaxation spectra of viscoelastic materials. In this paper, the non-linear Giesekus constitutive equation is used to model all of above mentioned effects for unsteady blood flow in a capillary. This model is also able to spot the effect of the third invariant of shear rate in stress tensor and elongational effect. The constants of Giesekus model are obtained by implementing an optimization process on experimental data of rheological tests. Here, the governing equations are solved using finite volume method based on the pressure implicit with splitting of operators (PISO) algorithm. The numerical simulations are also performed via Newtonian and Oldroyd-B model and the results are compared with Giesekus constitutive equation. 2. Methods 2.1. Viscosity of blood Viscosity of blood depends on different parameters such as hematocrit, shear rate and diameter of vessel [11, 20, 24, 34]. Hematocrit of normal human is 45% [12] and its high shear viscosity for vessels larger than 20 m is about 0.0032 pa.s [11, 34]. It is possible to consider the blood as a dilute viscoelastic solution in which the Plasma is the Newtonian solvent and other components play the role 0.001 of viscoelastic additive. Based on this assumption, the viscosity ratio of blood is ηηps = 0.0022 [21, 22,
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
111
Fig. 1. Diagram of blood viscosity versus the shear rate for three rheological models and experimental data [46].
30], where ηs is the plasma viscosity (as the solvent) and ηp is the viscosity of other components (as the viscoelastic additives). In this paper, the constants of Giesekus constitutive equation are calculated optimally based on the non-linear viscosity of blood (see Fig. 1). Here, these constants are obtained using the least square of error between the experimental data and viscometric response of Giesekus model in steady shear test [46]. The results of optimization process indicated that α ≈ 0.5 and η0 = 0.0063pa.s are the best approximations. These are close to the results of similar study done by Hakim et al. [13] to find the constants of Giesekus model. In order to obtain the whole diagram of viscometric function of a material, specific kinds of rheometer should be used at any range of shear rates [7]. However, the former researchers reported only the viscosity of blood using one kind of viscometer [15, 20, 46]. The experimental data which was used in this study is obtained using Couette viscometer [46]. This viscometer is not suitable for very small shear rate [7], so noticeable deviation is expected in this region. According to Fig. 1, it is clear that the Newtonian constitutive equation is not a suitable model for simulating the blood flow at small shear rates where the changes of viscosity are significant, but it has an acceptable approximation at large enough shear rates. The quasi-linear viscoelastic models (such as UCM, second-order and Oldroyd-B models) predict the elastic material modules of blood flow but they considered a constant viscosity for all shear rates in steady shear flows [5]. Therefore, it would be expected that the non-linear constitutive equations simulate the blood flow in capillaries better than the Newtonian or quasi-linear models.
2.2. Governing equations The governing equations consist of the continuity and momentum equations are presented as: ∇ · v = 0
ρ
∂v + v · ∇v = ∇ · T ∂t
(1a) (1b)
where v = (v r , v θ , v z ) is the velocity vector, T , t and ρ are the Cauchy stress tensor, time and density of fluid, respectively. Assuming the body force to be zero and the flow to be incompressible, T is the sum of hydrostatic pressure and extra stress:
112
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
T = −pI + τ
(2)
where τ, p and I are extra stress, hydrostatic pressure and identity matrix. Cylindrical coordinate system (r, θ, z) is used to analyze the flow in a straight circular pipe, where the flow is axisymmetric (∂ ∂θ = 0 , v θ = 0) [47]. Using the axisymmetric condition, the governing equations are simplified as follows: 1 ∂ ∂v z (rv r ) + =0 r ∂r ∂z ∂v r 1 ∂ ∂v r ∂v r ∂ = + vr + vz (rTrr ) + Tzr ρ ∂t r ∂r ∂r ∂z ∂z
∂v z ∂v z ∂v z ρ + vr + vz ∂t ∂r ∂z
=
1 ∂(rTrz ) ∂ + Tzz r ∂r ∂z
(3a) (3b)
(3c)
2.3. Boundary conditions In this study, the inlet velocity is uniform and changes with time. Unlike the previous works, instead of using a cosine function for temporal changes of inlet velocity, the real data of Tsukada [42] is used. There are many experimental data which report velocity profile at some specific time [8, 10, 36] or velocity at some specific point for an interval time [42], but there is no information regarding the changes of mean velocity with time, which is needed in this study for inlet velocity. Therefore, an approximate relation was derived between mean and maximum of velocity profiles [37]. After that in order to have mean velocity versus time, this relation was applied to temporal data of Tsukada [42] which reports changes of maximum velocity in arteriole. At the outlet, the flow is assumed to be fully developed. Also zero gradient condition is used for the velocity and stress components in flow direction. Furthermore, no-slip condition is used for the velocity at walls. On the axis of symmetry, vr and radial gradient of vz are zero. Pressure has zero gradients across all boundaries. The stress field is known at the inlet to satisfy the prescribe inlet velocity. The zero value for τrz and τrθ , and zero radial gradient for other stress components are applied at the symmetry boundary. Also, the zero gradients are used for stress field at walls. 2.4. Initial condition The Giesekus constitutive equation has high order non-linear terms and the convergence of numerical simulation is highly depended on the initial condition. To overcome this problem, the initial condition for the velocity is set as the inlet velocity, while for pressure and stress are set as zero [37, 45]. Due to the periodic behavior of flow field, the effect of initial condition on the final solution is completely waned after some cycles. 2.5. Constitutive equation There are some families of quasi-linear and non-linear viscoelastic models which are obtained based on the different approaches in rheology [2]. The Oldroyd-B constitutive equation is a well-known quasi-linear model. The model is derived by replacing the time derivative of stress by upper convected derivative in Jeffreys model. The Oldroyd-B model is suitable for dilute viscoelastic solutions and its extra stress is usually derived by considering the Newtonian behavior for solvent and UCM model for viscoelastic additives:
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
113
τ = τs + τp
(4a)
τs = ηs γ˙
(4b)
τp + λ1 τp(1) = η0 γ˙
(4c)
where τ is the extra stress tensor, λ1 is the relaxation time, η0 is the zero shear rate viscosity, γ˙ is the shear rate tensor, ‘s’ is the solvent contribution (Newtonian) and ‘p’ is the polymeric contribution (viscoelastic). Here, τ(1) is the upper convected time derivative of stress [5] and can be expressed as τ(1) =
D τ − {(∇v)† . τ + τ. (∇v)} Dt
(5)
D is the substantive derivative, ∇v is the tensor of derivative for the fluid and superscript “†” where Dt means transpose of the tensor. In a steady shear flow test, this model presents the following viscometric functions:
η = η0
(6a)
1 = 2η0 (λ1 − λ2 )
(6b)
2 = 0
(6c)
where 1 and 2 are the coefficients of the first and second normal stress differences and λ2 is the retardation time. The model is linear because the viscosity and first normal stress difference coefficient are revealed to be constant. Also the second normal stress coefficient is constant and zero. In order to have a non-linear model, a non-linear term of stress should be added to the constitutive equation. This modification is implemented in Giesekus model [5]: τp + λ1 τp(1) + α
λ1 {τp . τp } = ηp γ˙ ηp
(7)
The Giesekus model contains four parameters: a relaxation time (λ1 ), the solvent and polymeric contribution of zero shear rate viscosity (ηs and ηp ) and mobility factor (α). The origin of the term involving α can be associated to the anisotropic Brownian motion and/or anisotropic hydrodynamic drag on the constituent polymer molecules [5]. In a steady shear flow test, Giesekus model for a simplified case (λ2 = 0, α = / 0, 1 and γ˙ → ∞) present the non-linear viscometric functions: η∼
(1 − α)/α
√
η0 λ1 γ˙
2 η0 λ1 (α(1 − α))1/4 ˙ 3/2 α (λ1 γ) η0 λ1
2 ∼ − ˙ 2 (λ1 γ)
1 ∼
(8a) (8b) (8c)
Unlike the Oldroyd-B model, the viscometric functions of Giesekus model shows a fractional dependency on shear rate which can be fitted appropriately on experimental data. It should be noted that this formulation is defined for polymeric solutions and extra stress is the sum of polymeric and Newtonian contributions. The polymeric part of solution causes the viscoelastic properties. This is the same as blood in which the RBCs plays the viscoelastic role while Plasma (and other components with small molecular weight) acts as the Newtonian counterpart. The Giesekus constitutive equation can be used in various modes to model the relaxation spectra. Normally, viscoelastic
114
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
fluids have a spectrum of relaxation time instead of only one mode [5]. Some dominant modes can be used if more accurate simulation is required. 2.6. Numerical method In this study, finite volume method is used to solve the governing equations. Algorithm of solution is chosen as Pressure Implicit with Splitting of Operators (PISO). PISO algorithm is a pressure-velocity calculation algorithm, which has one predictor step and two corrector steps and may be assumed as an extended version of Semi-Implicit Method for Pressure Linked Equations (SIMPLE) [43]. It’s possible to use more than two corrector steps for a rapid convergence. All of flow parameters are discretized using central differences, except for the convection terms which are approximated by the linear-upwind differencing scheme (LUDS). This is a generalization of the well-known upwind differencing scheme (UDS), where the value of a convected variable at a cell face location is given by its value at the first upstream cell center. In the linear-upwind differencing scheme, the value of that convected variable at the same cell face is given by a linear extrapolation based on the values of the variable at the two upstream cells. In general, LUDS is the second-order accurate, as compared with first-order accuracy of UDS, and thus, the problem of numerical diffusion would be reduced by using it [29]. The time discretization is carried out using the second order Crank Nicholson scheme. It should be noted that all of the governing equations can be written in the form of the general transport equation as follows: ∂ ∂t
∂ Aφ) + (AUj φ ∂xj
∂ = ∂xj
∂φ ∂xj
+ Sφ
(9)
where the working variable φ is a component of a vector or tensor and even an unknown constant, U is velocity and t is time. The coefficients A and have different meanings for different working variables and Sφ is a source term which includes all the terms that cannot be accommodated in the convective and diffusion terms. This has different contents for different equations. The flow domain is divided into a set of non–overlapping control volumes ∇V with bounding surface area [44]. By Integrating Equation (9) over any control volume and a time increment δt and using the divergence theory, we have: A
V (φp − φpn ) + δt
(Auj φ −
∂φ )nj dA = S¯ φ ∂xj
(10)
where the subscript p refers to the grid point, superscript n denotes the value evaluated at time–level n and S¯ φ is the volume integral of source term Sφ . This can be linearized as: S¯ φ =
Sφ dV = S¯ c + S¯ p φp
(11)
where S¯ c does not explicitly depend on φ and S¯ p is the coefficient of φp , which is made negative to enhance the numerical stability of the discretized equation system. Using a proper spatial variation approximation scheme, the final discretized equation, which relates φp to its neighboring grid point values, can be symbolically expressed in a general form for each control volume: ap φp =
nb
anb φnb + S¯ c + ap0 apn
(12)
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
115
where ‘nb’ means neighborhood and ap0 = A
V , ap = anb + ap0 − S¯ p δt nb
(13)
It is then solved by PISO iterative algorithm. The system of linear equations for pressure is solved by Multi-Grid (MG) algorithm. The convergence criterion of numerical solution is reaching to the residual less than 5 × 10−8 . The system of linear equations for the velocity and stress is solved using Cholesky algorithm with the residual of 5 × 10−7 . The maximum time interval (δt) is set as 5 × 10−8 s for the Giesekus model, 10−6 s for the Oldroyd-B model, and 10−5 s for the Newtonian model. The shorter value of δt for Giesekus model is related to the higher order non-linearity of this constitutive equation. The vessel is simulated as a straight pipe with circular cross section. Due to the axisymmetry condition, a 2D rectangular mesh was used. The ratio of duct length to its radius is considered equal to 40 and its radius is 25 m. The reference of cylindrical coordinate system is located at the left bottom corner of rectangular. The length of pipe is considered long enough to satisfy the fully developed conditions. A structured uniform mesh with rectangular cells is used, where 15 nodes in the r direction and 150 nodes in the z direction are considered. The density of blood is 1053 kg.m–3 [16, 30], its viscosity at large shear rates (η∞ ) is set to 0.0032 pa.s [11, 12, 30, 34], and the relaxation time of blood is assumed to be 0.008 s [20]. 3. Results and discussion 3.1. Grid study and validation In order to check the grid dependency of numerical solution, we used four types of meshes with different grid numbers. Characteristics of these grids are: M1 (5∗ 10 nodes), M2 (10∗ 100 nodes), M3 (15∗ 150 nodes) and M4 (20∗ 200 nodes). In order to have independent results from mesh size, these grids are tested for a known flow and finally one of them is selected to be used for other numerical simulations. The selected grid has an acceptable accuracy while it does not have excessive runtime. The flow of Oldroyd-B model with constant inlet velocity is chosen for the test case. Here, the maximum value of main flow velocity in steady fully developed condition is used to check the accuracy of different grids. The exact analytical solution of this problem indicates that the ratio of maximum value of axial velocity to inlet velocity is 2.0 for steady fully developed Oldroyd-B fluid flow in a straight pipe [5]. The solution of Oldroyd-B flow is obtained by considering the zero value for α parameter of Giesekus constitutive equation in our CFD code. The error of maximum fully developed axial velocity according to the analytical solution and runtime of simulation, in order are: M1 (26.65%, 69 min), M2 (3.9%, 175 min), M3 (0.75%, 238 min) and M4 (0.35%, 491 min). It’s clear that M1 and M2 have considerable errors while the accuracy of M3 and M4 is acceptable (less than 1%). Because the runtime of M3 is less than the half of M4 so M3 is selected as the reference mesh grid structure. It is also possible to check the order of CFD simulation for non-Newtonian flows based on the grid numbers [25–28]. Applying the Richardson extrapolation on data of these four types of mesh, the order of numerical solution is obtained around 2. 3.2. Numerical simulation In this section, the effect of non-linearity of viscoelastic model on periodic blood flow is studied and its results are compared by the results of Newtonian and quasi-linear models which are used in previous studies. The numerical simulation indicated that the Newtonian and Oldroyd-B cases run about 18.4 and
116
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
Fig. 2. Velocity profile in two specific times: (a) Diastole (minimum velocity) and (b) Systole (maximum velocity).
6.2 times faster than the Giesekus case, respectively. In other words, the CFD simulation of unsteady blood flow by Giesekus model is slower and more unstable which is related to the non-linearity of this model. Figure 2 shows the systolic and diastolic velocity profiles (maximum and minimum value of velocity in a period) for different rheological models. According to the figure, the velocity profile of Giesekus model is flatter than the Newtonian and Oldroyd-B fluids. This type of velocity distribution is attributed to the shear thinning behavior of Giesekus model that provides larger values of viscosity in the vicinity of central region of vessel and consequently the small value of shear rate. The Maximum to mean ratio of velocity (a parameter that shows flatness of velocity profile) is approximately 1.5 for blood in experimental studies [10, 36]. This ratio is exactly 2 for Newtonian fluid flow [5]. According to the Fig. 2, this parameter is exactly 2 for Newtonian model (maximum velocity is twice of mean (inlet) velocity) and for the Giesekus model is closer to experimental value so the velocity profile of Giesekus model is shown to fit better with the experimental observations [10, 36] of blood flow in vessels in comparison to other models. According to the study of Popel and Johnson [33], the velocity profile becomes flatter in vessels smaller than 32 m, where the size of RBCs becomes more significant fraction of vessel diameter. Mean difference between results of Giesekus and Oldroyd-B model is 2.1% in Fig. 2a and 1.09% in Fig. 2b. Mean difference between results of Giesekus and Newtonian model is 4.12% in Fig. 2a and 2.24% in Fig. 2b. Figure 3a and 3b shows systolic and diastolic velocity contours for Giesekus model. Here it’s shown again that velocity become fully developed so quickly. According to these figures, it seems that selection of 40 as the value of geometry aspect ratio looks overestimated. Developing region of shear stress of Giesekus model is shown in Fig. 3c and 3d. It’s clear that the developing length is longer for stress and suitable length is selected for the pipe. The stress has a maximum in the upper corner of the entrance. In this region, the inlet velocity faces the wall and has to become zero; therefore high shear rate exists in this region. Figure 4 compares developing length for velocity and stress in centerline for systolic flow. Here, it is shown that for having fully developed flow in outlet, the shear stress is the determinant
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
117
Fig. 3. Contours of velocity for blood flow using Giesekus model: (a) Diastole and (b) Systole. Contours of shear stress for blood flow using Giesekus model: (c) Diastole and (d) Systole.
Fig. 4. Comparing developing length of velocity with shear stress of Giesekus model in systole (maximum velocity).
parameter and not the velocity. Therefore, the blood flow in short vessels has not enough length to get to the fully developed situation and this point should be taken into account in CFD simulations. The history of the Giesekus, Oldroyd-B and Newtonian velocity at the center of pipe in fully developed conditions (at downstream) and also inlet velocity (at upstream) are shown in Fig. 5. The fully
118
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
Fig. 5. History of inlet velocity and velocity in center of vessel for three models.
Fig. 6. History of wall shear stress in center of vessel for three models.
Fig. 7. Profile of Newtonian part, viscoelastic part and total shear stress of Giesekus model in systole.
developed velocities are shown to be larger than the inlet velocity because the growth of boundary layers is completed in this region. According to this figure, the maximum value of Newtonian velocity is exactly twice the size of inlet velocity and is larger than the velocity derived by Giesekus or Oldroyd-B constitutive equations. Algorithm of changes and time of minimum and maximum values are the same for three models. Mean difference between results of Giesekus and Oldroyd-B model is 2.96% and between Giesekus and Newtonian model is 6.93%. One of the most important parameters in studying the cardiovascular systems is the stress of walls of vessel, which is significant in many diseases [8]. Shear stresses on the walls are shown in Fig. 6. The velocity gradients near the walls for the three models are very close, but it’s seen that the wall shear stress is smaller for the Giesekus case. Mean difference between results of Giesekus and Oldroyd-B model is 8.89% and between Giesekus and Newtonian model is 19.32% and seems significant.
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
119
Fig. 8. Diagrams of shear rate and velocity profile versus the radius of vessel for systole.
This event interestingly match with the characteristics of blood flow in small vessels: the presence of RBCs in center of vessel and existence of a layer of plasma near to the wall [3]. In this scenario, blood is divided into two independent parts: RBCs (as viscoelastic additives) and plasma (as a Newtonian solvent). The concentration of RBCs is very low in vicinity of walls so the Newtonian part of stress is the dominant term in this region. However, the Newtonian part of Giesekus wall shear stress (dominant part) is smaller than the pure Newtonian case (wall shear rate is almost equal but viscosity of Newtonian part of Giesekus case is smaller than pure Newtonian case), and the viscoelastic part of stress has a more concentrated distribution in the vicinity of center line of vessel due to presence of RBCs (see Fig. 7). In other words, the viscoelastic portion of Giesekus shear stress obtains a small value around the wall because the shear rate obtains its maximum value and the viscometric function of this constitutive equation is a fractional shear-thinning model. Therefore, it is expected that the total wall shear stress to be smaller in the Giesekus fluids. Figure 8 illustrates the diagrams of shear rate in a cross section versus the radius of vessel. This graph is presented for systole in fully developed region. As can be seen, the shear rate is in order of 103 (s–1 ) at the region near the wall so the Giesekus model can be more reliable while the linear and quasi-linear models can not be suitable choices in this range of shear rate. 4. Conclusions In this paper, the suitability of linear and non-linear constitutive equations in simulation of the unsteady blood flow in capillaries is investigated. It is shown that the Newtonian model is more suitable in the large shear rate cases. On the other hand, at small shear rates, the Giesekus model can more accurately estimate the unsteady blood flow behavior. The Giesekus model gives a flatter profile of velocity, which has revealed a better agreement with experimental observations [10, 33, 36]. The CFD simulations indicate that the results of Oldroyd-B model stay between the Newtonian and Giesekus models. The wall shear stress in the Giesekus model coincides with the other studies concerning with the algorithm of RBCs movement in small vessels. As it is shown the Newtonian model cannot describe stress in small vessels correctly so it recommends utilizing the complicated viscoelastic models. References [1] M. Anand and K.R. Rajagopal, A shear-thinning viscoelastic fluid model for describing the flow of blood, Int J Cardiovasc Med Sci 4 (2004), 59–68. [2] H.A. Barnes, J.F. Hutton and K. Walters, Rheology series - volume 3: An introduction to rheology. 1st ed. Elsevier, Amsterdam (1989), 152–156.
120
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
[3] O. K. Baskurt, H. J. Meiselman, Blood rheology and hemodynamics, In Seminars in thrombosis and hemostasis (Vol 29). New York, stratton intercontinental medical book corporation, c1974-(2003), 435–450. [4] O.K. Baskurt, M. Hardeman, M.W. Rampling and H.J. Meiselman, Handbook of hemorheology and hemodynamicsvolume 69 biomedical and health research. IOS Press, Amsterdam (2007). [5] R.B. Bird, R.C. Armstrong and O. Hassager, Dynamics of polymeric liquids- Volume 1: fluid mechanics. Wiley (1987). [6] J. Chen and L.U. Xi-Yun, Numerical investigation of the non-Newtonian pulsatile blood flow in a bifurcation model with a non-planar branch, J Biomech 39 (2006), 818–832. [7] R.P. Chhabra, J.F. Richardson, Non-newtonian flow and applied rheology – Engineering application. 2nd ed. Elsevier (2008), 5–8. [8] S. Chien, R.G. King, R. Skalak, S. Usami and A.L. Copley, Viscoelastic properties of human blood and red cell suspensions, Biorheology 12 (1975), 341–346. [9] S. Cho, B. Namgung, H.S. Kim, H.L. Leo and S. Kim, Effect of erythrocyte aggregation at pathological levels on NO/O2 transport in small arterioles. Clin Hemorheol Microcirc 59 (2015), 163–175. [10] S. Einav, H.J. Berman, R.L. Fuhro, P.R. DiGiovanni, S. Fine, J. D. Fridman, Measurement of velocity profiles of red blood cells in the microcirculation by laser doppler anemometry (LDA). Biorheology 12 (1975), 207-210. [11] R. F˚ahræus and T. Lindqvist, The viscosity of the blood in narrow capillary tubes, Am J Physiol 96 (1931), 562–568. [12] A.C. Guyton and J.E. Hall, Textbook of Medical Physiology. 11th ed. Elsevier, Saunders (2006), 162–169. [13] S. Hakim, J. Morshedian, M. R. G. Narenji, P. K. Nia, Rheological modelling of Caspian pony blood. Iran Polym J 10(5) (2001), 293–303. [14] Z. Keshavarz-Motamed and L., Kadem, 3D pulsatile flow in a curved tube with coexisting model of aortic stenosis and coarctation of the aorta. Med Eng Phys 33(3) (2011), 315–324. [15] S. Kim, A study of non-Newtonian viscosity and yield stress of blood in a scanning capillary-tube rheometer. PhD thesis, Drexel University (2002). [16] R.E. Klabunde, Cardiovascular physiology concepts. Lippincott, Williams & Wilkins (2011). [17] D.N. Ku, Blood flow in arteries, Annu Rev Fluid Mech 29 (1997), 399–434. [18] M.X. Li, J.J. Beech-Brandt, L.R. John, P.R. Hoskins and W.J. Easson, Numerical analysis of pulsatile blood flow and vessel wall mechanics in different degrees of stenosis, J Biomech 40 (2007), 3715–3724. [19] D.S. Long, M.L. Smith, A.R. Pries, K. Ley and E.R. Damiano, Microviscometry reveals reduced blood viscosity and altered shear rate and shear stress profiles in microvessels after hemodilution, Proc Natl Acad Sci 101 (2004), 10060–10065. ¨ [20] J.A. Long, A. Undar, K.B. Manning and S. Deutsch, Viscoelasticity of pediatric blood and its implications for the testing of a pulsatile pediatric blood pump, ASAIO J 51 (2005), 263–566. [21] M.A. Moyers-Gonzalez and R.G. Owens, A non-homogeneous constitutive model for human blood Part II. Asymptotic solution for large P´eclet numbers, J Non-Newtonian Fluid Mech 155 (2008), 146–160. [22] M.A. Moyers-Gonzalez, R.G. Owens and J. Fang, A non-homogeneous constitutive model for human blood Part III. Oscillatory flow, J Non-Newtonian Fluid Mech 155 (2008), 161–173. [23] B. Namgung, H. Sakai and S. Kim, Influence of erythrocyte aggregation at pathological levels on cell-free marginal layer in a narrow circular tube, Clin Hemorheol Microcirc (2014), doi: 10.3233/CH-141909 [24] J. Newman, Physics of the Life Sciences. 1st ed. Springer, (2008), 231–247. [25] M. Norouzi, M.H. Kayhanm, M.R.H. Nobari and M.K. Demneh, Convective heat transfer of viscoelastic flow in a curved duct, World Acad Sci Eng Technol 56 (2009), 327–333. [26] M. Norouzi, M.H. Kayhani, C. Shu and M.R.H. Nobari, Flow of second-order fluid in a curved duct with square cross-section, J Non-Newtonian Fluid Mech 165 (2010), 323–329. [27] M. Norouzi, M.R.H. Nobari, M.H. Kayhani and F. Talebi, Instability investigation of creeping viscoelastic flow in a curved duct with rectangular cross-section, Int. J Non-Linear Mech 47 (2012), 14–25. [28] M. Norouzi, M.M. Shahmardan and A. Shabani Zahiri, Bifurcation phenomenon of inertial viscoelastic flow through gradual expansions, Rheol Acta 54(5) (2015), 423–435. [29] P.J. Oliveira, F.T. Pinho and G.A. Pinto, Numerical simulation of non-linear elastic flows with a general collocated finite-volume method, J Non-Newtonian Fluid Mech 79 (1998), 1–43. [30] R.G. Owens, A new microstructure-based constitutive model for human blood, J Non-Newtonian Fluid Mech 140 (2006), 57–70. [31] R. Ponalagusamy, R. Tamil Selvi and A.K. Banerjee, Mathematical model of pulsatile flow of non-Newtonian fluid in tubes of varying cross-sections and its implications to blood flow, J Franklin Inst 349 (2012), 1681–1698. [32] G. Pontrelli, Blood flow through a circular pipe with an impulsive pressure gradient, Math Models Methods Appl Sci 10 (2000), 187–202. [33] A.S. Popel and P.C. Johnson, Microcirculation and hemorheology, Annu Rev Fluid Mech 37 (2005), 43–69.
A. Shariatkhah et al. / Numerical simulation of blood flow through a capillary
121
[34] A.R. Pries and T.W. Secomb, Microvascular blood viscosity in vivo and the endothelial surface layer, Am J Physiol Heart Circ Physiol 289 (2005), H2657–H2664. [35] H.A.G. Rojas, Numerical implementation of viscoelastic blood flow in a simplified arterial geometry, Med Eng Phys 29 (2007), 491–496. [36] D.S. Sankar and K. Hemalatha, A non-Newtonian fluid flow model for blood flow through a catheterized artery—Steady flow, Appl Math Model 31 (2007), 1847–1864. [37] D.S. Sankar and U. Lee, Mathematical modeling of pulsatile flow of non-Newtonian fluid in stenosed arteries, Commun Nonlinear Sci Numer Simul 14 (2009), 2971–2981. [38] A. Shariatkhah, Numerical Investigation of Unsteady Viscoelastic Flows inside the Closed Channels. MSc thesis, Shahrood University of Technology (2013). [39] Y. Sugii, S. Nishio and K. Okamoto, Measurement of a velocity field in microvessels using a high resolution PIV technique. Ann NY Acad Sci 972 (2002), 331–336. [40] G.B. Thurston, Viscoelasticity of human blood, Biophys J 12 (1972), 1205–1217. [41] G.B. Thurston, Rheological parameters for the viscosity viscoelasticity and thixotropy of blood, Biorheology 16 (1979), 149–162. [42] K. Tsukada, H. Minamitani, E. Sekizuka and C. Oshio, Image correlation method for measuring blood flow velocity in microcirculation: Correlation ‘window’ simulation and in vivo image analysis, Physiol Meas 21 (2000), 459–471. [43] H.K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics- THE FINITE VOLUME METHOD. 2nd ed. Pearson (2007), 193-197. [44] S.C. Xue, N. Phan-Thein and R.I. Tanner, Numerical study of secondary flows of viscoelastic fluid in straight pipes by an implicit finite volume method, J Non-Newtonian Fluid Mech 59 (1995), 191–213. [45] S.C. Xue, R.I. Tanner and N. Phan-Thien, Numerical modeling of transient viscoelastic flows, J Non-Newtonian Fluid Mech 123 (2004), 33–58. [46] Y.L. Yeow, S.R. Wickramasinghe, Y.K. Leong and B. Han, Model-independent relationships between hematocrit, blood viscosity, and yield stress derived from Couette viscometry data, Biotechnol Progr 18 (2002), 1068–1075. [47] P. Yue and J.J. Feng, A general criterion for viscoelastic secondary flow in pipes of noncircular cross section, J Rheol 52 (2008), 315–332.
Appendix The components of stress for the Giesekus constitutive equation in axisymmetric cylindrical coordinate system ( ∂θ∂ = 0, vθ = 0) are given as follows:
∂τrr ∂τrr ∂τrr ∂v r ∂v r + vr + vz −2 τrr + τzr τrr + λ1 ∂t ∂r ∂z ∂r ∂z
λ1 2 ∂v r 2 +α τrr + τrz = 2ηp ηp ∂r
∂τzz ∂τzz ∂τzz ∂v z ∂v z τzz + λ1 + vr + vz −2 τrz + τzz ∂t ∂r ∂z ∂r ∂z
∂v z λ1 2 2 +α τrz + τzz = 2ηp ηp ∂z
∂τrz ∂τrz ∂τrz ∂v r ∂v r τrz + λ1 + vr + vz − τrz + τzz ∂t ∂r ∂z ∂r ∂z λ1 ∂v z ∂v r +α τrr τrz + τrz τzz = ηp + ηp ∂z ∂r τzr = τrz
(A1)
(A2)
∂v z ∂v z − τrr + τrz ∂r ∂z
(A3) (A4)
