Abstract
Shearinduced migration of red blood cells (RBCs) is a wellknown phenomenon characterizing blood flow in the small vessels (micrometre to millimetre size) of the cardiovascular system. In large vessels, like the abdominal aorta and the carotid artery (millimetre to centimetre size), the extent of this migration and its interaction with secondary flows has not been fully elucidated. RBC migration exerts its influence primarily on platelet concentration, oxygen transport and oxygen availability at the luminal surface, which could influence vessel wall disease processes in and adjacent to the intima. Phillips' shearinduced particle migration model, coupled to the Quemada viscosity model, was employed to simulate the macroscopic behaviour of RBCs in four patientspecific geometries: a normal abdominal aorta, an abdominal aortic aneurysm (AAA), a normal carotid bifurcation and a stenotic carotid bifurcation. Simulations show a migration of RBCs from the nearwall region with a lowering of wall haematocrit (volume fraction of RBCs) on the posterior side of the normal aorta and on the lateralexternal side of the iliac arteries. A marked migration is observed on the outer wall of the carotid sinus, along the common carotid artery and in the carotid stenosis. No significant migration is observed in the AAA. The spatial and temporal patterns of wall haematocrit are correlated with the nearwall shear layer and with the secondary flows induced by the vessel curvature. In particular, secondary flows accentuate the initial lowering in RBC nearwall concentration by convecting RBCs from the inner curvature side to the outer curvature side. The results reinforce data in literature showing a decrease in oxygen partial pressure on the inner curvature wall of the carotid sinus induced by the presence of secondary flows. The lowering of wall haematocrit is postulated to induce a decrease in oxygen availability at the luminal surface through a diminished concentration of oxyhaemoglobin, hence contributing, with the reported lowered oxygen partial pressure, to local hypoxia.
1. Introduction
Blood is a biphasic fluid composed of a fluid phase (plasma) and a solid phase composed of blood cells such as red blood cells (RBCs, 6–8 µm in diameter), white blood cells (WBCs, 7–80 µm in diameter) and platelets (PLTs, 2–3 µm in diameter). Blood shows remarkable nonNewtonian effects primarily caused by RBC–RBC interactions and RBC–protein interactions: the aggregation and disaggregation of RBCs and their deformation dictate the value of viscosity [1,2]. A wide body of experimental evidence shows that flowing suspensions of rigid and deformable particles exhibit particle migration [3–10] with a general trend showing migration from regions of higher shear rate to regions of lower shear rate. RBCs also show this behaviour and their migration to the centre of micro and small vessels (micrometre to millimetre size) and the consequent segregation of PLTs near the wall has been the subject of intense study in the past years, both numerically [11–14] and experimentally [15–17]. Different approaches have been proposed to model RBC migration. For example, in [11] a mesoscopic approach employing the immersed boundary method was used to simulate the motion of RBCs, modelled as deformable capsules, flowing through twodimensional channels 20–300 µm in height. The results showed migration of RBCs perpendicular to the wall and formation of a cellfree layer. An extension to the threedimensional case of the same mesoscopic approach is presented in [18]. In [19], the immersed boundary method was also employed to model RBC migration in twodimensional channels of 6–12 µm in height; cell migration from the vessel wall was also observed. In [20], the dissipative particle dynamics method was applied to blood flow in microtubes ranging from 10 to 40 µm in diameter, and cell migration away from the wall to the tube centre was reported. A different approach is represented by continuum models, which are able to simulate larger domains for longer time scales. An example of this approach can be found in [12], where blood flow in vessels with diameter ranging from 40 to 100 µm was investigated using the shearinduced migration model proposed in [21]. Particle migration away from the wall was also observed in this case. RBC distribution influences blood viscosity, and therefore the velocity profile, and also a series of biological activities such as oxygen distribution in the lumen and the scavenging of nitric oxide (NO) by haemoglobin in the RBCs [22]. Arterial wall hypoxia has been proposed to be a contributing factor to atherosclerosis and intimal hyperplasia [23–25], and altered NO transport in arteries has been postulated to induce atherogenesis [22]. Evaluation of the adimensional wall oxygen consumption rate, the Damköhler (Da) number, and comparison with the adimensional mass transport coefficient, the Sherwood (Sh) number, led to postulate that oxygen supply to the arterial wall is fluidlimited [25], therefore underlining the importance of mass transport from and to the wall. It has been shown [26] that the effect of disregarding the presence of haemoglobin, and therefore its possible lowering in concentration due to RBC migration, is a drastic decrease in oxygen transport to the wall. Therefore, given the importance of RBC concentration in the aforementioned series of biological activities, this study aims at evaluating RBC distribution as a result of the synergy between shearinduced migration and secondary flows in four patientspecific geometries, a normal abdominal aorta, a fusiform abdominal aortic aneurysm (AAA), a normal carotid bifurcation and a stenotic carotid bifurcation, and to discuss possible effects on oxygen transport.
2. Material and methods
2.1. Geometry representation and discretization
Four different patientspecific geometries were investigated: Case (A) a normal abdominal aorta, Case (B) a fusiform AAA, Case (C) a normal carotid bifurcation, Case (D) a stenotic carotid bifurcation. The patientspecific geometries were reconstructed from computer tomography angiography (CTA) scans acquired at Karolinska University Hospital, Stockholm. Standard CTA scans, details in table 1, of the abdominal aorta and of the carotid bifurcation were obtained with a 64slice CT machine (Lightspeed VCT, General Electric). The luminal geometries were reconstructed using active contour (deformable) models with the diagnostic software A4clinicsResearch Edition (VASCOPS GmbH, Graz, Austria). For more details, please refer to [27]. The investigated geometries are shown in figure 1, the basic determinants are reported in tables 2 and 3 and the bounding box sizes are reported in table 4.
2.2. Mathematical models and solution procedure
2.2.1. Fluiddynamical model
Particle migration in inhomogeneous shear flows can be captured by modelling the effects of spatially varying interparticle interaction frequency and spatially varying viscosity. This work is based on Phillips' model [21] for shearinduced migration of rigid particles in concentrated suspensions. Blood is considered to be a concentrated monomodal suspension of solid particles, RBCs, in a Newtonian incompressible fluid, plasma. RBCs are modelled in a Eulerian sense through their volume/mass fraction. RBC density is different from the density of plasma but since in the present cases the gravitational field is neglected, buoyancy is absent.
The equations reported below are equivalent to the extension of Phillips' model [21] to the case of solid and fluid phases with different densities developed in [28]. The commercial software COMSOL Multiphysics, used to solve the set of equations, employs a formulation based on the slip velocity u_{slip} between the fluid and the solid phase [29], while the shearinduced migration model proposed in [21,28] is formulated in terms of particle fluxes. This difference in formulations required some manipulation of the original equations with the slip velocity u_{slip} being formulated as a function of the particle fluxes. Indicating with ϕ and c_{s}, the volume and mass fractions of the solid phase, with ρ_{f} and ρ_{s} the constant densities of the fluid and solid phase, respectively, with u_{slip} = u_{s} − u_{f} the slip velocity between the fluid and the solid phase and with u_{migr} = (1 − c_{s})u_{slip} the RBC migration velocity, the continuity equation (the full derivation is given in appendix A) reads [28,29] 2.1where the mixture velocity (massaveraged velocity) is given by 2.2with the mixture density (kg m^{−3}) defined as 2.3The momentum equations take the form [28,29] where u_{slip} is expressed as a function of the particle flux J_{s} (kg m^{−2}s^{−1}) as 2.5with the particle flux defined according to Phillips' model [21] 2.6where denotes the scalar shear rate with D = grad_{S}l being the symmetric part of the velocity gradient tensor l [30], and K_{c} and K_{μ} phenomenological constants of order unity. In this work, K_{c} was set equal to 0.41 and K_{μ} to 0.62 as used in, for example, [12,21].
Finally, the transport equation for the solidphase volume fraction reads 2.7where the solidphase velocity is given by 2.8A list of all variables and parameters used in the model is reported in table 5.
The model also assumes that the Péclet number , with a being the particle radius, the scalar shear rate introduced above and D the Brownian diffusion coefficient, is large due to a very small Brownian diffusion—which is therefore neglected in this model. Both phases share the same pressure field [29].
RBC migration leads to the formation of a plasma skimming layer, a thin layer of around 3 µm adjacent to the wall devoid of RBCs [31]. The continuum approach to model RBC migration can capture the macroscopic haematocrit distribution but not the plasma skimming layer, which is inherently a microscopic phenomenon. For this reason, the wall haematocrit distribution presented in this work must be considered as the haematocrit value at the inner boundary of the plasma skimming layer.
Boundary conditions in the form of an inlet mass flow rate and an outlet pressure waveform, taken from [32] for the normal abdominal aorta and from [33] for the AAA, were applied. For the carotid simulations, flow rates taken from [34] were applied at the common carotid artery (CCA) inlet and at the external carotid artery (ECA) outlet (figure 2). At the internal carotid artery (ICA) outlet, a constant pressure was applied. Additionally, for all cases, a no viscous stress condition, τn = 0, was applied at the pressure outlets. At time t, the prescribed mass flow rate m(t) is related to the inlet/outlet velocity by 2.9with S and n defining the inlet/outlet surface area and the outward normal unit vector, respectively. The noslip boundary condition, u_{w} = 0, was applied at the walls. Regarding the solid phase, at the inlet a fixed volume fraction of ϕ = 0.45, representing a normal haematocrit value, was set for all cases. The use of a pure convective boundary condition at the outlet induced severe oscillations of the volume fraction of the solid phase, forcing to impose a fixed volume fraction value of ϕ = 0.45; for an analysis regarding this assumption refer to the §4. At the wall, a noflux condition, −n · ϕu = 0, was applied. In the case of rigid wall and incompressible fluid, the imposition of a pressure wave at the outlet is not strictly necessary since given the mass flow rate the pressure differential (Δp) is calculated. This means that a constant pressure can also be applied, as done for the carotid artery case. However, in order to get a meaningful pressure distribution (amplitude and phase), a physiological pressure wave needs to be used. A correct pressure boundary condition is also required in cases like fluid–structure interaction simulations. For a detailed explanation of different boundary conditions, the reader is referred to [35].
2.2.2. Constitutive modelling of blood
Blood has complex rheological properties involving shearthinning, thixotropy and viscoelasticity [36,37]. Even in large arteries, the dependence of viscosity from the shear rate has been shown [38] to produce substantial differences in the flow pattern compared with a Newtonian viscosity model. In the present context, the viscosity is considered a function of shear rate and haematocrit and blood is modelled as a structured fluid [39]. To this end, the Quemada viscosity model [40] was used to represent this dependency. Hence, the viscosity reads 2.10where q is an empirical parameter, k_{0} and k_{∞} the intrinsic viscosities for and , respectively, and a critical shear rate related to the critical shear stress under which the structure is broken. The values used in the present work are taken from [14] and summarized in table 6. The viscosity behaviour as a function of shear rate and haematocrit is shown in figure 3.
2.2.3. Finiteelement discretization and numerical algorithms
Lagrange P_{2}P_{1} elements were employed to discretize, with quadratic interpolation, the velocity field and, with linear interpolation, the pressure field. Lagrange linear elements were used for the volume fraction of the solid phase. The Galerkin method was used to discretize the equations with anisotropic diffusion with δ_{ani} = 0.25 for the momentum transport and δ_{ani} = 0.5 for the dispersed phase transport [29].
2.2.3.1. Timeadvancing algorithm
The generalizedα method [41] with the numerical parameter ρ_{∞} = 0.75 was used the as timeadvancing algorithm. An adaptive time stepping approach was employed [29] to ensure a proper resolution of the transient flow field. The equations were solved using an affine invariant form of the damped Newton method [42] with the use of a segregated approach where the first segregated step solved for the velocity and pressure fields while the second segregated step solved for the solidphase volume fraction. The PARDISO direct solver was used to solve the arising linear system of equations.
All computations were performed on a 64bit Dell machine equipped with 4 four cores (16 cores in total) Intel(R) Xeon(R) CPU E78837 2.66 GHz with 512 GB of RAM, with Windows Server 2008 R2 Enterprise.
3. Results
3.1. Model validation
To validate the code used in this work, a series of comparisons with previous published results were conducted; a representative one is reported below.
The used benchmark corresponds to the pipe flow studied in [43]: it consists of a pressuredriven axial flow of a monomodal suspension through a 122 cm long pipe with radius 2.54 cm. A uniformly mixed suspension of spherical particles of 3178 µm in diameter with initial volume fraction of 0.50 was set in motion under an inlet mass flow rate of 0.3 kg s^{−1} with a constant inlet volume fraction of 0.5. Viscosity followed the relation proposed by Krieger [44] 3.1where ϕ_{max} is the maximum packing volume fraction. Table 7 summarizes the parameters and data of the simulation. Only one quarter of the tube was simulated with a quadrilateral mapped mesh with a boundary layer mesh consisting of 10 layers of increasing thickness. The mesh consisted of 21 200 elements, resulting in a total of 582 322 d.f. It is important to underline that for this case the steadystate solution is only a function of the inlet volume fraction, the maximum packing concentration and the ratio K_{c}/K_{μ}. Figure 4 shows the volume fraction and the streamwise velocity in the radial direction at the pipe outlet. Good agreement with the results of [43] is achieved, even though the present solution shows a more pronounced migration with lower values of volume fraction at the pipe outer wall and a smaller slope at the centre of the tube. The normalized velocity field also showed slight discrepancies in line with the volume fraction distribution. Several factors might be the cause of this discrepancy; a possible cause can be the difference in the stabilization method employed. This set of equations is known to be prone to instability [45] and it has been found that the consistent stabilization method Galerkin least squares [46] was not able to stabilize the equations, i.e. resulting in pronounced oscillations in the solution. For this reason, anisotropic diffusion was used instead. Anisotropic diffusion is an inconsistent stabilization method, which might lead to a slightly different solution than the one of the original governing equations. For a discussion regarding the presence of the cusp in the solidphase volume fraction curve at the pipe centre refer to [43].
3.2. Role of secondary flows
Secondary flows play a major role in RBC migration. In order to elucidate their role, before tackling the more complex patientspecific cases, a comparison between the steady and laminar flow in a straight pipe, where no secondary flows are present, and the steady and laminar flow in a curved pipe where welldefined secondary flows, in the form of Dean vortices, are present was conducted. The data relative to the two pipes, representative of a CCA, are reported in table 8. While the Reynolds number used here ensures that the flow is laminar, the value of the Dean number ensures the presence of two stable Dean vortices in the crosssection of the curved pipe.
In figure 5, a comparison between the wall haematocrit in the straight pipe and in the curved pipe, inner and outer curvature walls (see figure 6 for nomenclature), is reported. The straight pipe's wall haematocrit values lie in between the ones of the curved pipe, which shows higher values for the outer curvature side and lower values for the inner curvature side. This result shows the effect of secondary flows, which consists of moving the solid phase from the inner curvature side and accumulate it on the outer curvature side. The following mechanism is postulated to take place: at first, the strong shear layer in the nearwall region induces particle migration thus slightly lowering the haematocrit on both sides, then secondary flows kick in and convect particles away from the inner curvature side region further lowering the haematocrit value. This mechanism is responsible for the haematocrit pattern observed in the patientspecific cases (§3.3).
In figure 6, the solidphase velocity vectors are plotted on top of the solidphase volume fraction showing the Dean vortices in the crosssection of the curved pipe. The wall haematocrit is lower on the inner curvature side owing to the convective motion of the solid phase induced by the Dean vortices.
3.3. Patientspecific arteries
3.3.1. Normal aorta
The normal aorta studied in this work presents a prominent anterior–posterior curvature which leads the posterior wall, the inner curvature wall, to experience a slightly larger migration compared with the anterior wall. The minimum wall haematocrit, of around 0.42–0.43, was observed around t = 0.2 s (peak systole). A slightly larger migration is observed in the iliac arteries, particularly in the right iliac, where a minimum haematocrit of approximately 0.41 is observed at t = 0.30 s (figure 7a). The lateralexternal wall, the inner curvature wall, of both iliac arteries consistently shows a lower haematocrit than the lateralinternal wall, the outer curvature wall (figure 7a). Notably, the low haematocrit spots on the inner curvature wall of the iliac arteries are one of the locations where atherosclerotic lesions are usually observed [47]. RBC migration is caused by two phenomena: the shear layer formed in the nearwall region inducing initial RBCs migration from the wall to the core flow and the secondary flows directed away from the inner curvature wall enhancing RBC transport from the inner curvature wall to the outer one (figure 7b). The time scale of RBC migration is longer compared with the fluiddynamical scales; this can be observed in figure 7a where, despite ample variations in the flow field, the wall haematocrit distribution is fairly constant throughout the cardiac cycle.
3.3.2. Fusiform aneurysm
In the fusiform AAA, no important RBC migration is observed. The proximal portion of the aorta on the inner curvature wall (figure 8) presents wall haematocrit values of around 0.43. The aneurismatic bulge, instead, presents wall haematocrit values of around 0.45, the prescribed value, throughout the cardiac cycle indicating that the shear rate values reached in the boundary layer are not sufficient to induce RBC migration (figure 8). It is important to underline here that even with the presence of strong secondary flows, the presence of shearinduced migration is indispensable to create a lower haematocrit region near the wall which can then be convected. In this case, despite the presence of secondary flows, no sufficient nearwall shearinduced migration, except for the proximal region and a few spots in the iliac arteries, is observed and therefore the secondary flows have no effects on the haematocrit value. Two spots of decreased haematocrit are observed on the lateralexternal wall of both iliac arteries throughout the cardiac cycle (figure 8), consistent with what observed in the normal aorta case.
3.3.3. Normal carotid bifurcation
In the normal carotid bifurcation, a clear RBC migration is seen at the outer wall, the inner curvature wall, of the ICA sinus and on the inner curvature wall of the CCA (figure 9). Wall haematocrit values at these locations show approximately a 10% decrease with respect to the normal value of 0.45. Secondary flows convect RBCs away from the CCA inner curvature wall and shift the low wall haematocrit region from the centre of the CCA towards the ICA sinus. This is illustrated in figure 10, where the RBC inplane velocity vectors at t = 0.30 s are shown together with the wall haematocrit at t = 0.40 s. A timeshift of 0.10 s was adopted since the haematocrit distribution evolves from the RBC velocity field history. In the proximal CCA region (figure 10c), the secondary flows induce a velocity field that moves RBCs from the inner curvature wall towards the outer curvature wall; in the middle CCA (figure 10b) and in the carotid sinus (figure 10a) the secondary flows induce a velocity field that moves RBCs towards the ICA sinus. The origin of this complex flow behaviour is the presence of streamwise vortices running through the whole carotid length (figure 10). Similar vortical structures (VSs) are present in the other analysed cases (not shown). It is worth noting that the pattern of low wall haematocrit correlates with sites where atherosclerotic plaque formation is reported [49,50].
3.3.4. Stenotic carotid bifurcation
The stenotic carotid bifurcation presents low haematocrit values in the stenotic throat caused by the high values of that lead to a pronounced RBC migration (figure 11). The wall haematocrit in the throat region reached values of approximately 0.35, a decrease of around 22% with respect to the normal value. In the carotid sinus, a streak of low haematocrit values on the outer wall, the inner curvature wall, caused by the secondary flows, is observed throughout the cardiac cycle (figure 11). The inhomogeneities in wall haematocrit observed in the ICA and ECA are likely the result of the turbulent nature of the flow in those particular regions (see the Discussion section). The spot of low haematocrit present on the ICA sinus inner wall (figure 11) only at t = 0.1 s is caused by the strong impinging systolic jet originating in the throat.
4. Discussion
Oxygen availability to the underlying vessel wall is linked to wall haematocrit values [26]; considerations on the Sh and Da numbers led to postulate that oxygen supply to the arterial wall is fluidlimited [25], hence the need to evaluate RBC migration in large vessels. Consequently, this work investigated the shearinduced migration of RBCs and its relation with secondary flows in the abdominal aorta and in the carotid artery using Phillips' model [21]. The present results show the importance of secondary flows in the RBC migration process in the macrovasculature. In regions of marked curvature like the CCA, the ICA sinus, the iliac arteries and for this particular case also the abdominal aorta, secondary flows convect RBCs from the inner curvature wall to the outer curvature wall. Shearinduced migration induces the initial lowering in wall haematocrit values while secondary flows enhance or dampen the migration (figure 13). In the stenosis, instead, the drastic reduction in haematocrit is due to the strong shear layer formed in the throat. A correlation emerged between regions of RBC migration and regions prone to atherosclerosis, namely the external wall of the iliac arteries, the external wall of the ICA sinus and the carotid stenosis, being that strong abdominal aortic curvature is not a common occurrence. The present results show that shearinduced migration, and by direct consequence the overall transport, of RBCs loses strength during the enlargement process of the abdominal aorta that leads to the formation of AAAs (compare the normal aorta case with the AAA case), while during the progression of atherosclerosis, and therefore stenosis, RBC migration increases its magnitude, possibly inducing or strengthening hypoxic conditions in the stenotic area, which in turn might lead to further atherosclerosis progression. The decreased wall haematocrit leads to a reduced oxyhaemoglobin availability in the nearwall region which can result, following [26], in a reduced Sh number thus reinforcing hypoxic conditions. The present results are in line with what was observed in [51], where oxygen mass transport in a carotid bifurcation was simulated considering only free oxygen dissolved in plasma. A lowered Sh number on the outer wall of the carotid sinus, the site where atherosclerotic plaques are usually observed, was reported. This decrease in Sh number was attributed to the secondary flows established in the carotid sinus region, which were directed from the outer wall (inner curvature wall) to the inner wall (outer curvature wall) of the carotid sinus, convecting oxygen away from the former to the latter. Also in [52], a decrease in Sh number along the outer wall of the carotid sinus due to recirculating separated flow was observed. This led the authors to postulate that mass transfer of various molecules (e.g. oxygen, NO, mitogens) can be impaired at the outer wall of the carotid sinus, and therefore relating to regions of intimal thickening and atherosclerosis. The same curvature effect was also reported in [53] for a coronary artery model.
Decades of intense research on the role of wall shear stress (WSS) in atherosclerosis were able to unravel its effects on the endothelial cell layer (e.g. [54]). Despite the clear role played by WSS in atherosclerosis, no conclusive evidence of a direct link between low and oscillatory WSS and this pathology [55] has surfaced. WSS, on the other hand, is a manifestation of a more global phenomenon, the flow field as a whole. Focusing on wallmetrics such as WSS, oscillatory shear index (OSI) and relative residence time (RRT) led to the conclusion that curvatureinduced secondary flows have a beneficial flow stabilizing effect [56]. On the other hand, this work shows the negative aspects of these secondary flows, which cannot be captured by wallmetrics alone. In this view, coupling wallmetrics (WSS, OSI and RRT) with global flow field structures (secondary flows) in the analysis process might prove beneficial. A possible way to do so is to improve the haemodynamically inspired geometric variables approach [57], by taking into consideration global flow structures in the definition of the geometric variables. Future studies should therefore analyse the colocalization of regions of low/high WSS and low/high mass transport of oxygen and other important atherogenic molecules such as LDL.
Infrarenal aorta and CCA flow rates, ICA/ECA and iliac flow partitions are patientspecific. On the other hand, proximal CCA and bifurcation geometry are the primary responsible for flow characteristics [58,59], meaning that the carotid flow features investigated in this work are strongly geometry dependent. The same assumption is likely to hold valid also for aortic flows. For this reason, the same secondary flows characteristics are most likely to be observed also under different boundary conditions. The open points and limitations of this study that should be addressed in future research are presented below.
4.1. Shearinduced migration model

