Two-dimensional Simulation of Motion of Red Blood Cells with Deterministic Lateral Displacement Devices

Deterministic lateral displacement (DLD) technology has great potential for the separation, enrichment, and sorting of red blood cells (RBCs). This paper presents a numerical simulation of the motion of RBCs using DLD devices with different pillar shapes and gap configurations. We studied the effect of the pillar shape, row shift, and pillar diameter on the performance of RBC separation. The numerical results show that the RBCs enter “displacement mode” under conditions of low row-shift (∆λ < 1.4 µm) and “zigzag mode” with large row shift (∆λ > 1.5 µm). RBCs can pass the pillar array when the size of the pillar (d > 6 µm) is larger than the cell size. We show that these conclusions can be helpful for the design of a reliable DLD microfluidic device for the separation of RBCs.


Introduction
Human blood consists of plasma and mainly red and white blood cells and platelets [1]. Through continuous circulation, blood provides oxygen and nutrients, removes metabolic waste from tissue, regulates body pH and temperature, in addition to other biological functions. RBCs are the most critical component, and have a direct impact on hemodynamics and hemorheology. RBCs are highly deformable due to their biconcave, seedless, and highly flexible membrane [2]. Researchers has recognized the deformability of RBCs as an inherent indicator of diseases such as diabetes and malaria. The separation of deformable RBCs through Lab-on-a-chip techniques, which is based on the formers' mechanical properties, is becoming an essential process for medical research and clinical disease diagnosis [3][4][5]. Further understanding of the dynamic behavior of the RBCs in a microchannel, and exploring new methods by which to efficiently separate RBCs from plasma, are urgent to the development of biomedical engineering.
In medical research such as liquid biopsies and cytopathology, researchers have begun to focus on the direct isolation of target cells from whole blood. In 2018, a new method for continuously focusing and separating biological particles directly with shear-induced diffusion from whole blood was successfully demonstrated; the method effectively combines the inherent complexity of blood with the migration from the flow inertia while causing little pollution [6]. They also demonstrated a new multi-flow micro-fluidic (MFM)system for unlabeled circulating tumor cells (CTCs) from the peripheral blood of patients, and proposed a promising alternative method by which to realize CTC capture. Their results show that the method will offer the possibility of achieving chemical individuality treatment [7].
Moreover, several techniques for cell separation have been proposed, such as label-free discrimination and fractionation of cell populations [8], fluorescence-activated cell sorting (FACS), magnetic-activated cell sorting (MACS), and centrifugation separation. Fluorescence-activated cell

Simulation Methods and Models
We adopted the COMSOL ® Multiphysics software as the computational tool, which is relatively well-suited to the calculation of microspheres manipulation [21]. We optimized the simulation code of [22], in which the separation trajectories in a DLD device are accurately predicted using the finite element method. Furthermore, we studied the velocity field, surface stress, trajectories, and velocity of red blood cell movement. During the simulation, it was shown that when the height of the microspheres/cells is consistent with the thickness of the device, the geometric nature of the microspheres/cells can be simplified to two-dimensions [20]. Also, compared to the 3-D model, the 2-D model can greatly reduce the number of calculations. Additionally, the friction between the microspheres and the device wall can be neglected. Figure 1 shows a typical schematic of a DLD device. The microfluidic channel has one inlet on the left and two outlets on the right; the separated particles/cells will flow out of one of the outlets. To enhance the analysis, we chose two different pillar shapes (circle and the triangle) as the calculation cases. The detailed dimensions of the pillar posts are shown in Figure 1. We adopted the COMSOL ® Multiphysics software as the computational tool, which is relatively well-suited to the calculation of microspheres manipulation [21]. We optimized the simulation code of [22], in which the separation trajectories in a DLD device are accurately predicted using the finite element method. Furthermore, we studied the velocity field, surface stress, trajectories, and velocity of red blood cell movement. During the simulation, it was shown that when the height of the microspheres/cells is consistent with the thickness of the device, the geometric nature of the microspheres/cells can be simplified to two-dimensions [20]. Also, compared to the 3-D model, the 2-D model can greatly reduce the number of calculations. Additionally, the friction between the microspheres and the device wall can be neglected. Figure 1 shows a typical schematic of a DLD device. The microfluidic channel has one inlet on the left and two outlets on the right; the separated particles/cells will flow out of one of the outlets. To enhance the analysis, we chose two different pillar shapes (circle and the triangle) as the calculation cases. The detailed dimensions of the pillar posts are shown in Figure 1.  Table 1 lists the channel length, width, the pillar diameter, the height, and the obstacle spacing of the DLD device. Table 2 lists the physical properties of the fluid and the RBCs. Two working fluids were tested.   Table 1 lists the channel length, width, the pillar diameter, the height, and the obstacle spacing of the DLD device. Table 2 lists the physical properties of the fluid and the RBCs. Two working fluids were tested.  In this two-dimensional simulation, we need to calculate an incompressible liquid flow by solving the Navier-Stokes equations [23,24] and the continuity equation:

