Enhanced discrete phase model for multiphase flow simulation of blood flow with high shear stress

Abstract

Shear stress is often present in the blood flow within blood-contacting devices, which is the leading cause of hemolysis. However, the simulation method for blood flow with shear stress is still not perfect, especially the multiphase flow model and experimental verification. In this regard, this study proposes an enhanced discrete phase model for multiphase flow simulation of blood flow with shear stress. This simulation is based on the discrete phase model (DPM). According to the multiphase flow characteristics of blood, a virtual mass force model and a pressure gradient influence model are added to the calculation of cell particle motion. In the experimental verification, nozzle models were designed to simulate the flow with shear stress, varying the degree of shear stress through different nozzle sizes. The microscopic flow was measured by the Particle Image Velocimetry (PIV) experimental method. The comparison of the turbulence models and the verification of the simulation accuracy were carried out based on the experimental results. The result demonstrates that the simulation effect of the SST k-ω model is better than other standard turbulence models. Accuracy analysis proves that the simulation results are accurate and can capture the movement of cell-level particles in the flow with shear stress. The results of the research are conducive to obtaining accurate and comprehensive analysis results in the equipment development phase.

Keywords

Computation fluid dynamics multiphase flow blood-contacting devices shear stress discrete phase model particle flow turbulence model particle image velocimetry experiment

Introduction

Blood-contacting devices include blood pumps, prosthetic heart valves, blood filters. The internal flow of these devices is complex, and flow with shear stress usually has the risk of causing hemolysis.^1,2,3 Hemolysis is a phenomenon in which red blood cells are damaged in flow. Excessive shear stress in the flow cause deformation and rupture of red blood cells. Hemolysis can cause red blood cells to lose the ability to transport oxygen, which seriously endangers the life of the disease.⁴

In blood-contacting devices, the shear stress in the flow is a core factor leading to hemolysis. In the flow with shear stress, the shear stress due to the velocity gradient causes the hemolysis risk of the device.^5,6,7 The hemolysis test can directly obtain the hemolysis index of the device, but this method has several defects. First, the hemolysis experiment is to measure the hemolysis index of the device, but this does not capture the connection between flow and hemolysis, so it is difficult to support the development and improvement of the device structure. Second, the results of hemolysis experiments are also affected by many inevitable factors, which often mask the resulting characteristics of the device structure.⁸ Also, cost and time are also factors to consider. Therefore, in the research and design stage, the CFD method is usually used to simulate the internal flow of the equipment. In this way, device flow rationality and hemolytic damage are estimated.^9,10 The design requirements of blood-contacting devices are much higher than those of ordinary equipment, which requires CFD simulation to have sufficient accuracy.

CFD simulation is an important method for blood flow analysis in devices. The simulation results are used to evaluate and improve devices.^11,12 For example, in the research and development of blood pumps, researchers use CFD methods to analyse and evaluate the flow field state and hemolysis level of different structures, to promote the optimization of the structure and performance of blood pumps.^13,14,15,16 However, most current studies have simplified blood into a single-phase flow liquid. Although the results of this simplified method are considered to be reliable. But the method of multiphase flow is also one of the directions of research. Blood is a typical liquid-particle multiphase flow. It is mainly composed of plasma, red blood cells, white blood cells, and platelets. The purpose of hemolysis analysis is to predict the damage of red blood cells during flow, thereby evaluating, and improving the structure of the device. Therefore, liquid-particle multiphase flow simulation is not only closer to real physical characteristics, but also more in line with the needs of research purpose. This will help to obtain more accurate and comprehensive analytical data.

At present, the multiphase flow simulation method of blood is not perfect. Multiphase flow models and experimental validation are major issues. In this study, an enhanced multiphase flow simulation method is proposed. The multiphase flow model is based on DPM, and the model is analysed and improved based on the multiphase flow characteristics of blood. In the experimental verification, we designed nozzle models to simulate the flow with shear stress, varying the degree of shear stress through different nozzle sizes. The results show that the simulation effect of the SST k-ω model is better than other standard turbulence models. Verification proves that the simulation results are excellent and can enable flow analysis of cell-level particles in the blood flow with shear stress. The research results will help to achieve more accurate and comprehensive flow analysis during the development phase of the equipment.

Materials and methods

Nozzle model

The shear stress in the flow is due to the velocity gradient. The nozzle can produce flow with shear stress similar to blood-contacting devices, which are used to verify simulation accuracy.^17,18,19,20 In the nozzle, the fluid velocity accelerates in the tapered section of the pipe. Then, the fluid is ejected from the nozzle jet position to the shear stress region. After the injection, the fluid velocity will rapidly decelerate in the shear stress region. After spraying, the fluid velocity will rapidly decrease in the shear stress region. The velocity reduction gradient will produce shear stress, which will cause hemolysis damage when it is too large. This flow with shear stress is the leading cause of hemolysis risk in blood-contacting devices.

