Numerical Simulation on Solid–Liquid Multiphase Flow Including Complex-Shaped Objects with Collision and Adhesion Effects Using Immersed Boundary Method

Thrombosis has become more serious sickness since the number of people who die from the thrombosis is the second largest after a cancer. The cause of thrombosis is blockage of blood flow due to a thrombus. For instance, cerebral infarction and cardiac infarction are typical diseases. An accurate prediction of these diseases by advanced medical equipment is necessary to realize a society that people can lead a long healthy life. In these days, the balloon treatment is given as the most general treatment of the thrombosis due to the arteriosclerosis. However, this treatment is reconsidered due to the risk to be the restenosis. In contrast, a stent treatment can rapidly decrease the risk of restenosis by detaining the stent at the occluded part (Fig. 1.1). There are two kinds of stents, bare metal stent (BMS) and drug eluting stent (DES) [2]. The BMS is a stent made only metals and the DES is a stent that liquates out the drags to suppress a cell proliferation. However, in recent years, a stent thrombosis developed by detaining a stent has been considered as a problem. The DES’s risk that restenosis happens is smaller than the BMS’s one. Actually, the restenosis can be decreased by using the DES. Therefore, the DES has possibility to become a late stent thrombosis which the BMS has few. The late stent thrombosis developed from a delay of a reendothelialization at the stent is a sickness observed over 30 days later after the detainment of a stent. The stent thrombosis can be prevented by the supplement of the antiplatelet agents. The development of a stent to prevent the stent thrombosis without drags is required because the antiplatelet agent has side effect according to the patient’s constitution.

For this reason, it is important to examine the condition and understand how to form the thrombus. In these days, a numerical simulation has been used to express the blood flow because it is difficult to reproduce the bloods flow using experiments. There are few researches who research on the thrombogenesis. The multiphase flow analysis is essential to consider the interaction between blood cells and blood plasma. In the previous study, Tomita et al. [3] developed the adhesion model and examined the flow including virtual platelets and an obstacle as the medical stent using one-way coupling method. As a result, they suggested that the platelet tends to deposit behind the medical stent. However, this study cannot capture the interference between the blood cells and the blood plasma. Sugiyama et al. [4] carried out the blood flow analysis by their novel numerical method for massively parallel computing of fluid-structure interaction problems. For this study, it was confirmed that the number of adhesive platelets increases with increasing of the number density of the platelet and decreasing the constriction size. Their study also treats a molecular dynamics calculation to solve an adhesion phenomenon in detail. Therefore, long time calculation by large scale computer is necessary. Accordingly, their results are limited in the initial stage of the thrombogenesis. Then, in this study, we developed a flow solver to solve an incompressible solid–liquid multiphase flow. Furthermore, a multiphase channel flow including a lot of objects is simulated to investigate how to form the thrombus.

In this study, the flow in the blood vessel is solved with Cartesian mesh and immersed boundary method. The blood cells are represented by using immersed boundary method based on level set method. The phenomenon of the adhesion between the blood cells and the vessel walls obtained in the present method is compared with the realistic adhesion model proposed by Tomita et al. [3]. Moreover, the phenomenon of the collision between the blood cells or a blood cell and a vessel wall is simulated by the developed algorithm for complex-shaped objects, for instance, red blood cells. As a result, collision and adhesion of the platelets or red blood cells are successively captured using the proposed algorithm.

In this study, the blood is treated as an incompressible Newtonian fluid. The Knudsen number based on the size of red blood cells is calculated in eqn (2.1).

The Kn, k_B, T, P and d_p mean the Knudsen number, the Boltzmann’s constant, a human body temperature, a diameter of a water molecule, total pressure in a blood vessel and a diameter of a platelet. These parameters were set in eqn (2.2).

From eqns (2.1) and (2.2), the Knudsen number is calculated in eqn (2.3).

From eqn (2.3), the blood flow is not complete continuation region. However, in this study, the blood flow is treated as a completely continuum flow as well as some present studies. The governing equations of the fluid are two-dimensional incompressible Navier–Stokes equations (eqn (2.4)) and the continuity equation (eqn (2.5)). Furthermore, the equation of motion (eqn (2.6)) and the angular momentum equation (eqn (2.6)) are coupled to solve the object motion in this study. These equations are non-dimensionalized by the freestream quantities and the representative length. Using the Einstein summation convention, non-dimensionalized governing equations are:

where i=1, 2, j=1, 2 k=3. The variables x, t, u, Re, p, L, I, ω, u_m, F and m are Cartesian coordinates, time, fluid velocity, Reynolds number based on the freestream quantities and characteristic length, pressure, angular momentum, moment of inertia, angular velocity, object velocity, force affected to the objects, non-dimensionalized time and mass of object. Also, i and j is represented according to a sum terms.