The DLD Model
where ρ ∂u f luid ∂t represents the unsteady inertia force (N/m 3 ), 12 represents the nonlinear inertia force (N/m 3 ), and F represents the volume force. ρ represents the fluid density (kg/m 3 ), p is the pressure (Pa), () represents the gradient operator, I is the unit diagonal matrix, µ f is the dynamic fluid viscosity (Pa·s), and d z represents the channel height (mm). When gravity or other volumetric forces are not considered, F = 0. u fluid represents the fluid velocity field (m/s). In a microfluidic system, the flow rate is small and the Reynolds number (R e << 100) is expressed as where l is the characteristic length of the rectangular channel (l = 109 µm), and U is the average velocity.

Boundary Conditions for Fluid-solid Interaction (FSI)
In a fluid-solid coupling boundary setting, the solid boundary is affected by the viscous force and flow pressure. Therefore, the velocity of the fluid can be used to derive the displacement rate of change of the solid microspheres. u fluid = u w (4) where Γ is the force on the solid boundary, which includes the fluid pressure and viscous resistance, and n denotes the outward normal bound. At the fluid-solid coupling boundary, the fluid velocity, u fluid , is equal to the rate of change of the solid displacement, u w . u solid represents the solid displacement field (m/s), u w represents the displacement rate of change of the solid microspheres and satisfies the no-slip condition at the channel wall, and σ is the Cauchy stress.

Initial Conditions
The fluid flow in the channel is fully developed laminar flow, and is driven by the pressure difference. To ensure that the nonlinear solver has the best possible convergence at the beginning, we set the inlet velocity to be non-constant. According to the parabolic characteristics of the laminar flow velocity distribution, we used the following formula to calculate the normal inlet velocity: where u 0 is the average flow velocity at the inlet (m/s), H is the channel height (mm), and Y is the value of the y-coordinate of the center of the microsphere/RBCs. At the exit, a fixed boundary condition, also known as the Dirichlet boundary condition, was used to determine the value of the pressure.

Red Blood Cell Model
Due to their non-expandable nature, RBCs can preserve their area and arc length in two dimensions. The deformability of cells depends on the elasticity of the cell membrane, the viscosity of the cytoplasm, and the applied flow rate [3,25]. Therefore, for the simulation of the cell screening process, it was necessary to consider the deformability of RBCs. Young's modulus, which is the modulus of elasticity detected by atomic force microscopy, is a physical quantity used to characterize the surface properties of cells in the biological field. It is determined only by the physical properties of the material itself; the larger Young's modulus, the greater the rigidity of the cells, and the more difficult it is to deform them.

Linear Elastic Material
If the rigid microspheres underwent small levels of deformation and were subjected to a low load, their displacement and deformation satisfy the linear elastic momentum conservation (see the Equation (11)) and the governing equilibrium equation [23,24]. In this simulation, the Young's modulus E = (0.25 ± 0.08) kPa, the Poisson's ratio v = 0.3 for RBCs.
Equation (10) represents Newton's equation of motion; Equation (11) is the linear elastic stress-strain law, and Equation (14) represents the strain-displacement (compatibility) equation. F v is the boundary force and J as the stiffness matrix. ε is the infinitesimal strain tensor (we also use ε to present the row shift fraction in the discussion section); finally, C is the stiffness matrix.

The Capillary Number
To consider the factors that affect cell deformation, such as cell membrane, cytoplasm, and the flow rate of the surrounding fluid, we introduced a dimensionless number, capillary number Ca, by combining the characteristics of these aspects [3].
The bending stiffness k b of a healthy RBC is 1 × 10 −19 and the effective radius R eff is 2.5 µm [1,3]. The maximum flow velocity U max between the pillars is from 10 µm/s to 10 mm/s, and the post-gap G = 1 µm. For healthy RBCs, the capillary number Ca range is [0.0375, 375], while the Ca range of pathological RBCs will increase nearly 10 times due to their stiffness, which means that the value is [0.0038, 37.5] [26]. Therefore, the deformation of RBCs can be judged by the capillary number Ca; the larger the capillary number Ca, the more easily the RBCs are deformed.

