Continuum Scale Non Newtonian Particle Transport Model for Haemorheology -- Implementation and Validation

We present a continuum scale particle transport model for red blood cells following collision arguments in a diffusive flux formulation. The model is implemented in FOAM, in a framework for haemodynamics simulations. Modern mechanistic rheology models are implemented and tested. The model is verified against a known analytical solution and shows excellent agreement for high quality meshes and good agreement for typical meshes as used in vascular flow simulations. Simulation results for different size and time scales show that migration of red blood cells does occur on physiologically relevany timescales on small vessels below 1 mm and that the haematocrit concentration modulates the non-Newtonian viscosity. This model forms part of a multi-scale approach to haemorheology and model parameters will be derived from meso-scale simulations using multi-component Lattice-Boltzmann methods. The code, haemoFoam, is made available for interested researchers.

blood rheology is dominated by the interaction of cells, with a multitude of models having been proposed to account for such meso-scale effects as deformation, aggregation, and rouleaux formation which underline emergent macroscopic flow properties like concentration dependant viscosity and shear thinning. The authors are currently developing a multi-scale approach, explicitly modelling meso-scale effects using Lattice Boltzmann Models (LBM), in which erythrocyte mechanics are fully resolved, while describing the macro-scale rheology using particle transport modelling and quasimechanistic non-Newtonian rheology models. The latter will eventually be parameterised using LBM data. Here, we present the continuum mechanical part of the modelling approach, which allows the simulation of realistic vessel geometries and complex flow patterns.

Particle Migration Model
In a high particle load suspension like blood, many types of mechanical interactions between particles and carrier fluid take place. Mesoscale modelling, using the multi-component Lattice-Boltzmann Method 1 , which has widely acknowledged facility for Lagrangian particulate flows 7,16,17,19,10 is employed to describe these interactions and the dynamics of the collision in detail.
As with direct numerical simulation in turbulence modelling, finite computational resource means that detailed explicit particulate models are limited to small volumes containing relatively few particles in their simulation domain (an the order of magnitude of hundreds to thousands at the time of writing). To address the much greater scales of medical significance, it is, therefore, necessary to develop macro-or continuum scale models, encapsulating the integral effect of these interactions without explicitely resolving them. Crucially, these models must be amenable to parameterisation using mesoscale data, such as 4 . Currently, the models which have been proposed for this task can, roughly, be divided 11 into suspension balance models 18,23 and diffusive flux models 26 .
Suspension balance models use an Euler-Euler mixture modelling approach, where the carrier fluid and the particle load are represented as separate species with a transport equation (typically convection-diffusion) and physical transport properties for each species, while in the diffusive flux models, the suspension is modelled as a single species with the particle volume fraction being modelled as a scalar property, which influences the bulk transport properties.
Our macroscopic model is a particle transport model after Phillips 26 and follows the collision arguments by Leighton and Acrivios 18 . It describes the particle migration based on the gradients of shear strain, concentration and viscosity. The local concentration of haematocrit is then used to establish the local effective viscosity.
A detailed treatise on the rationale behind the compression arguments can be found in Leighton and Acrvios, and Phillips 18,26 , we only give a brief outline at this point.
The transport of haematocrit is dominated by advection -following the bulk flow -variations in concentration are evened by diffusive processes, and the migration within the bulk is driven by a migration pressure. This migration pressure is the result of two phenomena: (1) spatial variation of collision (interaction) frequency, and (2) spatial variation of viscosity.

Spatial variation of collision frequency
Particles that are moving relative to each other in neighboring shear surfaces will experience collisions. The frequency of these collisions is proportional to the shear rateγ, the particle concentration φ, and the particle collision radius a. In a field of constant concentration and constant shear,γφ = const, the collisions are in equilibrium either side of the shear surface, and no net migration will occur. In the presence of gradients of shear rate or concentration, the imbalance of collisions will lead to a "migration pressure" down the gradient. This collsion driven migration pressure can be described as a function of a∇(γφ)). Using a proportinality factor of K c and assuming a displacement proportional to the particle radius a, the migratory flux N c due to variations in collision frequency can be expressed as (using the chain rule):

