Comparison of Newtonian and non-Newtonian flows in a two-dimensional carotid artery model using the lattice Boltzmann method

Numerical modelling is a powerful tool for the investigation of human blood flow and arterial diseases such as atherosclerosis. It is known that near wall shear is important in the pathogenesis and progression of atherosclerosis. When modelling arterial blood flow it is generally assumed that blood is Newtonian. In this paper, blood flow is modelled in a realistic two-dimensional carotid artery geometry in order to investigate this assumption and its effect on the measurement of near wall shear. The assumption is tested in stenosed and unstenosed geometries and the non-Newtonian blood is modelled using the Carreau–Yasuda model. It is found that the velocity and shear fields, particularly near the walls of the geometries, exhibit small differences in general (<5%) between Newtonian and non-Newtonian models, even in the stenosed geometry with peak differences of 13.6%. Thus, when using numerical modelling to study the haemodynamic influences on atherosclerotic progression, we can safely neglect the non-Newtonian nature of blood.


Introduction
Atherosclerotic cardiovascular disease is a leading cause of morbidity in the industrialized world (Murray andLopez 1996, Caro 2001). There is a body of evidence that suggests a correlation between atherosclerosis, regions of low blood flow velocity, rotational flow and low shear stress near the walls of arteries (Malek et al 1999, Asakura and Karino 1990, Gnasso et al 1997. Therefore, the study of the haemodynamic properties of the blood flow in these regions of the artery can lead to a greater understanding of atherosclerosis and its dependence on flow parameters. However, accurate measurements of quantities of interest, such as shear stress, are difficult to make in vivo, thus numerical simulation becomes a valuable investigative tool. The lattice Boltzmann model (LBM) (Chen and Doolen 1998, Succi 2001, Wolf-Gladrow 2000 uses a simplified kinetic equation to simulate fluid flow and has been applied to many general problems including turbulence (Cosgrove et al 2003), magnetohydrodynamics (Chen et al 1991) and multiphase flows (Shan and Chen 1993). The LBM has also been applied to the simulation of blood flow (Artoli et al 2002, 2003, 2006a, Li et al 2004, Fang et al 2002, Krafczyk et al 1998, Bellerman and Sloot 2001, Boyd et al 2004, Ouared and Chopard 2005, Yi et al 2005 and has been shown to be suitable for modelling a number of features which are important to arterial haemodynamics. In arterial modelling, blood is in general assumed to be a Newtonian fluid. It is held that this assumption is an acceptable approximation for larger vessels, such as the carotid artery (McDonald 1960, Quarteroni et al 2000. However, given the dependence of atherosclerosis on near wall shear, it needs to be asked whether this assumption remains valid for applications studying atherosclerotic progression. LBM modelling of non-Newtonian blood flow is an active area of research. Artoli et al have studied the properties of non-Newtonian blood flow using LBM techniques under a limited number of parameter conditions (Artoli et al 2006b, Artoli andSequeira 2006). Differences between Newtonian and non-Newtonian flow in simple oscillatory flow in a two-dimensional pipe system have also been studied previously (Boyd et al 2007). It was found that differences depended on the Womersley number (α = r √ ωρ/η) and the Reynold number, defined using the stokes layer thickness Re δ = √ 2U 2 0 ρ/ωη , where r is the pipe radius, ω is the angular frequency, ρ is the density, η is the fluid viscosity and U 0 is the peak velocity. Differences in the velocities were seen to increase with decreasing α, and increase slowly with decreasing Re δ . Differences in shear increased with decreasing Re δ and also increased towards an α value of around 4. These differences varied smoothly over Womersley and Reynold's numbers, indicating that any small variations in artery size or pulse profiles in the present study would not significantly change the results. Boyd et al (2007) (using a regular geometry and oscillatory flow) indicated larger differences between Newtonian and non-Newtonian flows for arteries similar or smaller in size to the brachial artery.
This paper will examine the Newtonian blood assumption in an accurate two-dimensional carotid artery model with a particular emphasis on how Newtonian and non-Newtonian velocity and shears differ near the walls of the geometry. The assumption will be examined for an unstenosed and stenosed artery geometry in order to determine whether the presence of a stenosis has any effect. Newtonian carotid artery simulations using the LBM model have been shown (Boyd et al 2005) to be in good agreement with previous numerical (Perktold et al 1991) and physical (Ku et al 1985) models.
In this paper, a two-dimensional model of the human carotid artery with a simulated stenosis growth is presented. The non-Newtonian properties of the fluid will be modelled using the Carreau-Yasuda (C-Y) model (Abraham et al 2006) and the results will be compared to analogous Newtonian flows. The simulations are performed using rigid walls. This is a good approximation in the region of the stenosis (Steinke et al 1994) which is the area of interest here in terms of measuring the wall shear stress. Healthy regions of the artery will exhibit a level of compliance which is not considered. Even in a large artery like the carotid this will have some influence on the calculated velocity and shear stress field; however, given the two-dimensional nature of the simulation and the fact the stenosis shape, size and position varies considerably between patients, the aim of this paper is to investigate how the velocity and shear vary between Newtonian and non-Newtonian models and in the presence of stenosis growth, rather than determining detailed flow fields.

The lattice Boltzmann method
The lattice Boltzmann method (Chen and Doolen 1998) has recently been developed as an alternative method for simulating a range of fluid flow using a modified molecular dynamics approach. In the LBM particle distribution functions, f i (x, t) at point x at time t, are confined to move synchronously on a regular lattice. The distribution functions interact on the lattice in a way that conserves mass, momentum and ensures that the fluid is isotropic and Galilean invariant. Here i labels the lattice link the distribution function is on. The lattice used in this paper is the D2Q9, shown in figure 1. The evolution of the distribution functions on the lattice is governed by the discrete Boltzmann equation (Qian 1992, Chen andDoolen 1998) where for the D2Q9 lattice, see figure 1, e 0 = (0, 0), (i = 0), and i is the collision operator. The fluid density ρ and velocity u can be calculated directly from the distribution functions at each node by ρ = i f i and ρu = i f i e i , respectively. It is assumed that the distribution functions f i can be expanded formally around a local equilibrium distribution such that where ε is a small parameter often taken to be the Knudson number, f The collision operator i is given by the Bhatnagar-Gross-Krook approximation as (Bhatnagar et al 1954, Chen andDoolen 1998) where τ is the relaxation time and f eq i (x, t) is the equilibrium value of the distribution function (Qian 1992, Chen andDoolen 1998). The equilibrium form of the distribution function in two dimensions for the D2Q9 lattice is given by (Qian et al 1992) where w 0 = 4/9, w i = 1/9 for i = 1,2,3,4 and w i = 1/36 for i = 5,6,7,8. The relaxation time τ is related to the viscosity η by The LBM reproduces the Navier-Stokes equation in the nearly incompressible limit and is second-order accurate in the body of the fluid (Chen and Doolen 1998). The stress tensor for an incompressible fluid with pressure p is given by where δ αβ is the Kronecker delta and S αβ = (∇ β u α + ∇ α u β )/2 is the strain rate tensor. It can be shown (Artoli 2003) that S αβ can be calculated locally at each node in the LBM as The f (1) i terms are usually calculated as part of the velocity calculations in the LBM algorithm. Thus, calculating shear in this manner is efficient since it removes the need to calculate derivatives of the velocity. Furthermore, the shear is calculated locally, which is particularly advantageous if the LBM is being implemented in parallel. This method for implementing shear-dependent non-Newtonian flows has been shown to be second-order accurate for a simple power law flow (Boyd et al 2006).
A sub-grid accurate extrapolation boundary scheme (Guo et al 2002) is used to implement the artery geometry in the model. This boundary scheme retains the second order nature of the LBM and is well suited to modelling stenosis growth since it enables the shape of the artery wall and the stenosis to be modelled at a resolution greater than that provided by the underlying lattice (Boyd et al 2004(Boyd et al , 2006. At the entry of the artery a predetermined profile was implemented. This boundary condition was implemented by setting the distribution functions at the entry equal to their equilibrium values, calculated from equation (7), for the desired velocity and density. The profile was applied uniformly across the width of the entry except for a boundary layer region of width approximately 1 mm over which the velocity reduced to zero (Boyd et al 2004). It can be seen from figures 2(a) and (b), which show the carotid artery geometry and the region of the artery presented respectively, that the entry position of the computational artery is a significant distance from the measurement zone. Small variations in the entry parameters, such as the width of the boundary layer, were not found to make a significant difference to the presented results.
The unknown distribution functions at an exit site x were found from a linear extrapolation, based on (Neal 2002)

The non-Newtonian fluid model
The Carreau-Yasuda model (Abraham et al 2005) for the shear thinning behaviour of blood is also commonly used in haemodynamical simulations Veneziani 1997, Gijsen 1998). In this model, the apparent viscosity is given by where a, n and λ are empirically determined constant parameters. The parameters a and n are dimensionless, the parameter λ has dimensions of time.
The shear rate is determined directly from the second invariant of the strain tensor, D II = l α,β=1 S αβ S αβ , asγ = 2 √ D II . Thus using equation (10) the strain, shear rate and apparent viscosity can all be found locally at each grid point. The apparent viscosity is introduced into the LBM by replacing the fixed value of τ in equation (8) with the apparent relaxation time obtained from equation (12).
The main advantage of the Carreau-Yasuda model over other non-Newtonian blood models is that it is continuous for allγ 0. Blood is a shear thinning fluid corresponding to n < 1; thus, we note that limγ →0 η(γ ) = η 0 and limγ →∞ η(γ ) = η ∞ , indicating that at high shear rates, the fluid acts like a Newtonian fluid with viscosity η ∞ , whereas at low shear rates, the fluid acts like a Newtonian fluid with viscosity η 0 . The parameters a, n and λ control how the fluid behaves in the non-Newtonian regime between these two asymptotic viscosities. The continuity of this model at low shear rates allows for an easier implementation in numerical modelling schemes.
In this paper, parameter values of η 0 = 0.1600 Pa s, η ∞ = 0.0035 Pa s, λ = 8.2 s, a = 0.64 and n = 0.2128 were used, obtained from Abraham et al (2005). The density of blood was taken to be ρ = 1 × 10 3 kg m −3 . A dimensionless number Re CY = u 0 Lρ/η ∞ that is analogous to the Reynold's number was also defined.

Methods
Figure 2(a) shows the carotid artery geometry that was used for the simulation, namely the region around the carotid bifurcation. This artery geometry has been used in previous studies (Boyd et al 2005). The common (CCA), internal (ICA) and external (ECA) carotid arteries are indicated in this figure. The unstenosed and stenosed geometries are shown in figure 2(b). The stenosis was placed in a region commonly associated with atherosclerosis in the carotid artery (Caro et al 1971, Zarins et al 1983, Gnasso et al 1997. Figure 2 Holdsworth et al (1999), was implemented at the base of the artery.
Velocity profiles were taken along three lines S 1 , S 2 and S 3 in the artery geometry, see figure 3(b). Simulations were run until the following stop criterion was satisfied: where u(i, t) is the velocity magnitude given by u = u 2 x + u 2 y 1 2 at the ith node along S n at time t, S n is the length of S n in lattice units, T is the period of the pulse cycle and ε was taken to be 1 × 10 −7 . This stop criterion was checked at every 100th of a period, once it was satisfied across each S n simultaneously, the simulations were run for another period and data were output. In general, this convergence criterion was met within 3 to 4 periods. The velocity results will be examined first.

Velocity
The velocity difference between the Newtonian and C-Y velocity was then calculated at every 100th of the period using the following metric: where u N (x, t) and uN (x, t) are the Newtonian and non-Newtonian velocities at time t. N is the number of nodes in the simulation andû N is the peak Newtonian velocity at time t. Velocity profiles were then taken across lines in the unstenosed and stenosed carotid artery geometries (see figure 2(b)) in order to observe differences in the Newtonian and non-Newtonian flows locally in regions of interest in the carotid artery.
Very little difference occurred between the VFFs of the two models at the times examined. The largest difference occurred in the CCA during the low velocity portion of the pulse cycle. The non-Newtonian VFFs in general exhibit lower velocities, but for the majority of the pulse cycle they are very close to the corresponding Newtonian VFFs. Similar results were seen in the stenosed case (data not shown). Figure 4 shows the difference, see equation (14), in velocities between the two models over the whole of the pulse period for the unstenosed and stenosed geometries. The carotid artery velocity waveform is indicated by the dotted line.
There is a sharp decrease in the differences around the peak velocity of the flow. This is the time when the highest shear rates will be present, and the blood viscosity will be approximately Newtonian. The largest differences occur at instances of low velocity in the pulse cycle, with the peak differences occurring at the times of rapid velocity change in the pulse cycle. The difference values for the unstenosed artery are in general smaller than those of the unstenosed artery. The peak difference value in the unstenosed case is V = 0.0232, with an average difference over the period of 0.0141. The peak difference in the stenosed case is V = 0.0198 with an average difference over the period of 0.0133.
In order to observe where the differences occurred in the artery, velocity profiles were taken across the lines indicated in figure 2(b). These profiles were taken at time t = 0.72 s corresponding to the vertical line in figure 4. At this time the differences between the velocities are large, however the same trends were observed throughout the pulse period. Figures 5(a)-(d) show these velocity profiles across the CCA, ECA, ICA and RPB, respectively. It is observed that the largest differences occur in the central region of the profiles, where the non-Newtonian profiles tend to be flatter than the corresponding Newtonian profiles. The profiles agree closely near the walls of the artery, except in the right side of ICA, figure 5(c), where differences occur.
The largest differences are seen in the ICA and RPB. This is true in general for the majority of the pulse cycle. The largest differences are also confined to the central region of the artery during the whole pulse cycle. The peak velocity in the ICA of the stenosed artery, figure 5(c), is larger than that observed in the unstenosed artery. The peak velocities of the profiles in the RPB are similar, but the profile in the stenosed artery, figure 5(d), exhibits some broadening and a shift in the position of the peak velocity towards the right artery wall.
These results support the assumption that Newtonian blood viscosity is a good approximation for the simulation of blood velocity in both the unstenosed and stenosed carotid arteries. Both Newtonian and non-Newtonian blood flow models exhibited physiological features that have been observed previously (Boyd et al 2005).

Shear
All shear calculations were made concurrently with the velocity simulations and were subject to the same stop criterion given by equation (13). The shear difference between the Newtonian and C-Y shear fields was then calculated at every 100th of the period using the following metric: where σ N (x, t) and σN (x, t) are the Newtonian and non-Newtonian |σ xy | shears at time t. N is the number of nodes in the simulation andσ N is the peak Newtonian |σ xy | shear at time t. Shear profiles were then taken across lines in the unstenosed and stenosed carotid artery (see figure 2(b)) geometries in order to observe differences in the Newtonian and non-Newtonian shears locally in regions of interest in the carotid artery.
Both Newtonian and non-Newtonian shears showed accurate physiological characteristics that have previously been observed in a similar Newtonian carotid artery model (Boyd et al 2005). These include persistent regions of low near wall shear along the outer walls of the ICA and ECA near the base of the bifurcation as well as high shear around the bifurcation itself (data not shown). Figure 6 shows the difference, see equation (15), in shears between the two models over the whole of the pulse period. The velocity waveform is indicated by the dotted line, the vertical line for comparison. The shear difference of the pulse period shows the same features as the velocity difference, see figure 4, with a sharp decrease in S around the peak velocity of the flow and increase during low velocity periods of the pulse cycle. Peak differences occur at the times of rapid velocity change in the pulse cycle. The peak difference value in the unstenosed artery is S = 0.0161 with an average difference over the period of 0.0104. In the stenosed geometry the peak difference value was S = 0.0172 and the average difference over the period was 0.0109. Shear profiles were also taken across the lines indicated in figure 2(b) at time t = 0.72 s.
Figures 7(a)-(d) show these shear profiles across the CCA, ECA, ICA and RPB, respectively, for the unstenosed and stenosed, carotid artery geometries. It is observed that the shear profiles match closely in general. The largest differences are seen in the central region ICA for the unstenosed case and in the central regions of the ICA and ECA in the stenosed case, figure 7(c), and near the outer wall of the ICA and the RPB. Similar differences were observed for the majority of the pulse cycle.
Shear profiles in the CCA and ECA of the stenosed artery agree with the corresponding profiles in the unstenosed geometry, figures 7(a) and (b), respectively. Shear near the outer wall of the ICA of the unstenosed artery is smaller in magnitude than the corresponding region of the stenosed artery. Conversely, the shear on the right-hand wall of the RPB shows a sharp increase between the unstenosed and stenosed geometries, corresponding to the region of high near wall shear on the boundary of the implemented stenosis. The near wall shear on the right wall of the RPB shows a small difference between the Newtonian and non-Newtonian models.  figure 7(c) that relatively large shear differences occur there at t = 0.72 s. In both the unstenosed and stenosed cases, the near wall velocity is significantly smaller than the centreline velocity. The largest differences between the Newtonian and non-Newtonian velocities are seen at times of low pulse velocity in the unstenosed and stenosed cases. In the stenosed artery, a large difference is also seen during the time of the steep decrease in the velocity of the driving flow.
Similar results were seen in the near wall shear profile. However, we note that in general the near wall shear is larger than the corresponding shear in the body of the flow, see figures 7(a)-(d).
The relative differences in the near wall velocity and shear were defined as respectively, where u N (t) and σ N (t) are the Newtonian velocity and shear, uN (t) and σN (t) are the non-Newtonian velocity and shear andû N (t) andσ N (t) are the peak Newtonian velocity and shear in the unstenosed and stenosed geometries, respectively. δ u and δ σ for the right wall of the ICA are shown in figures 9(a) and (b), respectively. δ u and δ σ show similar behaviour over the period of the pulse cycle with peak unstenosed differences from 0.28T − 0.36T . Differences in the unstenosed case were generally higher than in the corresponding unstenosed case. The unstenosed near wall shear difference in particular was greater than in the unstenosed case for the majority of the period. It is observed that in general, velocity and shear differences were less than 5% with a brief peak at 13.6%. These percentages are taken as percentages of the peak velocity and shear differences and do not reflect the point-wise difference.

Discussion
These results confirm that the assumption of Newtonian blood viscosity is also valid for modelling the shear fields in a two-dimensional carotid artery model. Shear fields and profiles showed small differences throughout the pulse cycle in both the unstenosed and stenosed artery geometries. Near wall shear in particular also generally showed small differences, indicating that a Newtonian blood viscosity model is also appropriate for studying the near wall shear properties in the presence of atherosclerosis in the carotid artery.
In general, measuring near wall shear is difficult to accomplish accurately in vivo. Being able to model the near wall shear accurately is important, particularly when studying the pathogenesis of atherosclerosis, which is known to be geometrically focal in regions of low near wall shear stress (Malek et al 1999). Accurate modelling can lead to better disease prevention in patients who are at risk of atherosclerosis and can also increase our knowledge of the cellular mechanisms that promote the progression of atherosclerosis. It is known that endothelial cells can react to changes in shear in a variety of ways (Tardy et al 1997), including changes in morphology (Helmlinger et al 1991), proliferation and migration behaviour (Tardy et al 1997, White et al 2001 and modulation of the synthesis and secretion of humoral factors (Galbusera et al 1997, DePaola et al 1999. There is also evidence to suggest that the shape of stenosis plaques and the near wall flow properties play an important role in determining whether a plaque will rupture under the stress produced by changing flow conditions (Richardson et al 1989). It has been suggested that irregularly shaped plaques are more likely to rupture and lead to acute coronary syndromes such as ischaemia (Rothwell et al 2000). Accurate modelling of these physiological situations also has benefits in stent design, arterial graft placement and the planning of angioplastic procedures.

Conclusion
In this paper, the common assumption that a Newtonian approximation of blood viscosity is accurate for a large artery such as the carotid artery has been tested using the LBM. Particular attention has been given to differences in near wall velocity and shear which are known to be significant in the progression of atherosclerosis. This assumption was tested in both unstenosed and stenosed carotid artery geometries. It was found that the differences between the Newtonian and non-Newtonian flows were generally small and occurred mostly in the central regions of the artery. Thus when studying the haemodynamic influences on atherosclerotic progression, we can safely neglect the non-Newtonian nature of blood.