The Mesh and Independence Verification
After creating a mesh of the domain and discretizing the equation, we use the Any Lagrangian-Eulerian (ALE) method to describe the interface between the fluid and the RBCs. The geometry of the domain can be defined by moving the mesh, and the deformation and movement of the boundary are the same. COMSOL can calculate the new grid coordinates of the channel region based on the moving boundary and mesh smoothing of the structure, and solve the Navier-Stokes equations in the moving coordinate system in the mesh and fluid domain. The microsphere meshes in the domain move and deform freely with the simulated motion, and are automatically divided.
Also, we need to conduct grid independence verification to select the appropriate mesh. The most basic method to change the mesh size; we will choose this mesh scale when the difference between the results of two consecutive mesh scales reaches a point at which it can be ignored. Table 3 gives the corresponding degrees of freedom for solving the various equations at different grid scales.  Figure 2a plots the rate of change of displacement in the X direction for circular microspheres with t = 0 to t = 0.3 seconds at different grid scales. The grid is very close and the error is small in terms of regularity, refinement, ultra-fineness and extreme refinement. To consider the computational efficiency and the solver convergence, the ultra-thinning ratio is used for subsequent simulations. According to the COMSOL mesh quality detection method, the meshing is consistent with the rule that the mesh quality is close to 1, meaning that the mesh quality is good.

Results and Discussion
To simplify the simulation, we consider the RBCs to be microscopically deformed microspheres by ignoring their deformability in order to study the deflection effect of DLD devices on different  Figure 2b is a meshing diagram and shows enlarges the microsphere area. The mesh is divided into free triangle meshes, and the ultra-fine meshing method is selected. Inside the obstacle, the moving mesh will change with the deformation of the obstacle; in the outer boundary of the watershed, the deformation in all directions is set to zero. Therefore, the initial mesh at t = 0 s is not evenly divided, but it can be seen that the mesh is evenly distributed around the microspheres. The mesh is looser on the fluid domain, while the fluid-solid boundary has a denser and smaller cell grid.

Results and Discussion
To simplify the simulation, we consider the RBCs to be microscopically deformed microspheres by ignoring their deformability in order to study the deflection effect of DLD devices on different sized microspheres. Three kinds of microspheres with different diameters were used for numerical simulations; the trajectories of the microspheres flowing through the DLD microchannel were analyzed.

The Trajectories of Microspheres
In a DLD device, the microspheres will separate according to the arrangement of the pillars under the laminar flow at low Reynolds number.
From the Figure 3, there are two main modes in which the rigid microspheres move between the gaps of pillar arrays. It is interesting to note that the position of the pillars will change the trajectories of the microspheres, leading their trajectories into a "zigzag mode"; however, the trajectories return to the original state after three rows of obstacles. When the radius is greater than the width of the channel, the microspheres adopt "displacement mode", and may not be able to enter channel 1 of the fluid in the gap, thus behaving differently in the pillar array [27].

Results and Discussion
To simplify the simulation, we consider the RBCs to be microscopically deformed microspheres by ignoring their deformability in order to study the deflection effect of DLD devices on different sized microspheres. Three kinds of microspheres with different diameters were used for numerical simulations; the trajectories of the microspheres flowing through the DLD microchannel were analyzed.