For this study, the model structure has been modified to produce varying degrees of shear stress. Figure 1 shows the nozzle model in this study. The model is designed with three different nozzle sizes. The nozzles are 5, 6, and 7 mm in diameter, corresponding to three levels of shear stress at high, medium and low. Hereinafter referred to as 5 mm model, 6 mm model, and 7 mm model, respectively. The model has a circulating flow of 5 L/min, which produces a shear stress at the position of the dashed box in Figure 1. In the previous calculations of each model, the maximum velocity of the nozzles was about 4.3, 3.2, and 2.5 m/s, respectively. The maximum shear stress in the model flow field is approximately 400, 300, and 200 Pa, respectively. In this way, studies of various degrees of blood flow with shear stress can be achieved. In the experimental and simulated flow analysis, the main focus is on the shear stress region, centerline and radial cuts lines (Line 1, Line 2) of Figure 1. Among them, Line 1 and Line 2 are 10 and 20 mm from Original Point, respectively. The original point is the center of the nozzle jet.

Figure 1.

Nozzle model.

DPM for blood flow with shear stress

Blood is mainly composed of plasma, red blood cells, white blood cells, and platelets, which is a typical liquid-particle multiphase flow. Among them, red blood cells account for more than 95% of blood cells, and the movement and damage of red blood cells is the focus of related equipment research.^21,22 In this study, blood was considered to be composed of plasma and red blood cells, ignoring the non-Newtonian properties of the blood. Blood is a non-Newtonian fluid under normal flow conditions. Non-Newtonian fluids have a non-linear relationship between shear stress and shear strain rate, which is common in biological fluids. However, in some studies of implanted devices, it was found that when the shear rate in the flow field is large, the blood will behave like Newtonian fluid.²³ In a high shear rate flow field, the viscosity of blood will stabilize within a small range. Therefore, in the study of blood shear damage, blood is generally regarded as Newtonian fluid, and this study also uses this fluid setting.

The red blood cells are in the shape of a double concave ellipsoid in a static state, and gradually deform into an ellipsoid shape in the flow with shear stress. In this study, based on the current research progress and the limitations of computing software, we simplify the particle to a sphere with a size close to that of red blood cells. In studies that do not consider microscopic cell deformation, this simplification is commonly adopted and is considered to be reliable.²⁴ Moreover, the goal of this multiphase flow calculation is to obtain the flow of particles, because the particle size of the cell particles is small, and the density is close to the liquid phase, which makes the tracing ability of the particles very high.²⁵ In this multiphase flow, the effect of particle shape on particle flow is basically negligible.

This multiphase flow model is based on the discrete phase model (DPM). The model processes the liquid phase into a continuum by solving the Navier-Stokes equation, while tracking the discrete phase particles through the calculated flow field. The calculations of the continuous phase and the discrete phase are in the Euler reference frame and the Lagrangian reference frame, respectively. Based on the calculated targets, multiphase flow model and current hardware computing capabilities, we have the following simplifications and assumptions in particle motion calculations: All particles are simplified to be spherical with the same size, and the particles are used as the moving mass point in the calculation. Some details of the flow around the particles are neglected, such as vortex shedding, flow separation and boundary layers. local properties of the dispersed phase are predicted from spatial averaging over particle trajectory. The model predicts the trajectory of a discrete phase particle by integrating the force balance on the particle. This force balance equates the particle inertia with the forces acting on the particle, and can be written as

\frac{d u_{p}}{dt} = \frac{u - u_{p}}{τ_{r}} + \frac{g (ρ_{p} - ρ)}{ρ_{p}} + F

(1)

where F is an additional acceleration (force/unit particle mass) term. ρ and ρ_p are the densities of the continuous and discrete phases, respectively. u and u_p are the velocity of the liquid phase and the particle phase, respectively. The first item on the right is the drag force per unit particle mass, τ_r is the droplet or particle relaxation time.

The difference in density between the particle phase and the liquid phase of the blood is less than 5.0%, which is quite different from the general particle flow. When the density of the fluid is close to the particle, the influence of the virtual mass force and the pressure gradient on the particle balance is not negligible. Equation (1) incorporates additional forces F in the particle force balance that can be important under special circumstances. The virtual mass force is the force required to accelerate the fluid around the particle. This force can be written as

F_{vm} = C_{vm} \frac{ρ}{ρ_{p}} (u_{p} \nabla u - \frac{d u_{p}}{dt})

(2)

where C_vm is the virtual mass factor with a default value of 0.5. Instructions for the choice of particle shape: For micron-sized particles such as red blood cells, the flow drag is mainly affected by particle size and density, and the difference in fine shape has little effect on flow.^26,27 At the same time, in order to better analyze and demonstrate the effectiveness of simulation and models through PIV experiments, we hope to ensure that the simulation and experimental conditions are close. For the most critical tracer particles in the experiment, we mainly ensure that the particles are as close as possible to the red blood cells in terms of density and size. Therefore, considering the current calculation conditions, simulation methods, experimental techniques and research status, we simplified the red blood cells as spherical particles.

In addition, there will be additional forces due to the pressure gradient in the fluid. This force can be written as

F_{p} = \frac{ρ}{ρ_{p}} u_{p} \nabla u

(3)

Numerical simulation