The governing equations of the fluid are discretized by the finite difference method and solved on the equally spaced Cartesian mesh. Fractional step method is applied for the coupling of the velocity and the pressure field [5]. The convection terms are evaluated by the hybrid scheme of the first-order upwind and second order central difference scheme. Here, the pres- sure distribution is determined by solving the Poisson equation using successive over-relaxation method.

In this study, object boundary is represented by immersed boundary method based on level set method and ghost cell method. This method was set to refer to the previous study by Mittal et al. [6]. The immersed boundary method gives two advantages; saving computational cost to make the body fitted grid along with the object surface and recreating the computational grid for handling, moving or deforming boundary. The level set function is used as a signed minimum distance from the surface of objects in the present method. The boundary conditions at the object surfaces are given by ghost cells. Whole cells are classified into three types, fluid cells, object cells and ghost cells (Fig. 2.1). Ghost cells exist between fluid cells and object cells, and work as guard cells with giving the boundary condition. The quantities of the ghost cells are determined by means of the quantities at the image point. Here, the image point means the edge of probe extended from the ghost cells to the fluid cells which is normal to the surface. Also, in this study, the probe length between the surface and the image point is set to be 1.45Äx. The x is the width of grids. It is based on the length of the longest edge of the cells in the two-dimensional case. The quantities at the image point are determined by the bilinear interpolation by using the quantities of the cells vicinity of the image point. The velocity at the ghost cells is determined to satisfy the non-slip condition in the object surfaces. Also, the pressure at the ghost cells is determined by Neumann condition based on the value of the image point.

We proposed new adhesion algorithm to simulate the collision and interaction of the plate- let-platelet and platelet-vessel wall. This adhesive strength determined by using Kelvin Voigt model is enforced on the immersed boundary. This algorithm is based on the adhesion model proposed by Tomita et al. [3]. The amoeba-formed projection grows from the surface of activated platelets which form sphere whereas non-activated platelets form disc-like shape. The activated platelet tends to adhere to the vessel wall due to the projection. The equation of motion for the spring-mass-damper system is given in eqn (2.8).

Here, k, c and F_a are spring constant, damping coefficient and adhesion force. The spring constant and the damping coefficient are determined according to the activation condition of the platelet. These conditions are shown below:

Here, l, v and S are grade of activation, the object velocity and shear rate, A–C and α–δ are constant. This constant is determined to express the actual phenomena of adhesion. As shown in eqn (2.9), the spring and damper coefficients become large with increasing the grade of activation. The system always becomes the state of overdamping. In the model proposed by Tomita et al. [3], the spring constant at the l = 5 is infinite value; however, in this study, large value is used to prevent the bankruptcy of calculation.

We propose a collision algorithm for complex-shaped object to analyze the red blood cells. In the past, a collision of cylinders is detected based on the positional relationship of center of gravity of the objects. Meanwhile, in the current algorithm, the collision of objects is deter- mined from the relationship of the contact of the object surfaces. When an image point of a ghost cell belonging into an object is located in the ghost cell of other object, it is detected as a collision (Fig. 2.2). The red arrows mean that the image point is in the ghost cell. The collision force which is calculated by the collision angle, the object mass, and the object velocity is assigned to the cell at the colliding surface. Then, the momentums to the normal and tangential direction are exchanged and preserved, respectively.

As a variation, calculation of channel flow was conducted and compared with the analytical solution. The conceptual diagram is shown in Fig. 3.1. The D and L mean the channel height and length. The Reynolds numbers based on the channel height and a freestream were set to be 10 and 20. The number of grid point was set to be 800 × 400. Figure 3.2 shows the velocity distribution obtained by our code and analytical solution at x = L/2. Here, U_inf means the maximum velocity. The present velocity distribution is in good agreement with the analytical solution.

The flow around a circular cylinder was solved as next verification. The Reynolds numbers based on the freestream quantities and the diameter of the cylinder were set to be between 20 and 300. The grid width was set to be 0.1, 0.05, 0.025 and 0.0125. In other words, dc_c is 10, 20, 40 and 80. The dc_c means division number of the diameter of cylinders. The drag coefficient in each Reynolds number is shown in Fig. 3.3. The previous results by Ishikawa et al. [7] and Imamura et al. [8] were obtained by Lattice Boltzmann method with IBM. The other result by Rajani et al. [9] was obtained by implicit pressure-based finite volume method. The results by Panton et al. [10] and Zhang et al. [11] were measured by experiments. As shown in Fig. 3.3, the present result shows good agreement with the previous studies under the low Reynolds number. However, difference between present result and previous study becomes large in low resolution cases. It is caused from the lack of the mesh resolution to capture the boundary layer.

Figure 3.4 shows the influence of the grid resolution on each component of the drag coefficient at the Reynolds number of 300. It is normalized by the value of the boundary-fitted coordinate (BFC) grid shown by Takahashi et al. [12]. This result is obtained by means of IBM at the Mach number of 0.3. As shown in the Fig. 3.4, total drag coefficients, pressure drag coefficients, and viscous drag coefficient show good agreement in the same grid resolution. Also, the difference between the results of present and BFC decreases as the grid resolution increases.