The Trajectories of Microspheres
In a DLD device, the microspheres will separate according to the arrangement of the pillars under the laminar flow at low Reynolds number.
From the Figure 3, there are two main modes in which the rigid microspheres move between the gaps of pillar arrays. It is interesting to note that the position of the pillars will change the trajectories of the microspheres, leading their trajectories into a "zigzag mode"; however, the trajectories return to the original state after three rows of obstacles. When the radius is greater than the width of the channel, the microspheres adopt "displacement mode", and may not be able to enter channel 1 of the fluid in the gap, thus behaving differently in the pillar array [27].  Compared with the longitudinal rectangular microchannels described in reference [27], we used horizontal microchannels in the numerical simulation to analyze the trajectories of the microspheres. We studied the factors affecting the separation of microspheres in the DLD device, the row shifts fraction ε, and the diameter of the microspheres D. The row-shift fraction ε is expressed as As shown in Figure 4, ∆λ is row shift and λ is the center-to-center distance. There is a critical diameter D c between the "displacement mode" and the "zigzag mode". Based on the gap size and microsphere diameters, the critical diameter of the microspheres is [28], where G is the post-gap, indicating the distance between the two pillars. The result shows the occurrence of "displacement mode" at low row shifts (∆λ ≤ 2.5 µm) and "zigzag mode" for larger row shifts (∆λ ≥ 3 µm). This means that the separation trajectories of the microspheres is dependent on the row shifts in the presently-studied DLD device.
To investigate the effect of the complex pillar array on the separation of microspheres, we used different arrangements of pillars, i.e., by varying the row shift, the post-gap and the center-to-center distance between the two pillars. In Figure 4a, the flow pattern of microspheres 1.4 µm and 3 µm in diameter demonstrates "displacement mode" when ∆λ = 1.0 µm, D c = 1.95 µm, G = 4 µm, and ε = 0.1111. In Figure 4b, the trajectories of microspheres with the same diameter adopt "zigzag mode" when ∆λ = 4.0 µm, D c = 3.79 µm, G=4 µm, and ε = 0.4444. This indicates that the trajectories of the microspheres are related to the row shift, which is consistent with the conclusions found in the literature.
Compared with the longitudinal rectangular microchannels described in reference [27], we used horizontal microchannels in the numerical simulation to analyze the trajectories of the microspheres. We studied the factors affecting the separation of microspheres in the DLD device, the row shifts fraction ε, and the diameter of the microspheres D. The row-shift fraction ε is expressed as As shown in Figure 4, ∆λ is row shift and λ is the center-to-center distance.
There is a critical diameter Dc between the "displacement mode" and the "zigzag mode". Based on the gap size and microsphere diameters, the critical diameter of the microspheres is [28], where G is the post-gap, indicating the distance between the two pillars. The result shows the occurrence of "displacement mode" at low row shifts (∆λ ≤ 2.5 µm) and "zigzag mode" for larger row shifts (∆λ ≥ 3 µm). This means that the separation trajectories of the microspheres is dependent on the row shifts in the presently-studied DLD device.
To investigate the effect of the complex pillar array on the separation of microspheres, we used different arrangements of pillars, i.e., by varying the row shift, the post-gap and the center-to-center distance between the two pillars. In Figure 4a, the flow pattern of microspheres 1.4 µm and 3 µm in diameter demonstrates "displacement mode" when ∆λ = 1.0 µm , Dc = 1.95 µm, G = 4 µm, and ε = 0.1111. In Figure 4b, the trajectories of microspheres with the same diameter adopt "zigzag mode" when ∆λ = 4.0 µm, Dc = 3.79 µm, G=4 µm, and ε = 0.4444. This indicates that the trajectories of the microspheres are related to the row shift, which is consistent with the conclusions found in the literature.  Compared to [20], we found that the diameter of the microspheres affects their trajectories; we also selected different diameters of the microspheres in our simulation to verify this result. Three different diameters of microspheres were released from the same position in order to analyze the trajectories when staggered pillars were used. From top to bottom, the microspheres have diameters of 1.0 µm, 0.4 µm and 0.2 µm with the ∆λ = 1.5 µm, Dc = 1.88 µm, G = 3 µm, and ε = 0.1875 respectively, and the trajectories of the three microspheres is in "zigzag mode". When ∆λ = 3.5 µm, Dc = 2.82 µm, G = 3 µm, and ε = 0.4375, it was found that microspheres with diameters of 0.1 µm, 0.2 µm, 0.4 µm, 1 µm and 2 µm have different displacement rates, but that the trajectories are the same in "zigzag mode". Therefore, in this arrangement of pillars in Figure 5, the diameter of the microspheres does Compared to [20], we found that the diameter of the microspheres affects their trajectories; we also selected different diameters of the microspheres in our simulation to verify this result. Three different diameters of microspheres were released from the same position in order to analyze the trajectories when staggered pillars were used. From top to bottom, the microspheres have diameters of 1.0 µm, 0.4 µm and 0.2 µm with the ∆λ = 1.5 µm, D c = 1.88 µm, G = 3 µm, and ε = 0.1875 respectively, and the trajectories of the three microspheres is in "zigzag mode". When ∆λ = 3.5 µm, D c = 2.82 µm, G = 3 µm, and ε = 0.4375, it was found that microspheres with diameters of 0.1 µm, 0.2 µm, 0.4 µm, 1 µm and 2 µm have different displacement rates, but that the trajectories are the same in "zigzag mode". Therefore, in this arrangement of pillars in Figure 5, the diameter of the microspheres does not change their trajectories; therefore, we suspect that this indicates a new law between "zigzag mode" and "displacement mode" regarding row shifts in these pillar arrangements.
To explore the critical row shifts between the "displacement mode" and the "zigzag mode" in this pillar arrangement, the same diameters, i.e., 1.0 µm, 0.4 µm and 0.2 µm, of the microspheres are in ∆λ = 1.0 µm, D c = 1,55 µm, G = 3 µm, and ε = 0.125 and another in ∆λ = 1.4, D c = 1.82 µm, G = 3 µm, and ε = 0.175 again in Figure 6. The results show that they all lead to "displacement mode", and that in the arrangement of the pillars, the "displacement mode" shifts at the low row shifts (∆λ ≤ 1.4 µm), and the "zigzag mode" shifts at the larger row shifts (∆λ ≥ 1.5 µm). With the trajectories of microspheres changing from "zigzag mode" to "displacement mode", the critical diameters become smaller, and the row shift fraction becomes larger.
Micromachines 2019, 10, x FOR PEER REVIEW 9 of 15 not change their trajectories; therefore, we suspect that this indicates a new law between "zigzag mode" and "displacement mode" regarding row shifts in these pillar arrangements. To explore the critical row shifts between the "displacement mode" and the "zigzag mode" in this pillar arrangement, the same diameters, i.e., 1.0 µm, 0.4 µm and 0.2 µm, of the microspheres are in ∆λ = 1.0 µm, Dc = 1,55 µm, G = 3 µm, and ε = 0.125 and another in ∆λ = 1.4, Dc = 1.82 µm, G = 3 µm, and ε = 0.175 again in Figure 6. The results show that they all lead to "displacement mode", and that in the arrangement of the pillars, the "displacement mode" shifts at the low row shifts (∆λ ≤ 1.4 µm), and the "zigzag mode" shifts at the larger row shifts (∆λ ≥ 1.5 µm). With the trajectories of microspheres changing from "zigzag mode" to "displacement mode", the critical diameters become smaller, and the row shift fraction becomes larger. The microspheres showed periodic variations in this DLD device. We mainly studied microsphere separation in "zigzag mode" using staggered pillars. In this section, we considered RBCs as being rigid, and simulated the microspheres with the same parameters as those of the RBCs. We believe that this arrangement can effectively predict the trajectory of RBCs during separation. We not change their trajectories; therefore, we suspect that this indicates a new law between "zigzag mode" and "displacement mode" regarding row shifts in these pillar arrangements. To explore the critical row shifts between the "displacement mode" and the "zigzag mode" in this pillar arrangement, the same diameters, i.e., 1.0 µm, 0.4 µm and 0.2 µm, of the microspheres are in ∆λ = 1.0 µm, Dc = 1,55 µm, G = 3 µm, and ε = 0.125 and another in ∆λ = 1.4, Dc = 1.82 µm, G = 3 µm, and ε = 0.175 again in Figure 6. The results show that they all lead to "displacement mode", and that in the arrangement of the pillars, the "displacement mode" shifts at the low row shifts (∆λ ≤ 1.4 µm), and the "zigzag mode" shifts at the larger row shifts (∆λ ≥ 1.5 µm). With the trajectories of microspheres changing from "zigzag mode" to "displacement mode", the critical diameters become smaller, and the row shift fraction becomes larger. The microspheres showed periodic variations in this DLD device. We mainly studied microsphere separation in "zigzag mode" using staggered pillars. In this section, we considered RBCs as being rigid, and simulated the microspheres with the same parameters as those of the RBCs. We believe that this arrangement can effectively predict the trajectory of RBCs during separation. We The microspheres showed periodic variations in this DLD device. We mainly studied microsphere separation in "zigzag mode" using staggered pillars. In this section, we considered RBCs as being rigid, and simulated the microspheres with the same parameters as those of the RBCs. We believe that this arrangement can effectively predict the trajectory of RBCs during separation. We analyzed the velocity and surface stress of different microsphere shapes (circular, elliptical, triangular and diamond) and the effect of the shape of the pillars.