Numerical simulations were performed using the commercial software ANSYS CFX 17.0 software (ANSYS, Inc., Canonsburg, PA, USA). In the simulation, the flow calculation of the particles is achieved by the improved DPM. To verify the simulation accuracy through experiments, the particle phase and liquid phase settings in the simulation are consistent with the experiments, which are air-glass microspheres and glycerin solution, respectively. Specific analysis of the particulate phase and the liquid phase will be described in detail later. The boundary conditions are the velocity inlet and pressure outlet, and the inlet flow and outlet pressure are set to 5 L/min and 0 Pa, respectively. The particle inlet is 0.01 kg/s. The calculated convergence residual is 10⁻⁶. In the calculation, the key parameters such as inlet pressure, outlet pressure, inlet flow, and outlet flow are monitored to ensure that the values of these parameters are stable when the calculation reaches convergence.

Turbulent flow is highly complex, and it is difficult to achieve direct calculations like laminar flow. The reason is that there are nearly infinite multi-scale vortex flows and strong nonlinearities in turbulence, which makes it difficult to solve the turbulence problem by theoretical experiments and simulations. In CFD simulations, the calculation of turbulent flow usually needs to be assisted by the turbulence model. The turbulence model is based on Reynolds Averaged Navier-Stokes (RANS) equations. The model averages the physical quantities of the flow field and solves the time-averaged governing equations. The flow appears as a distinct turbulent state in the injection region. For the simulation of turbulent flow, the turbulence model is a crucial factor affecting the accuracy of the results.²⁸ For flow conditions, the applicability of each turbulence model is different. There is no uniform theory or method to determine the applicability of the turbulence model. In the case of high precision requirements, it is best to choose the appropriate turbulence model through experiments.

The 3D computational grid of the nozzle model was established by ANSYS ICEM. The grid is one of the factors that have the greatest influence on the calculation results in the simulation. The density, type and scheme of the grid largely depend on the geometric structure of the computing object. In the calculation of three-dimensional flow, commonly used grid types include hexahedral grid, unstructured grid and mixed grid. Compared with the use of structured hexahedral grids, the use of unstructured grids usually requires more cells. However, the establishment of a hexahedral grid has higher technical requirements for the operator, and will increase the time cost of research. In the case of insufficient skills, the poor quality hexahedral grid will affect the accuracy of the calculation results. Therefore, unstructured grids or mixed grids are often used in many blood devices with complex geometric structures. In the simulation of three-dimensional flow, when the density and quality of the grid are sufficient, it is generally considered that the unstructured grid can meet the research requirements.^23,29

Quantity is one of the most important factors of the grid. The problem of this grid independence analysis is that when the grid is established densely, the amount of calculation will increase and the calculation cycle will be longer. Without affecting the results, waste of computing resources and time should be avoided. Secondly, as the grid is encrypted, rounding errors caused by computer floating-point operations may also increase. This calculation uses an unstructured grid, the grid type is Tetra/Mixed, and the grid creation method is Robust (Octree). The purpose of this grid independence analysis is to ensure that the number of grids will not affect the calculation results. Because the input of the simulation is the inlet flow rate, the key output of the calculation result is the inlet and outlet pressure. In the verification, it is necessary to ensure that the number of grids will not change the pressure of the inlet and outlet. This verification method has been widely adopted in related studies.²³ Table 1 shows the results of this grid independence analysis, taking the results of the 6 mm model as an example. Table 1 contains the calculation results of five groups of different grid number models. The comparative analysis shows that when the number of grids reaches 4 million, continuing to increase the number of grids will not affect the key calculation results. Through the verification of the independence of the number of grids, the number of grids is determined to be about 4.5 million, and the quality of the grids is all above 0.6.

Table 1.

Verification of the independence of the number of grids.

Number of grids	About 1 million	About 2 million	About 4 million	About 6 million	About 8 million
Global element seed size	0.87	0.65	0.52	0.45	0.41
Lowest grid element quality	About 0.30	About 0.55	About 0.60	About 0.60	About 0.60
Convergence residual	10⁻⁴	10⁻⁶	10⁻⁶	10⁻⁶	10⁻⁶
Pressure difference between nlet and outlet (Pa)	About 3900	About 4200	About 4335	About 4335	About 4335

In this study, the effects of each turbulence model are compared and evaluated based on the results of PIV experiments. The choice of turbulence models for comparison is mainly based on literature, books and our previous experience. The turbulence models for this comparison include k-ε, k-ω, and SST k-ω. Among them, the k-ε model is one of the most widely used turbulence models in CFD, and its coefficients are given by empirical formulas, which have a good calculation effect on turbulent flow. The advantage of the k-ω model lies in the calculation of the wall boundary layer and free shear flow. And the SST k-ω model is a combination of k-ε and k-ω models. The above turbulence models are the most commonly used in related simulations.^29,30

PIV experiment

The PIV experiment was used to measure the internal flow of the model, and the accuracy of the simulation was verified by experimental results. PIV is a transient, multi-point, contactless flow velocity measurement method and is considered to be the most reliable flow measurement method.^31,32 This study used the PIV experimental system produced by TSI. The image collector of the PIV system is a Zyla 5.5 high-speed camera. The maximum acquisition frequency of this camera is 100 Hz. And the camera’s resolution is 2560 × 2160 pixels, which is approximately equivalent to 5 MP. The laser has a maximum rate of 30 Hz and maximum energy of 100 mJ, which is a two-pulse laser. The material of the experimental model is made of transparent plexiglass, and a rectangular water tank that reduces the effect of light refraction is placed outside the model. Before the experimental measurement, an appropriate amount of tracer particles are injected, and the focus and size calibration of the field of view is completed. In the measurement, the laser illuminates the flow field, and the tracer particles reflect the laser light. The high-speed camera takes two frames of shots at 100 μs intervals and calculates the velocity by particle displacement in post-processing.