Spatial variation of viscosity
The displacement of particles after a collision is moderated by viscous effects. In a constant viscosity field the displacement is isotropic and thus balanced with no net migration effects. In a viscosity gradient, the displacement will be less damped in direction of the lower viscosity, leading to a net migration effect down the viscisity gradient.
The displacement velocity is proportional to the relative change in viscosity over a distance that is of order a: a(1/µ)∇µ. With the displacement frequency scaling withγφ, and a proportionality factor of K µ , the migratory flux due to viscosity gradient can be described as (flux is proportional to φ): The scalar transport equation for haematocrit, φ, is then (neglecting molecular diffusion, Brownian motion), where D/Dt is the total differential: with a, particle radius,γ, shear strain rate magnitude, µ, dynamic viscosity, K c and K µ , collision parameters.
Typically, the viscosity is µ = f (γ, φ), which makes the last source term non-linear, which can, in turn, make the solution of this transport equation difficult.
Previous attempts to solve this problem analytically or implement this type of migration model in a numerical model used linearisation of this source term, which involves the derivative of µ in bothγ and φ, and thus limits the model to a specific viscosity model, for which it has been implemented 20,6 . Our current implementation deals with the nonlinear viscosity source term in a way that leaves the viscosity gradient term intact and is thus agnostic to the rheology model used.

Rheology Models
It is obvious from the third RHS term in equation 4, that the particle transport strongly depends on the rheology model it is coupled with. This model implementation aims to be independent of the rheology model. The draw-back of this approach is that errors present in the rheology model, which influence the particle transport, cannot be calibrated out with the parameters of the migration model alone, but the combined set of model parameters will need to be found for any new rheology model that is to be implemented.
Typically, only the shear thinning effects are taken into account, when modelling the non-Newtonian properties of blood in CFD. Common models are of the Carreau and Casson types (REF). In these models, the haematocrit concentration is only used as a bulk parameter in the parametrisation, if at all. Our framework, incorporating the transport of haematocrit, allows the rheological model to take the local particle concentration into account when calculating the local, effective viscosity.
The rheology models that have been implemented and tested in this study are the concentration dependent Krieger-Dougherty model 15 , the Quemada model 27,28,29 with modification by Das 9 (and a new parameter set, which avoids the singularity problem commonly associated with this model), an extended Krieger model, accommodating shear thinning and aggregation effects 12 , a Casson model with haematocrit dependence following Merril et al. 21,9 , and a modified Carreau type model, proposed by Yeleswarapu 32 . All model parameters have been fitted to the experimental data of Brooks 3 ( Figure 1).

Krieger-Dougherty Model
The traditional Krieger-Dougherty model 15 was developed to describe the rheology of high volume ratio suspensions of rigid spherical particles. Rigid, spherical particles do not exhibit shear-thinning behaviour, so the Krieger-Dougherty model is only dependent on the haematocrit concentration φ. It shows a singularity for φ = φ * , where φ * is the haematocrit concentration for which the suspension does stop to behave like a fluid. For rigid spheres φ * = 0.68 15 , while for blood it can go up to φ * = 0.98, which is ususally attributed to the deformability of the erythrocytes 12 .
The parameter n = kφ * is often set to n = 2, but more commonly to the high shear limit of n = 1.82 for φ * = 0.68 28,24 , which is also the value used in this work to allow comparison with the results from Phillips and others 26,20,6 . µ P is the Newtonian viscosity of the liquid phase (plasma).
In this study the Krieger-Dougherty model is not used as for modelling blood viscosity but as a reference model for verification and validation.

Quemada Model
The Quemada model is based on "optimisation of viscous dissipation" 27 . In its original form it is formulated as a Newtonian, concentration dependent viscosity: with k being related to the packing concentration and (for the high shear limit) given as: In this form it is closely related to the Krieger-Dougherty model (eq. 5).
In its non-Newtonian form k is expressed as 28,29 : where k 0 and k ∞ are the intrinsic viscosities at zero and infinite shear, respectively, andγ c is a critical shear rate.
The shear rate magnitudeγ is defined aṡ with D, the symmetric part of the velocity gradient tensor.
Different parameter fits have been proposed for k 0 , k ∞ ,γ c . Cokelet 8,21 proposed: Das 9 noted that Cokelet's parameter set causes the viscosity to be non-monotonous over haematocrit concentration for low shear, and exhibits singularities for zero shear. Das changed the parameter fit for k 0 to which results in a monotonous behaviour for low shear (the lowest shear measured in the Brooks dataset is aroundγ = 0.15 s −1 ), but still shows a singularity for φ = 80.4%. While this is outside the haematocrit values typically encountered in clinical practice, it can still pose a problem if cell migration is taken into account, which will concentrate cells in the core region. In order to overcome this problem, a new parameter set, based on Das's formulation, is derived in this work, which does not show a singularity. Figure 2 shows viscosity over shear rate for low shear rate (γ = 0.15 s −1 ) and zero shear rate. While all the curves show a good fit with the data, the new parameter set does show monotonous behaviour throughout and no singularity below the critical haematocrit.