The Surface Stress on the Microsphere
The deformability of RBCs is mainly due to stress on the surface. Next, we will discuss the stress induced by the flow. RBCs can be damaged or even rupture under complex loading forces, including viscous resistance and the pressure of the fluid. Therefore, it is necessary to analyze the surface stress on the microspheres in order to determine whether the microspheres will remain intact in DLD devices with a pillar structure [29].
The Young's modulus causes elastic deformation by affecting the level of stress on the microspheres. The surface stress of the microspheres is scalar, calculated from a stress tensor, and can be used to determine whether the microspheres remain intact when subjected to the loading force. In Figure 7a, when t = 0 s, at 0.126 s, 0.156 s, 0.171 s, and 0.19 s, the position of the microspheres of different shapes is shown in different colors. The maximum and minimum stress points are also presented (Figure 7b,c). When the microspheres are far from the pillars in the DLD device, the stress distribution is almost constant and uniform. The surface stress increases with increasing pressure. Circular and triangular microspheres are always symmetrical when they are in movement, so the maximum surface and minimum surface stress values are close to each other. Because of their different long and short axes, and the fact that microspheres collide with the pillars in the DLD device, the surface stress of elliptical and rhomboid microspheres will change abruptly. Under the low elasticity of the microspheres, surface stress is insufficient to cause apparent deformation. analyzed the velocity and surface stress of different microsphere shapes (circular, elliptical, triangular and diamond) and the effect of the shape of the pillars.