v_{x} = \frac{x (t) - x (t + Δ t)}{Δ t}

(4)

v_{y} = \frac{y (t) - y (t + Δ t)}{Δ t}

(5)

Where v_x and v_y are the instantaneous velocity calculation results in the x and y directions of the particle, respectively. Δt is the shooting interval of two frames of images, x(t)−x(t + Δt) and y(t)−y(t + Δt) are the displacement amounts of the two directions in the Δt time, respectively. Because the Δt of the PIV experiment is usually around 10^-4 s, the calculation result can be approximated as an instantaneous velocity. In the post-processing, the image is divided into regions, and 5–10 particles are guaranteed in each region, and the average velocity value of these particles is taken as the velocity result of the region.

This PIV experiment uses the 2D2C method, and the experimental result is the flow in the 2D plane. However, the simulation of this study is a 3D flow calculation, and the PIV experiment will lack the flow velocity perpendicular to the plane. For this study, this difference will not affect the comparative analysis of experiment and simulation, because the flow velocity perpendicular to the measurement plane is very small. For example, in the simulation results of the 6 mm model, the velocity range of the analysis plane is about 0.2 to 3.0 m/s. The maximum value of the partial velocity perpendicular to the plane is about 0.03 m/s, so we ignore the velocity in this direction, and the 2D plane velocity results are used in the comparative analysis of experiments and simulations.

This experiment is different from the ordinary PIV experiment. Typically, the PIV experiment uses the best tracer of flow following. In the results, the particle velocity is considered to be equal to the liquid flow velocity. However, this study is directed to liquid-particle blood multiphase flow. Therefore, the selection of the tracer particles and liquids does not focus on flow followability, but rather makes the experimental liquid-particle multiphase flow properties closer to the blood.

In the blood, the diameter and density of red blood cells are about 8 μm and 1.1 kg/m³, respectively. In this regard, we use air-glass microspheres produced by TSI, which have an average particle size and density of 10 μm and 1.1 kg/m³, respectively. The fluid used in this experiment is a water/glycerin mixture (62/38) solution by percent mass. The specific gravity and dynamic viscosity of the glycerin solution are about 1.08 and 3.50 MPa·s, respectively. The dynamic viscosity of blood is related to the shear rate, and is usually divided into high shear rate, medium shear rate, and low shear rate in medicine. Table 2 shows the range of dynamic viscosity and shear rate. However, the shear rate in blood equipment is higher than usual and varies with flow. Therefore, a glycerin solution with a dynamic viscosity of 3.50 MPa·s is usually used in the research. This method is generally considered reliable in relevant studies.

Table 2.

Blood dynamic viscosity and flow shear rate.

Flow condition	Low shear rate	Medium shear rate	High shear rate	Blood device
Shear rate (1/s)	1–20	50–60	100–200	>200
Dynamic viscosity (mPa·s)	8.2–9.6	5.5–6.4	4.4–4.9	Assume about 3.5

Figure 2(a) shows the arrangement of the PIV experimental equipment. The laser is irradiated vertically downward from the top of the model, forming a slice at the axial position of the model, and the high-speed camera shoots perpendicularly to the laser light. Figure 2(b) is an original particle image.

Figure 2.

PIV experimental device and particle image: (a) PIV experimental device and (b) PIV original particle image.

Results

Velocity field analysis

Figure 3 shows the results of the PIV experiment and the simulation results of different turbulence models. The results contain three sets of models. In the results, the flow is expressed as an injection state with a high-velocity gradient. In the shear stress region, the flow velocity is higher in the middle of the pipe and gradually decreases with the flow direction. Also, the velocity near the wall of the pipe is low, and there is a weak backflow. The smaller nozzle jet size has a higher velocity, and the highest velocity is at the original point. Compared with experimental and simulation results, the maximum velocity of the simulation results of the 5 mm model, 6 mm model, and 7 mm model are about 4.10, 3.20, and 2.50 m/s, respectively, and the difference between the turbulence model results is less than 0.05 m/s. In the experimental results, the maximum velocity is 4.09, 3.22, and 2.45 m/s, respectively. The difference between the experiment and the simulation is less than 1.0%.

Figure 3.

PIV experiment and simulated velocity results (The results are the resultant velocity of vertical velocity and axial velocity, and the corresponding velocity plane is shown in Figure 2).

The difference between the results of each model is mainly reflected in the flow state and velocity distribution. The flow state of the k-ε model results is close to the experimental results, but the distance in the high-velocity region is slightly farther than the experimental results. The k-ω model is quite different from the experimental results, and its jet flow is more slender. Obviously, the results of the SST k-ω model are closer to the experimental results in both the flow state and the velocity distribution. Only the distribution of velocity will be more symmetrical and stable than the experimental results. In the comparison of results, the SST k-ω model shows better simulation accuracy. In the subsequent content, a more comprehensive and detailed analysis of the accuracy of the SST k-ω model results will be made.

