Numerical simulations of blood cell flow in non-Newtonian fluid

We investigated how non-Newtonian viscosity behavior affects the flow characteristics of blood cells. Our findings offer insight about how shear thinning affects the dispersion of liposome-encapsulated hemoglobin and red blood cells in blood. The lattice Boltzmann method was used for fluid calculations, and the rheological properties of the non-Newtonian fluid were modeled with power-law relationships. The deformable three-dimensional red blood cell model was applied. First, we investigated the effects of shear thinning on the flow behavior of single blood cell. Simulation results indicate that shear thinning promotes the axial concentration of red blood cells. Next, varied the hematocrit to see how mutual interference between blood cells affects flow. At low hematocrit, shear thinning clearly promotes the axial concentration of red blood cells. As the hematocrit increases, in contrast, mutual interference has a greater effect, which counteracts shear thinning so the red blood cell distribution resembles the distribution within a Newtonian fluid.


Introduction
Red blood cells (RBCs) play an essential role in supplying oxygen to all tissues in the body. When the diameters of RBCs and vessels are on the same order of magnitude, the rheological properties of RBCs greatly influence blood flow. RBCs are biconcave and diskshaped, and they deform passively under external force. This deformability has some effect on RBCs' tendency to concentrate along the axis of flow paths, which helps them to flow smoothly through narrow microvessels.
Several recent numerical studies have focused on the behavior of deformable RBCs in microvascular flows [1][2][3][4][5][6][7][8][9]. The motion of RBCs is described well with existing numerical techniques. However, most studies have focused on the flow of RBCs in the Newtonian fluid of blood plasma. Currently, liposome-encapsulated hemoglobin (LEH) is being promoted vigorously in Japan as synthetic alternative to RBCs [10]. Typical LEH particles are 200-250 nm in diameter, which is 1/30 smaller than a human RBC. A solution of plasma and LEH shows non-Newtonian behavior in some conditions [11]. To better understand this behavior, we have simulated the flow of surrounding plasma-LEH solution on the flow of erythrocytes in microvessels. The simulation results offer valuable insights into the mechanisms by which the flow of RBCs delivers oxygen through microvascular networks.

Numerical method
We constructed a three-dimensional microvascular model and investigated the flow behavior of RBCs in the Newtonian and non-Newtonian fluid. Figure 1 shows schematics of the microvascular model we employed, in which the vessel diameter (D) is 20 µm and the vessel length (L x ) is 40 µm. The mean velocity in the vessel U m was set to 1.5×10 -3 m/s. Three hematocrit levels were simulated, at 3.0 %, 6.0 % and 15 %. Blood flow was simulated using the lattice Boltzmann method (LBM) [12] in conjunction with the immersed boundary method (IBM) [13] to account for fluid-membrane interactions between the flow field and deformable RBCs. Guo's method was adopted to account for spatiotemporal changes in forces on the flow and to model unsteady flows accurately [14]. At the upstream and downstream boundaries, periodic boundary conditions with a pressure difference were imposed. The no-slip boundary condition was applied to the boundary nodes to represent that the flow velocity at the wall matches the velocity of the wall.
To model RBCs, the unstressed structure was given the biconcave shape proposed by Evans and Fung [15]. Each RBC membrane was modeled with 1,272 triangular elements. Membrane forces due to strain, bending, surface, and volume constraints were computed from the corresponding energy terms using the principle of virtual work [16]. Reference values were used to determine the parameters of the RBC membrane [17]. Viscosity inside the RBCs was assumed to be five times greater than the viscosity of plasma without LEH.
The surrounding fluid was modeled as both plasma without LEH (Newtonian) and plasma with LEH (non-Newtonian). For non-Newtonian fluid, we applied the following power-law model to describe the fluid viscosity, where µ is the apparent viscous coefficient, γ is the shear rate, and n is an index that distinguishes non-Newtonian and Newtonian fluids. µ 1 means the viscous coefficient at n = 1.0 (Newtonian fluid). Figure 2 plots the relationship between shear rate and the viscosity of the surrounding fluid that follows from Eq. 1. For the non-Newtonian fluid, we assumed n = 1 and µ 1 = 1.2 mPa･s. Figure 3 shows the distribution of viscosity in the microvessel. The viscosity of the Newtonian fluid was constant in the vessel. For the non-Newtonian fluid, in contrast, the high-viscosity fluid was concentrated around the vessel axis since the flow velocity was low near the axis. This high-viscosity region was approximately five times more viscous than the fluid near the vessel walls.

Single RBC
First, we calculated the flow of single RBC in Newtonian and non-Newtonian fluids. Figure 4 shows the trajectories of the center of gravity of the simulated RBCs. We considered two initial positions with y 0 set to 2.0 and 5.0 µm. The position of y 0 = 5.0 µm is the midpoint between the vessel centerline (y = 0 µm) and wall (y = 10 µm). As time passes, the y-coordinates of the RBCs changed. Shear thinning in the non-Newtonian fluid tended to promote the axial concentration of RBCs. In a non-Newtonian fluid, viscosity is higher at the center of the vessel. In this simulation, since the mean velocity was constant, the capillary number (Ca = µU m /DE s ) increases toward the vessel center, where E s is the strain modulus of the RBC membrane. A higher capillary number causes larger deformation of the RBC. The associated variation in the viscosity in the vessel encourages the RBCs to cluster around the vessel center.

Multiple RBCs
Next, we calculated the flow of multiple RBCs. Figure 5 plots the RBC flow for the three hematocrit values we simulated. Figure 6 shows the time average of the RBC distribution in the vessel. The horizontal axis gives the distance from the vessel centerline, and the vertical axis indicates the area ratio of RBC passing through the downstream cross-section. These results show that a cellfree layer formed due to the axial concentration of RBCs. This cell-free layer became thicker as the hematocrit decreased. In the case of Hct = 3.0 %, the RBC distributions in Newtonian and non-Newtonian fluids are very different. With small hematocrit values, shear thinning dominates the axial concentration of RBCs. As the hematocrit increases, mutual interference between RBCs has more of an effect, so the distribution of RBCs is less different between Newtonian and non-Newtonian fluids. When Hct = 6.1 %, the RBCs were more clustered near r = 4 µm in non-Newtonian fluid than in Newtonian fluid. To investigate this phenomenon in more detail, we defined the orientation angle, which is the angle between x-direction and the longitudinal direction of the RBC. We recorded the orientation angles

Conclusion
In this study, to clarify the flow behavior of RBCs in blood with added LEH, the effects of shear thinning on the flow in blood vessels were investigated using threedimensional microvasucular models. First, we investigated the flow behavior of single blood cell in Newtonian and non-Newtonian fluids. The simulation results indicate that shear thinning encourages RBCs to concentrate around the vessel axis. Next, we investigated the flow behavior of multiple blood cells by varying the hematocrit. These simulation results indicate that shear thinning has more of an effect at low hematocrit levels.
The orientation angle of the RBCs was finally investigated. The difference in the orientation angles in Newtonian and non-Newtonian fluids explains the variation in RBC distribution within microvessels. As the hematocrit increases further, the interactions between the RBCs have more of an effect, so the discrepancy in the distribution RBCs in Newtonian and non-Newtonian fluids nearly disappeared.