The Surface Stress on the Microsphere
The deformability of RBCs is mainly due to stress on the surface. Next, we will discuss the stress induced by the flow. RBCs can be damaged or even rupture under complex loading forces, including viscous resistance and the pressure of the fluid. Therefore, it is necessary to analyze the surface stress on the microspheres in order to determine whether the microspheres will remain intact in DLD devices with a pillar structure [29]. The Young's modulus causes elastic deformation by affecting the level of stress on the microspheres. The surface stress of the microspheres is scalar, calculated from a stress tensor, and can be used to determine whether the microspheres remain intact when subjected to the loading force. In Figure 7a, when t = 0 s,at 0.126 s, 0.156 s, 0.171 s, and 0.19 s, the position of the microspheres of different shapes is shown in different colors. The maximum and minimum stress points are also presented (Figure 7b,c). When the microspheres are far from the pillars in the DLD device, the stress distribution is almost constant and uniform. The surface stress increases with increasing pressure.

The Velocity of the Microspheres
According to the five moments shown in Figure 7 above, we found that the velocity of the microspheres is only related to the viscous resistance at the entrance region. Since the viscous resistance is negligible, we consider the initial velocity of the microspheres to be close to the inlet velocity of the fluids (u 0 = 12 µm/s).
Next, we observed the velocity of the liquid in the DLD device. High velocities were measured in the gaps between the pillars. From the velocity of the microspheres, it was observed that the initial velocity of the microspheres was close to the inlet flow velocity, and in the case of a low level of velocity resistance, the shear stress of the fluid on the microspheres was negligible. This confirms that the microspheres have excellent follow-up performance in DLD devices. The result shows that the pillar structure can increase the velocity of the microspheres.
From Figure 8, we can see when the microspheres pass theough the DLD device, their velocity increases, with the maximum velocity occuring in the gaps between the pillar posts. When flowing between the pillars, the velocity was lower; this is caused by the immense surface stress required to produce negative resistance. When the microspheres approach the outlet, the velocity will gradually converge, eventually becoming zero.
Micromachines 2019, 10, x FOR PEER REVIEW 11 of 15 maximum surface and minimum surface stress values are close to each other. Because of their different long and short axes, and the fact that microspheres collide with the pillars in the DLD device, the surface stress of elliptical and rhomboid microspheres will change abruptly. Under the low elasticity of the microspheres, surface stress is insufficient to cause apparent deformation.

The Velocity of the Microspheres
According to the five moments shown in Figure 7 above, we found that the velocity of the microspheres is only related to the viscous resistance at the entrance region. Since the viscous resistance is negligible, we consider the initial velocity of the microspheres to be close to the inlet velocity of the fluids (u0 = 12 µm/s).
Next, we observed the velocity of the liquid in the DLD device. High velocities were measured in the gaps between the pillars. From the velocity of the microspheres, it was observed that the initial velocity of the microspheres was close to the inlet flow velocity, and in the case of a low level of velocity resistance, the shear stress of the fluid on the microspheres was negligible. This confirms that the microspheres have excellent follow-up performance in DLD devices. The result shows that the pillar structure can increase the velocity of the microspheres.
From Figure 8, we can see when the microspheres pass theough the DLD device, their velocity increases, with the maximum velocity occuring in the gaps between the pillar posts. When flowing between the pillars, the velocity was lower; this is caused by the immense surface stress required to produce negative resistance. When the microspheres approach the outlet, the velocity will gradually converge, eventually becoming zero.