— Phillips' model was developed for rigid particles while RBCs are deformable; the effect of deformability and loss of symmetry of RBCs might affect RBC migration.

— The best values for K_{c} and K_{μ} could not be assessed due to a lack of experimental data; values from literature were used instead.

— Phillips' model was developed for onedimensional flows, its extension to threedimensional flows employs the scalar shear rate . The validity of this approach cannot be fully quantified owing to the lack of experimental data.
4.2. Uncertainty on the inlet ϕ value
The inlet distribution of ϕ has been considered uniform and equal to 0.45; in reality, a nonuniform distribution is expected, with lower values of ϕ at the wall. Approximation of the ϕ distribution as a constant introduces a development length (conceptually similar to the length needed to develop a parabolic flow profile) that might affect the downstream wall haematocrit distribution and possibly leading to higher values of ϕ than would occur in reality.
4.3. Numerical challenges
Phillips' model showed a tendency to develop oscillations in the values of the solidphase volume fraction, an issue already reported in literature (e.g. [60]). Anisotropic diffusion, an inconsistent stabilization method, was able to partially solve this issue: in the outlet region and in localized regions of the analysed geometries, the presence of uncontrolled oscillations was still observed. At the outlets of the domains, a fixed volume fraction has been imposed due to uncontrolled oscillations in its value if a pure convective boundary condition was employed. This assumption, while imposing a constraint on the volume fraction field, is not considered restricting based on the fact that the bulk flow presents a constant haematocrit value of ϕ = 0.45 while the differences in haematocrit are confined into a small boundary layer leading to an average value, over the outlet section, of ϕ close to 0.45. Moreover, as can be observed in figure 5, the effects of this boundary condition are segregated in a very narrow region close to the outlet, therefore not affecting the domain of interest. Lastly, the use of an inconsistent stabilization method is one possible cause of the discrepancy observed in the model validation example.
4.4. Turbulence
All simulations assumed laminar flow. While this assumption is justified in the larger vessels, such as the normal abdominal aorta and the AAA, where the Re is low [38], in the normal carotid artery and, in particular, in the stenotic carotid artery turbulent flow can occur [61,62], as observed in the wall haematocrit pattern for the stenotic carotid artery case (figure 11). Therefore, in future work, a proper treatment of the turbulent flow should be applied for these cases to evaluate the effect of turbulence on RBCs transport.
4.5. Oxygen transport
This work assesses low wall oxygen availability by locating regions of low wall haematocrit. Despite being shown to be a valid approach [26], a complete assessment of wall oxygen flux/consumption would require to also consider oxygen transport.
5. Summary
In conclusion, the work presented here analyses the combined effect of nearwall, shearinduced migration and secondary flows on the distribution of wall haematocrit. Shearinduced migration was modelled using a phenomenological model taking into account gradients in scalar shear rate, haematocrit and viscosity. The synergistic effect of these two flow phenomena leads to lowered haematocrit values in regions shown to be prone to atherosclerosis. These lowered wall haematocrit values were shown [26] to induce a lowered oxygen wall flux, which in turn might induce wall hypoxia.
Funding statement
This work has been financially supported by the Project grant no. 20104446 from the Swedish Research Council.
Acknowledgements
The authors would like to thank Dr Anders Ekerot from COMSOL AB for his invaluable technical support and Andrii Grytsan for providing a routine to import COMSOL data into Tecplot.
Appendix A
The continuity equation (2.1) can be derived as follow. First, equation (2.8) is derived starting from equation (2.2) and recalling the definition of the slip velocity u_{slip} = u_{s} − u_{f} so to obtain A1This equation gives A2Recalling the definition of the mixture density, equation (2.3), and substituting A3which gives equation (2.8) A4Starting now from the general continuity equation for a mixture [63] A5and manipulating the first l.h.s. term, one obtains A6In order to eliminate the time derivative, a particle conservation equation formulated using the mixture velocity can be employed. This equation takes the form A7which when substituted in equation (A 6) gives A8The particle flux J_{s} is the particle flux induced by the shearinduced migration, hence in its formulation u_{migr} should be employed. The flux takes then the form A9
Substituting back one obtains A10
The second l.h.s. term in equation (A 5) can be written as A11
Putting the two terms together, one finally gets the continuity equation (2.1) A12
The same approach can be used to recover the momentum equations (2.4).
 Received April 16, 2014.
 Accepted April 30, 2014.
 © 2014 The Author(s) Published by the Royal Society. All rights reserved.