Modified 5 parameter Krieger Model
Hund et al. 12 proposed and developed a quasimechanistic extension to the Krieger-Dougherty model.
Starting from the traditional formulation of the Krieger-Dougherty model: describing the haematocrit dependence, the shearthinning behaviour is introduced by a variable exponent n: where φ st is the threshold haematocrit concentration below which no shear-thinning is observed. Based on Brooks 3 , this threshold is around φ = 0.15, and n ∞ is modelled using a exponential dependency on φ: Hund's 12 shear-thinning exponent n st comprises contributions of red blood cell aggregation and deformability: where each component is described by a power law: with the empirical coefficient β and ν, and the nondimensional shear rate γ = 1 + (λγ) ν g , as defined by Carreau and Yasuda 30 , with a time constant λ, and ν g = 2. This formulation ensures finite n st at zero shear.
In the 5-component form the aggregation and deformation influences on the shear-thinning exponent are combined into a single power law, due to the limited data on these effects: The model proposed by Hund et al. allows for inclusion of the influence of large molecule concentration (proteins polysacharides, lipids), as well as fibrinogen, and temperature on the constitutive model. Due to a lack of data these are not included in the 5-parameter model.

Yeleswarapu-Wu Model
This model is based on a visco-elastic Oldroyd-B model developed by Yeleswarapu et al. 32,31 . In this study the visco-elastic effects are neglected, only the shear-thinning behaviour and haematocrit dependency are implemented. The shear-thinning behaviour follows a modified Carraeu-type model based on a mixture model by Jung et. al 14 .
The model is based on a mixture model and thus the viscosity is decribed as a function of plasma viscosity µ P and red blood cell viscosity µ rbc 31 : where the red blood cell viscosity is described as: where, in this implementation, k is a constant model parameter, and µ 0 and µ ∞ are modelled as third order polynomials of φ:

Casson-Merrill Model
The Casson model 5 is a classical non-Newtonian model in which the viscosity is modelled as: where µ ∞ is the Casson viscosity (asymptote at high shear rate) and τ 0 is the yield stress. The yield effect means that this model has a singularity at zero shear, leading to infinite viscosity. While there is an argument that blood does exhibit yield at slow time scales and low shear, this effect will typically make this type of model unsuited for numerical simulation within a generalised Newtonian approach with a local effective viscosity due to numerical instability.
For blood, Merill et al. gave the expressions for µ ∞ and τ 0 as 21,9 with the fitting parameters α and β.

Characteristics of Rheology Models
All viscosity model parameters were fitted to experimental data for varying levels of haematocrit in ADC plasma reported by Brooks 3 . While this data is for steady state shear only, it is still considered on of the best datasets for blood rheology data and is used in the majority of work on blood rheology. The parameters were fitted using a Levenberg-Marquardt least squares fit, implemented in Scientific Python (scipy), using the MINPACK library. Table 1 shows the parameter sets for the different models, figure 1 shows the comparison of model results and experimental data. All models show a  The fundamental equations for mass and momentum conservation were implemented using the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) 25 method for steady state, and the PISO (Pressure-Implicit with Splitting of Opera-tors) 13 and PIMPLE (combining PISO and SIMPLE) methods for transient simulations.
Discretision is typically second order in space and time. The code supports all discretisation methods that are supported in the FOAM library (currently foam-extend 4.0 and OpenFOAM 1912).
The haematocrit transport equation 4 is implemented as a scalar transport equation, solved outside of the SIMPLE loop. The Laplacians in φ are implemented implicitely (fvm) as diffusion terms, while the source terms inγ and µ are calculated explicitely (fvc).
For steady state (SIMPLE) and transient cases with the PIMPLE algorithm, underrelaxation is required, typically the underrelaxation factor that is required can be estimated from the order of magnitude of the ratio between collision radius. Stable simulation has been achieved for relaxation factors of 0.1 log(O(a/R)), e.g. a radius R = 50 µm and collision radius of 3.5 µm will require an underrelaxation factor of ≈ 0.1 with no underrelaxation for the final iteration. The PISO algorithm does not use underrelaxation and requires a time step to be estimated from the Courant number (Co < 1) for O(a/R) > 1, and a smaller time step calculated based on a Courant number scaled with the migration velocity.
The discretisation schemes used in the calculations presented in this paper are: second order Euler backward in time and second order (Gauss linear, and Gauss linear upwind for advective terms) in space, gradients are approximated using the least squares theme.
Rheology models are implemented as quasi-Newtonian, with calculation of local cell viscosity based on the shear rate and haematocrit value in the cell from the previous iteration/time step.
The new rheology models that are implemented at the time of writing are the standard Krieger-Dougherty, the modified 5-parameter Krieger, the Yeleswarapu-Wu, and the Quemada model.