The Effect of Pillar Shapes
We investigated the effect of circular and triangular pillar structures on the flow characteristics in DLD devices. Figure 9 shows that there are some zero-velocity areas between the two pillars in Figure 8. The velocity of the microspheres changes over time. Four kinds of microspheres with circular, elliptical, triangular, and diamond shapes were compared with the data of the reference [22], and a series of time points were selected to analyze their velocity: t = 0 s, 0126 s, 0.156 s, 0.171 s, and 0.19 s.

The Effect of Pillar Shapes
We investigated the effect of circular and triangular pillar structures on the flow characteristics in DLD devices. Figure 9 shows that there are some zero-velocity areas between the two pillars in these two shapes of pillar posts. We found that the fluid velocity in the circular pillar posts was much greater than that of the triangular ones. We thus optimize the shape of the pillars and analyzed the effect of velocity on the motion and stress levels of the microspheres.
The paraments in these DLD devices are all in ∆λ = 3.5 µm, D c = 2.82 µm, G = 3 µm, and ε = 0.4375 in Figure 10. The liquid gap between the microspheres and the pillars are clear, while the boundary is slightly weaker. We can see that the microspheres tend to move closer to the pillars in the triangular DLD device, and that the stress level is uniform. This indicates that triangular pillars can increase the levels of interaction with the microspheres, significantly affecting their trajectories.
Micromachines 2019, 10, x FOR PEER REVIEW 12 of 15 these two shapes of pillar posts. We found that the fluid velocity in the circular pillar posts was much greater than that of the triangular ones. We thus optimize the shape of the pillars and analyzed the effect of velocity on the motion and stress levels of the microspheres. The paraments in these DLD devices are all in ∆λ = 3.5 µm, Dc = 2.82 µm, G = 3 µm, and ε = 0.4375 in Figure 10. The liquid gap between the microspheres and the pillars are clear, while the boundary The paraments in these DLD devices are all in ∆λ = 3.5 µm, Dc = 2.82 µm, G = 3 µm, and ε = 0.4375 in Figure 10. The liquid gap between the microspheres and the pillars are clear, while the boundary is slightly weaker. We can see that the microspheres tend to move closer to the pillars in the triangular DLD device, and that the stress level is uniform. This indicates that triangular pillars can increase the levels of interaction with the microspheres, significantly affecting their trajectories. According to the setting of Young's modulus and Poisson's ratio, it was found from the simulation that the deformation of the microspheres was almost consistent with that of RBCs. To better observe RBC deformation in the DLD device with triangular pillars, we decided to combine the RBC model with the triangular pillars to attain the trajectory of the RBCs in the DLD device. And the Figure 11 shows the schematic diagram of the RBCs in 2-D and 3-D.

The Deformation of Red Blood Cells
We analyzed the deformability of RBCs; asphericity needs to be introduced her eto describe the deviation from the spheric geometry [28], where λ1 and λ2 are the square roots of two non-zero eigenvalues of the radius of the rotation tensor. The δ from 0 to 1 reflects the extent to which perfect RBCs or microspheres changed in terms of elongation. In the 2-D model, δ equalled 0.29. A DLD device with triangular pillars and with a gap of s = 3 µm and height of D = 5 µm was selected to analyze the level of deformation of RBCs. It was found that there was no clear flow layer between the RBC model and the pillars, indicating that the DLD device has considerable sensitivity. In Figure 12a, the RBCs stayed near the top of the pillar for a long time, and could not continue to move. To explore the relationship between the trajectories of RBCs and the pillar scale, we compared the size of the pillars with a gap of s = 10 µm and a height of D = 15 µm.  In Figure 12b, it can see that the RBCs are able to pass through the DLD device of the pillar structure when D > 6 µm. Simultaneously, by increasing the diameter and the gap of the pillars, the flow velocity of the pillars decreases from the top to the bottom. The velocity of the RBCs near the pillars is less affected by the flow field, and the deformation is negligible.
To characterize the separation model of RBCs, the function of row shift fraction ε and capillary number Ca was used to verify the trajectories of the microspheres. It was found that the boundary between the displacement model and the "zigzag model" was very close to a simple exponential function [19].