Figure 4 shows the 3D particle motion trajectory obtained by this simulation. Since the motion characteristics of the three sets of model results are close, the analysis uses the results of the 6 mm model as an example. In the results, the velocity of the particles flows in the range of approximately −0.5 to 3.0 m/s, which is transiently above 3.0 m/s near the nozzle original point. The flow of the model inlet and outlet is stable with an average velocity of around 1.0 m/s. The high-velocity region of the particle flow is located in the central area behind the original point and is obscured by the particle trajectory near the wall of the pipe in the figure, which can be clearly observed in Figure 4. Besides, there are low velocity and recirculation regions near the wall of the pipe. The particle trajectory shows the motion state of the particles in the flow field, and the simulation results have no obvious clutter and errors.

Figure 4.

Three-dimensional particle motion trajectory.

Velocity along the centerline

Figure 5 shows the velocity curve of the centerline, which is used to evaluate the accuracy of the simulation in the main direction. In the result, the data starts from the original point and takes 21 data points along the axis. The distance between each data point is 0.5 mm. (a), (b), and (c) of Figure 5 are the results of the 5 mm model, the 6 mm model, and the 7 mm model, respectively.

Figure 5.

Experimental and simulation results of the centerline: (a) 5 mm model, (b) 6 mm model, and (c) 7 mm model.

In the results, the trends of the velocity curves of the three groups of models are basically the same. There is a brief rise in velocity from the starting point to the second data point, and the velocity will remain high in the range of around 2 mm. Then, the velocity will decline with a basically stable trend. For each model, different nozzle diameters only change the velocity range of the model and do not affect the trend of the velocity curve. Further, in the 0 to 20 mm portion of the centerline, the matching between the experimental and simulated results is perfect, and the velocity difference is substantially below 1.0%. In the 30 to 40 mm section of the centerline, the difference between the experiment and the simulation increased slightly, but the maximum difference did not exceed 5.0%. Overall, the simulation has good accuracy on the centerline.

In order to further improve the reliability of the simulation results, the comparative analysis and selection of the turbulence model which is the most critical for simulation are carried out. Figure 6 is the comparison result. The comparative analysis is based on the axial velocity. The axis velocity can most directly reflect the accuracy of flow changes. The ordinate value is the velocity difference between simulation and experiment.

Figure 6.

The velocity difference between the simulation and experiment of each turbulence model on the axis: (a) 5 mm model, (b) 6 mm model, and (c) 7 mm model.

In the results, there are some differences between the results. The simulation results of SST k-ω model have the best match with the experiment, and the velocity difference is within 0.1 m/s. Especially in the range of 0 to 20 mm, the velocity difference is within 0.02 m/s. The other two models have significantly larger velocity differences. In addition, the velocity difference of each model gradually increases in the range of 20 to 40 mm. However, the velocity difference of the SST k-ω model has the smallest increase, and the maximum value of the velocity difference is only about 40.0% and 16.7% of the other two models, respectively.

Velocity along with radial cuts

Figure 7 shows the velocity along with the radial cuts, which is used to evaluate the accuracy of the simulation in the radial direction. In the result, the data starts from the axis line. Along the Line1 and Line2, take 8 data points in each direction, and the distance between each data point is 0.25 mm. (a), (b), and (c) of Figure 6 are the results of the 5 mm model, the 6 mm model, and the 7 mm model, respectively.

Figure 7.

Experimental and simulation results of the radial cut line: (a) PIV and CFD velocity curve of 5 mm nozzle model, (b) PIV and CFD velocity curve of 6 mm nozzle model, and (c) PIV and CFD velocity curve of 7 mm nozzle model.

In the results, the three groups of models have similar velocity characteristics, and the simulation is very close to the experimental results. The velocity has a maximum at the axis position, and the velocity value decreases rapidly toward both sides. The velocity reduction trend on both sides of the axis is almost symmetrical. The velocity is close to zero near the wall, and there is a backflow near the wall. The degree of reflow increases as the nozzle jet diameter decreases.

Compare the experimental and simulation results of each model. The simulation and experimental results are well matched, and the corresponding velocity curves are almost coincident. In comparison, Line 1 results are better matched than Line 2, and there is no significant difference between the different models. The results in two locations are worthy of attention. First, at the axial position, the Line1 simulation results are slightly higher than the experimental results, and the difference is about 5.0%. In Line 1, there is no such velocity difference. Second, there is a difference between the simulation and the experimental results in the recirculation region, which is present in both Line 1 and Line 2. However, the difference in velocity is less than 0.05 m/s. In general, the radial velocity also has good simulation accuracy, but there may be subtle simulation errors in the recirculation region.

Shear stress analysis

Figure 8 shows the shear stress results of the PIV experiment and simulation. The results include 3 sets of models (5 mm model, 6 mm model, and 7 mm model). The result of Figure 8 is the 2D plane shear stress distribution, which is calculated by the velocity gradient. Combined with the results of Figure 3, the shear stress is mainly distributed on both sides of the high-velocity region, and gradually decreases with the flow. The high shear stress region appears on both sides of the inlet, and the model with higher flow velocity also has higher shear stress. The stress distribution of the simulation results is more regular than the experimental results, but the overall distribution characteristics are consistent. Moreover, the experimental and simulated stress values are also in the same range. Comparing the results of different turbulence models, the difference in stress is similar to the difference in flow velocity in Figure 3. In Figure 3, the jet flow of the k-ω model results is longer, and the high shear stress region in Figure 8 is also longer. In theory, because the shear stress is based on the velocity gradient results, Figures 3 and 8 will reflect the associated characteristics. Combining the results of Figures 3 and 8, SST k-ω shows better calculation results.