Results
All results shown in this paper are for fully developed pipe flow, with periodic boundary conditions between outlet and inlet, with prescribed average velocity. The radius of the pipe varies be- tween 50µm and 5mm, to represent typical vessel diameters. The pipe length is two diameters.

Verification and influence of mesh type
The verification case for the implementation is a pipe of radius 50 µm, average velocity V = 0.0065 m s −1 . The rheology model used in the verification case is the Krieger-Dougherty model to allow comparison to the analytical solution 26 (no analytical solution available for the non-linear terms in the shear-stress and concentration dependent models). Model parameters for the Krieger-Dougherty model are K c = 0.41, K µ = 0.62, φ * = 0.68, and n = 1.82.
The simulation was performed for different meshes, Figure 3, (a) an axisymmetric (2D) wedge with 50 cells resolution in radial direction, (b) a hexahedral, block structured mesh -50 cells radial, and (c) a poyhedral mesh with boundary layer inflation with 60 cells across the diameter -this type of mesh is common in the simulation of vascular flow in patient specific geometries. The given resolutions were chosen based on a mesh convergence study and realistic mesh resolutions as typically used in vascular simulations. The migration model requires a mesh that is of similar resolution as meshes that aim at resolving wall shear stress (WSS) and WSS derived metrics. peaked analytical solution at the axis. The polyhedral three-dimensional mesh also shows good agreement, but the additional numerical diffusion blunts the profile at the axis, the concentration close to the wall is well represented.

Wall shear strain scaling
The parabolic velocity profile for a Newtonian flow is given as: where V is the average velocity.
Therefore, the velocity gradient in radial direction is: So the gradient at the wall (r = R) scales with V and R −1 . The velocity is, therefore, scaled with R, such that the wall velocity gradient is constant. The Reynolds number scales with R 2 . For the given values of R = 0.05, 0.5, 5 mm, V = 0.0065, 0.065, 0.65 m s −1 , the wall velocity gradient is constant atγ w ≈ 650 s −1 , to cover the significant three decades of shear strain magnitude for shearthinning non-Newtonian blood models.
The steady state particle distribution profile is independent of the length scale and the diameter ratio. It will only depend on the ratio of K c /K µ . Figure  5 shows steady state profiles for a range of diameters from 0.1 − 10 mm. The computational effort for the particle migration model, however, scales with R 2 /a 2 , with R, the vessel radius, and a, the particle collision radius. While the small diameter D = 0.1 mm case is fully converged after around 10 4 iterations, the D = 10 mm case requires 10 6 iterations. This corresponds to the diffusion timescales.

Kinematic and particle migration timescales
Blood flow with particle migration is governed by several different time scales for flow kinematics and particle migration. The timescale for the development of the velocity profile (kinematic timescale) is The timescales for the development of the particle migration profile can be derived from the particle migration flux diffusion terms as: The kinematic timescale scales with R 2 /ν, while the particle migration timescales scale with the square diameter ratio R 2 /a 2 , where R is the pipe radius, and a is the particle (collision) radius.
The kinematic viscosity, ν ≈ 3 · 10 −6 m 2 s −1 , while for an average collision radius of red blood cells of a = 3.5 µm, the particle migration diffusion coefficients are of the order of 10 −9 − 10 −11 m 2 s −1 . This means the particle migration happens on timescales that are three orders of magnitude greater than the kinematic timescales. Figure 6 shows the temporal development of the particle distribution and non-Newtonian velocity profile. The flows were initialised with a fully developed parabolic velocity profile and a uniform particle distribution of φ = 0.45 volume fraction. The 0.1 mm case has reached steady state conditions within 0.5 s, the 1.0 mm case shows significant particle migration after physiologically relevant times, while the 10 mm case does show only minimal migration after 10 s. It can be seen that temporal scaling follows the predicted R 2 /a 2 scaling factor.