The Deformation of Red Blood Cells
We analyzed the deformability of RBCs; asphericity needs to be introduced her eto describe the deviation from the spheric geometry [28], where λ 1 and λ 2 are the square roots of two non-zero eigenvalues of the radius of the rotation tensor. The δ from 0 to 1 reflects the extent to which perfect RBCs or microspheres changed in terms of elongation. In the 2-D model, δ equalled 0.29. A DLD device with triangular pillars and with a gap of s = 3 µm and height of D = 5 µm was selected to analyze the level of deformation of RBCs. It was found that there was no clear flow layer between the RBC model and the pillars, indicating that the DLD device has considerable sensitivity. In Figure 12a, the RBCs stayed near the top of the pillar for a long time, and could not continue to move. To explore the relationship between the trajectories of RBCs and the pillar scale, we compared the size of the pillars with a gap of s = 10 µm and a height of D = 15 µm.

The Deformation of Red Blood Cells
We analyzed the deformability of RBCs; asphericity needs to be introduced her eto describe the deviation from the spheric geometry [28], where λ1 and λ2 are the square roots of two non-zero eigenvalues of the radius of the rotation tensor. The δ from 0 to 1 reflects the extent to which perfect RBCs or microspheres changed in terms of elongation. In the 2-D model, δ equalled 0.29. A DLD device with triangular pillars and with a gap of s = 3 µm and height of D = 5 µm was selected to analyze the level of deformation of RBCs. It was found that there was no clear flow layer between the RBC model and the pillars, indicating that the DLD device has considerable sensitivity. In Figure 12a, the RBCs stayed near the top of the pillar for a long time, and could not continue to move. To explore the relationship between the trajectories of RBCs and the pillar scale, we compared the size of the pillars with a gap of s = 10 µm and a height of D = 15 µm.  In Figure 12b, it can see that the RBCs are able to pass through the DLD device of the pillar structure when D > 6 µm. Simultaneously, by increasing the diameter and the gap of the pillars, the flow velocity of the pillars decreases from the top to the bottom. The velocity of the RBCs near the pillars is less affected by the flow field, and the deformation is negligible.
To characterize the separation model of RBCs, the function of row shift fraction ε and capillary number Ca was used to verify the trajectories of the microspheres. It was found that the boundary between the displacement model and the "zigzag model" was very close to a simple exponential function [19]. In Figure 12b, it can see that the RBCs are able to pass through the DLD device of the pillar structure when D > 6 µm. Simultaneously, by increasing the diameter and the gap of the pillars, the flow velocity of the pillars decreases from the top to the bottom. The velocity of the RBCs near the pillars is less affected by the flow field, and the deformation is negligible.
To characterize the separation model of RBCs, the function of row shift fraction ε and capillary number Ca was used to verify the trajectories of the microspheres. It was found that the boundary between the displacement model and the "zigzag model" was very close to a simple exponential function [19]. ε(Ca) = 2.9 + 5.4e −1.72C a According to Equation 20, ε(Ca) = 3 µm, Ca = 2.1 within the "zigzag mode" and the trajectory of RBCs is "zigzag mode" in the row shift ε = 3 µm. ε(Ca) = 4 µm and Ca = 0.9 in "displacement mode" and when the trajectory of RBCs is in "displacement mode" in the row shift fraction ε = 4 µm. Following Ref. [19], we know that the softer and more rigid gradient of the RBC yield Ca ranges from 0.4 to 1.2.

Conclusions
We studied RBC separation with a DLD device using the finite element simulation. This research will provide an essential reference for further research of RBC separation using DLD devices in the future. The following conclusions were obtained: (a) The trajectories of the RBCs are related to the row shifts. We found that the RBCs enter "displacement mode" under conditions of low row-shift (∆λ < 1.4 µm) and "zigzag mode" with large row-shift (∆λ > 1.5 µm). (b) The velocity and the surface stress of the microspheres are related to the shape of the pillars. The triangular pillars will produce high velocities and the uniform stress on the microspheres, which will enhance separation. (c) By considering the deformation of RBCs, the row shift fraction ε and the capillary number Ca are used to determine the RBC separation mode.