Figure 8.

Shear stress results of PIV and simulation.

Discussion

In this study, the liquid-particle multiphase flow simulation method of blood flow with shear stress was analysed and improved, and the simulation accuracy was verified by PIV experiment. The model is based on DPM. In the composition of multiphase flow, since most of the components of blood cells are red blood cells, and the main purpose of blood flow with shear stress analysis is to acquire the movement and damage of red blood cells. Therefore, the blood is simplified by a two-phase fluid composed of plasma and red blood cells. Second, since the density of the blood liquid phase and the particle phase are close, the effect of the virtual mass force and the pressure gradient on the particle balance is not negligible. Therefore, the virtual mass force model and the pressure gradient model are added to the DPM model, which makes the calculation of the model more comprehensive and accurate.

Experimental verification is necessary to determine the effectiveness of the simulation method. The simulated flow data is the basis and source of data for calculating shear stress and hemolysis estimates. Therefore, the demonstration of simulation accuracy focuses on the accuracy of simulated flow data. For the characteristics of blood multiphase flow, the PIV experimental method was adopted. The nozzle model was improved and designed. Through different nozzle diameters, the model can produce varying degrees of shear stress. To better verify the multiphase flow simulation results, the air-glass beads and glycerol solution were used in experiments, and their physical properties were close to those of red blood cells and plasma. The flow results of each model were obtained through experiments.

A suitable turbulence model is a prerequisite for accurate simulation. Under the requirement of high precision, the experiment is the best method to analyse and select the turbulence model. In this paper, the turbulence models are compared and selected based on the results of the PIV experiment. In the results of each turbulence model, the flow is expressed as the injection state, which is consistent with the experimental results. The velocity range is basically consistent with the experimental results, and the maximum velocity difference is only 1.0%. In these respects, each model can achieve accurate simulations with little difference from each other. However, the distance in the high-velocity region of the k-ε model results is slightly farther than the experimental results. The high-velocity region of the k-ω model is significantly narrower and longer than the experimental results. However, the results of the SST k-ω model do not have these error problems, and its flow state and velocity distribution are well matched to the experimental results. In addition, the simulated distribution is more symmetrical and stable than the experimental results, which is due to the principles of experimentation and simulation. Simulation is an idealized calculation, and experiments have unavoidable instability. In general, the SST k-ω model is better than the other two commonly used turbulence models.

Further, the study conducted a more comprehensive and detailed analysis of the simulation accuracy. In the results of the axial velocity curve, the trends of the experimental and simulated velocity curves are basically the same. On the curve, the velocity rises briefly first, then the velocity will maintain a high-velocity region, and finally the velocity will drop steadily. Comparing the experimental and simulated velocity values, the velocity difference is substantially below 1.0% in the centerline 0 to 20 mm. In the 30 to 40 mm section, the difference between the experiment and the simulation increased slightly, but the maximum difference did not exceed 5.0%. The results reflect the accuracy of the simulated velocity changes in the flow direction. Also, there is better accuracy at a location close to the original point.As a result of the radial cut line velocity, the velocity at the axis position reaches the highest value, and the velocity values on both sides appear as asymmetrical and stable falling state. The velocity is close to zero on the wall, and there is a backflow near the wall. Comparing the results of simulations and experiments, their velocity curves are basically coincident.

However, there will be a slight error in the two locations. The first is the axial position. The velocity value of the simulation result will be slightly larger than the experimental result, and the difference will gradually increase with the flow direction. Both the results of Figures 5 and 6 reflect this feature. The second is the low-velocity recirculation region near the wall. Reflow is one of the flows that simulate difficult to calculate accurately. This simulation accurately reflects the position of the reflow. For the velocity value of the recirculation region, the difference between the simulation and the experiment is less than 0.05 m/s, and the difference has reached a good accuracy. In summary, the comparison between experiment and simulation reflects the excellent accuracy of the proposed simulation method.

Finally, the PIV experiment and simulated shear stress distribution are analyzed. The result of Figure 8 is the shear stress distribution in the 2D plane calculated by the velocity gradient. Because the shear stress result is calculated by the flow velocity gradient, the flow velocity result will show a strong correlation with the shear stress result, and the results of Figures 3 and 8 reflect this feature. For example, the high shear stress region of Figure 8 is located on both sides of the high flow velocity region of Figure 3, and there is almost no stress in the high flow velocity region, because the high shear stress is caused by the change of flow velocity. Overall, Figures 3 and 8 reflect that the SST k-ω model has better calculation results, and the model results have sufficient accuracy for the calculation of flow velocity and stress.