Variation of rheology model and collision parameter ratio
As is obvious from equation 4, the particle migration is strongly dependent on the viscosity model and the balance between collision and viscosity driven migration, as expressed in the model parameters K c and K µ .
While the magnitude of K c and K µ controls the magnitude of the migration pressures and thus the temporal response of the system, the concentration profile only depends on the balance between collision and viscosity driven fluxes. This balance is expressed by the ratio between the parameters K c /K µ . Figure 7 shows the haematocrit profiles as they develop for different viscosity models -Krieger-Dougherty (K-D), Quemada (Q), Yeleswarapu-Wu Compared to the verification K-D case with K-ratio of 0.66, it can be seen that a shift in the balance to higher influence of the collision frequency (higher K-ratio) steepens the profile, while a lower K-ratio, i.e. a shift of the balance to the resistive influence of the viscosity increase in the low shear region causes a flatter profile.
Comparing the different viscosity models clearly shows the main difference in the core region, where the strong variation in the low shear behaviour, discussed earlier, leads to a strong variation in the relative viscosity gradient (last term in equa-tion{eq:phiTransport}). It is obvious that there is a need for further study and comparison with experimental or meso-scale modelling data to find realistic parameters for each of the potential viscosity models. Especially the modified Krieger model (K5) Based on these preliminary results, the Quemada model with a K-ratio of between 0.5 and 0.6 seems to be the most promising candidate for a semimechanistic rheology model for blood.

Discussion
While previous implementations 20,6 of this class of model are limited by the fact that the viscosity term in equation 4 is linearised in the viscosity gradient with H, our implementation avoids this by implementing the non-linear term directly which allows to use different viscosity models. Our implementation also avoids the use of artificial stabilisation terms that lead to underestimation of RBC migration 2 .
The particle migration time scales with (a/R) 2 , where a is the RBC collision radius. This means that the particle migration is most relevant for small vessels of a diameter of 1 mm or lower, where the migration occurs on physiologically relevant timescales. For larger vessels, minor effects caused by a synergy of particle migration and secondary flows 2 .
The parameters for the migration model would need to be calibrated to experimental data. While such data is available, albeit scarce, for rigid particles in suspension, e.g. based on nuclear magnetic resonance measurements of particle profiles, the authors are not aware of any such data for soft vesicles, in particular RBCs. We therefore hope to use meso-scale models (MCLBM) modeling the cell scale interactions to derive integral diffusion and particle migration measures that can be used to fit the continuum model parameters.
It has to be noted that the implementation uses the magnitude of the shear in the particle flux formulation. As noted by Phillips 26 this assumes an essentially one-dimensional shear state, and isotropic response. This limits the application of the model to flow situations where the shear tensor is aligned with the flow and the main shear in radial direction. As with isotropic turbulence modelling the isotropic migration model will overpredict migration pressure in regions with high anisotropy, e.g. stagnation points, strong acceleration, or rotational shear. It is planned to implement an explicit formulation for a localised, anisotropic shear and migration pressure tensor, similar to approaches proposed by Miller 22 or Fang et al. 11 .

Software
The continuum-scale haemorheology framework was implemented in foam-extend, version 4.0/4.1, and OpenFoam, version 1912. The software (haemoFoam) is freely available to interested parties on github (TS-CUBED/haemoFoam). Please contact the author for testing and developer access.
haemoFoam is a modelling framework for vascular flow simulation based on FOAM, that is intended to cater for the particular requirements of haemodynamics, in particular with respect to WSS related phenomena like atherosclerosis. -Windkessel boundary conditions for outlets viscoelastic rheology models (e.g. Oldroyd B) -platelet transport low density lipoprotein (LDL) transport fluid-structure-interaction (FSI) for flexible vessel walls