The collision of complex-shaped objects was investigated with comparing the energy before and after the collision. In this case, the mass, the moment of inertia and the diameter of each object are same. The left object (Obj1) enforcing the initial translational velocity moves toward the right object (Obj2) as stationary one in this simulation (Fig. 3.5). In the flow simulation, the fluid dynamic forces are NOT considered to investigate the accuracy of the current collision algorithm. The change of the energies of each objects are summa- rized in Table 3.2. The energy conservation can be archived as shown in Table 3.2. Furthermore, the rotational motion is confirmed due to the law of action and reaction in Fig. 3.5(c).

The channel flow including moving circular cylinders was investigated for the comparison with the previous study by Kamada et al. [1]. Kamada et al. shows the number of adhesive objects increases with the shear rate by using multiple particle method. The number of grid points is 600 × 200, the volume fraction vf is 4% of the whole computational domain and the dc_c is 20 cells. Other calculation conditions are shown in Table 4.1. The sr means the shear rate at the wall. The Reynolds number is based on the channel height and freestream quantities. Figure 4.1 shows the instantaneous flow field with x-direction velocity distribution at t = 0.75. From Fig. 4.1, the objects are flown in the channel flow without any adhesion to the wall. It is caused that the y-direction velocity component is generated by only collision in this case. Sugiyama et al. shows that platelets do not adhere to the wall too. So there is a room for the reconsideration about a kind of method.

The channel flow including moving circular cylinders was investigated for the comparison the velocity distribution in the channel with the previous study by Skorczewski et al. [13]. The previous study by Skorczewski et al. confirmed that in the flow including platelets and red blood cells, the platelets move to the wall direction due to the red blood cells by using lattice Boltzmann method. We confirm this phenomenon too. In this case, the number of grid point is 400 × 200 and the dc_c is 20 cells. The boundary conditions are shown in Table 4.2. The Reynolds number is based on the channel height and freestream quantities. Figure 4.2 shows the velocity distribution of the channel flow. A good agreement is obtained with the previous result at the vf = 0 % from Fig. 4.2. However, in the case of the flow including objects, the error between the Skorczewski’s result and present result is generated. It is due to the difference of the flow including cylinders and red blood cells with deformation in the Skorczewski’s study.

The channel flow including an obstacle and moving circular cylinders was investigated to confirm the behavior of the moving cylinders in the channel flow around an obstacle. In this case, the number of grid points is 4000 × 1400 and the dc_c is 40 cells. The Reynolds number based on the diameter of the cylinder and volume fraction are set to be 5–50 and 30 %. Figure 4.3 shows the instantaneous flow field with x-direction velocity distribution. For Fig. 4.3, it is clear that the vertical velocity component toward the wall is generated behind the obstacle due to the vortex. In addition, the vortices generated behind the obstacle become larger with increasing the Reynolds number. Also, the velocity defect becomes also large. Furthermore, at the Reynolds number of 50, the large number of object adheres to the wall because a part of vortex structure separates to the downstream side due to the deposit of the objects on the obstacle. Figure 4.4 shows the relationship between the number of adhered objects (N_obj) and the distance from the center of the obstacle in x-direction (d_obs). For Fig. 4.4, it is considered that the number of adhered obstacles and the distance between adhered points become constant according as leaving from the obstacle due to the periodic vortex generation.

The channel flow including moving circular cylinders and dumbbell-like objects is investigated to confirm the behavior of the platelets and red blood cells in the channel flow. In this case, the number of grid points is 4000 × 1400, the dc_c is 40 cells and dc_r is 160 cells. The dc_r means the division number of the major axial of the red blood cells. The Reynolds number based on the channel height and volume fraction for the cylinders of the red blood cells were set to be 10 and 10%. Figures 4.5 and 4.6 show the instantaneous flow field with x-direction velocity distribution and the instantaneous flow field with x-direction velocity distribution. It is obtained that the velocity component to the wall direction is generated due to the disturbance that the red blood cells generate from Figs 4.5 and 4.6. In contrast, in the flow including only cylinders, the cylinders do not have y-direction velocity component from Fig. 4.1.

In this study, numerical method to handle the collision and adhesion based on immersed boundary method was developed and validated. In addition, we confirmed the behaviour of the objects in the various channel flows. The result showed that the method can be applied to a flow including various blood cell models with the collision and adhesion to simulate interaction of the blood cells and blood plasma. Moreover, it is confirmed that in the flow including only platelets, the platelets hardly adhere to the wall because the interaction between objects didn’t occur. Whereas, in the flow including platelets and red blood cells, the velocity component to the wall direction is generated due to the disturbance of the red blood cells. Furthermore, the platelets adhere to the wall because the vortex is generated behind the stent disturbs the flow. And at the high Reynolds number, the platelets adhesion can be observed in wide range due to the large separation and vortex generation.