There were several limitations to this study. First, the evaluation of the simulation accuracy is entirely dependent on the PIV experiment. Although the PIV experiment is considered to be the best flow measurement method, it may not be enough. We will improve the level of experimentation and equipment and more comprehensively evaluate the accuracy of results in future work. Secondly, this study used a nozzle to simulate blood flow with shear stress, which was not applied to actual medical equipment. Although the reliability of nozzle is generally recognized. However, in the next study, the method will be applied to the blood pump, thereby further verifying the application value of the research results. In addition, due to the limitation of calculation and experiment, the size and density of the particles were matched with the red blood cells in the study, but the particle concentration was not completely consistent with the blood. The problem of particle concentration is a limitation in this study, but we think this study still has significance and progress in some aspects. The purpose of this study is to establish a multiphase flow model that can predict the flow of micro-particles. The results are more concerned about the flow process and velocity state of the particles. Judging from the results and experiments, the goal of this study has basically been achieved. In the analysis of blood equipment, most researches treat blood as a simple liquid, and the flow of liquid is directly equivalent to the flow of red blood cells. The multiphase flow calculation in this study is an independent flow calculation of red blood cells. In comparison, this study has made progress. Because blood is a complex biological fluid, some simplifications and limitations are included in this calculation. Limited by existing software and technology, some problems are difficult to solve in a short time. In this regard, we will try to solve these problems through theoretical and software breakthroughs. In addition, Khoo et al. recently found that normal stresses are important in a centrifugal blood pump.³³ It is necessary to consider the combined effects of shear force and normal stresses in future studies. In the subsequent research, we will continue to improve the calculation method and experimental method, and gradually make the calculation closer to the real blood.

Conclusions

This research is aimed at multiphase flow simulation of blood flow with high shear stress, and the application object is blood-contacting devices such as blood pumps. The main purpose of the research is to establish and demonstrate a simulation method that can realize the analysis of micro-particle movement and stress in shear flow. Through the establishment of multiphase flow model and PIV micro-particle experiment, this research proposes an enhanced discrete phase model for multiphase flow simulation of blood flow with high shear stress, and completed the analysis, improvement and demonstration of the multiphase flow model. The main research contents and results are summarized as follows:

Research has established a simulation multiphase flow model for blood flow with high shear stress, which realizes the simulation analysis of cell-level micro-particle flow and stress. The established multiphase flow model is based on DPM. According to the multiphase flow characteristics of blood, dynamic calculations of the influence of virtual mass force and pressure gradient are added to the micro-particles. (2) For the turbulence model that is critical to the accuracy of the simulation results, a comparative analysis and selection of the applicability of the high shear stress model are carried out. The comparative analysis shows that the SST k-ω turbulence model shows better simulation applicability to blood flow with high shear stress. The SST k-ω turbulence model can avoid the problem of over-prediction in the high-velocity zone in the simulation results of other turbulence models. (3) Aiming at the flow field of blood flow with high shear stress, nozzle experimental models are designed to simulate the flow field state required for this research. The designed nozzle experimental model can produce different degrees of shear flow through different nozzle sizes, which can realize a comprehensive research and analysis of different degrees of shear stress. (4) Aiming at the experimental difficulties of cell-level microparticles, the micro PIV technology is adopted, and the flow of red blood cells are simulated by micron-level tracer particles. Based on the experimental results, the established multi-phase flow model was analyzed, improved, and accuracy demonstrated. PIV experimental results show that the established simulation has fine simulation accuracy and research reliability.

In summary, the DPM multiphase flow simulation established in this study can accurately realize the flow and stress analysis of cell-level microparticles in flow with high shear stress. The results of this research provide a multiphase flow simulation model and method that is more in line with the characteristics of the blood flow field and calculation requirements, which is conducive to obtaining more accurate and comprehensive analysis data and results in the equipment development stage.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by National Natural Science Foundation Project of China, grant number 51475477 and 31670999, and Postgraduate Independent Research and Innovative Project of Central South University, grant number 2018zzts021.

ORCID iDs

Zheqin Yu

Jianping Tan

Shuai Wang

Author biographies

Zheqin Yu is a lecturer at Changsha University of Science and Technology, engaged in the research of fluid machinery and has published 10 research papers on blood pumps.

Jianping Tan is a professor at Central South University. He is currently the principal of the School of Light Alloy Research at Central South University, engaged in the design and drive research of micro fluid machinery.

Shuai Wang is a lecturer at Henan University of Science and Technology, engaged in the research of blood pump design and blood injury mechanism.

References

Meyer

Malehsa

Bara

, et al. Implantation of rotary blood pumps into 115 patients: a single-centre experience. Eur J Cardiothorac Surg 2013; 43(6): 1233–1236.

Thamsen

Blümel

Schaller

, et al. Numerical analysis of blood damage potential of the HeartMate II and HeartWare HVAD rotary blood pumps. Artif Organs 2015; 39(8): 651–659.

Dasi

Simon

Sucosky

, et al. Fluid mechanics of artificial heart valves. Clin Exp Pharmacol Physiol 2009; 36(2): 225–237.

Rother

Bell

Hillmen

, et al. The clinical sequelae of intravascular hemolysis and extracellular plasma hemoglobin: a novel mechanism of human disease. JAMA 2005; 293(13): 1653–1662.

Leverett

Hellums

Alfrey

, et al. Red blood cell damage by shear stress. Biophys J 1972; 12(3): 257–273.

Paul

Apel

Klaus

, et al. Shear stress related blood damage in laminar couette flow. Artif Organs 2003; 27(6): 517–529.

Zhang

Taskin

Fang

, et al. Study of flow-induced hemolysis using novel coquette-type blood-shearing devices. Artif Organs 2011; 35(12): 1180–1186.

Herbertson

Olia

Daly

, et al. Multilaboratory study of flow-induced hemolysis using the FDA benchmark nozzle model. Artif Organs 2015; 39(3): 237–248.

Janiga

Thévenin

Computational fluid dynamics-based design optimization method for Archimedes screw blood pumps. Artif Organs 2016; 40(4): 341–352.

10.

Mozafari

Rezaienia

Paul

, et al. The effect of geometry on the efficiency and hemolysis of centrifugal implantable blood pumps. ASAIO J 2017; 63(1): 53–59.

11.

Apel

Paul

Klaus

, et al. Assessment of hemolysis related quantities in a microaxial blood pump by computational fluid dynamics. Artif Organs 2001; 25(5): 341–347.

12.

Gouskov

Lomakin

Banin

, et al. Assessment of hemolysis in a ventricular assist axial flow blood pump. Biomed Eng 2016; 50(4): 233–236.

13.

Zhang

Hydrodynamic and hemolysis analysis on distance and clearance between impeller and diffuser of axial blood pump. J Mech Med Biol 2016; 16(2): 1650014.

14.

Noor

Parker

, et al. Investigation of the characteristics of HeartWare HVAD and Thoratec HeartMate II under steady and pulsatile flow conditions. Artif Organs 2016; 40(6): 549–560.

15.

Zhou

Wang

, et al. Numerical analysis of an axial blood pump with different impeller blade heights. J Mech Med Biol 2012; 12(3): 1250045.

16.

Maruyama

Nishida

Tsutsui

, et al. The hemolytic characteristics of monopivot magnetic suspension blood pumps with washout holes. Artif Organs 2005; 29(4): 345–348.

17.

Hariharan

Giarra

Reddy

, et al. Multilaboratory particle image velocimetry analysis of the FDA benchmark nozzle model to support validation of computational fluid dynamics simulations. J Biomech Eng 2011; 133(4): 041002.

18.

Trias

Arbona

Massó

, et al. FDA’s nozzle numerical simulation challenge: non-Newtonian fluid effects and blood damage. PLoS One 2014; 9(3): e92638.

19.

Stewart

Paterson

Burgreen

, et al. Assessment of CFD performance in simulations of an idealized medical device: results of FDA’s first computational interlaboratory study. Cardiovasc Eng Technol 2012; 3(2): 139–160.

20.

Malinauskas

Hariharan

Day

, et al. FDA benchmark medical device flow models for CFD validation. ASAIO J 2017; 63(2): 150–160.

21.

Fedosov

Peltomäki

Gompper

Deformation and dynamics of red blood cells in flow through cylindrical microchannels. Soft Matter 2014; 10(24): 4258–4267.

22.

Quinn

Pivkin

Wong

, et al. Combined simulation and experimental study of large deformation of red blood cells in microfluidic systems. Ann Biomed Eng 2011; 39(3): 1041–1050.

23.

Fraser

Taskin

Griffith

, et al. The use of computational fluid dynamics in the development of ventricular assist devices. Med Eng Phys 2011; 33(3): 263–280.

24.

Cordasco

Yazdani

Bagchi

Comparison of erythrocyte dynamics in shear flow under different stress-free configurations. Phys Fluids 2014; 26(4): 041902.

25.

Yuan

Analysis on factors influencing following features of tracer particles in centrifugal pumps. J Mech Eng 2012; 48(20): 174–181.

26.

Hadad

Liberzon

Bernhaim

, et al. Characteristics of seeding particles for PIV/PTV analysis. In: 64th annual meeting of the APS division of fluid dynamics, 20–22 November 2011, vol. 64, abstracts id. H26.005. American Physical Society.

27.

Liu

Bai

, et al. Characteristics of tracer particles used in PIV experiment of pipe flow. J Liaoning Technol Univ 2016; 35: 1057.

28.

Mahdavi

Sharifpur

Meyer

JP.

Simulation study of convective and hydrodynamic turbulent nanofluids by turbulence models. Int J Therm Sci 2016; 110: 36–51.

29.

ANSYS Inc. ANSYS Workbench Release 16.0 Documentation. ANSYS Inc., 2018.

30.

Drešar

Rutten

MCM

Gregorič

, et al. A numerical simulation of HeartAssist5 blood pump using an advanced turbulence model. ASAIO J 2018; 64(5): 673–679.

31.

Zhi

Shi-He

Yang-Zhu

, et al. Experimental study on of dynamics of particles in the flow filed with intensive gradients. Acta Phys Sin 2014; 63(18): 188301.

32.

Slama

RBH

Gilles

Chiekh

, et al. PIV for the characterization of focused field induced acoustic streaming: seeding particle choice evaluation. Ultrasonics 2017; 76: 217–226.

33.

Khoo

Cookson

Gill

, et al. Normal fluid stresses are prevalent in rotary ventricular assist devices. Int J Artif Organs 2019; 41: 738–751.