A Numerical Study on the Erythrocyte Flow Path in I-Shaped Pillar DLD Arrays

: Erythrocyte enrichment is needed for blood disease diagnosis and research. DLD arrays with an I-shaped pillar (I-pillar) sort erythrocytes in a unique, accurate, and low-reagent method. However, the existing I-shaped pillar DLD arrays for erythrocyte sorting have the drawbacks of higher ﬂow resistance and more challenging fabrication. A two-dimensional erythrocyte simulation model and the arbitrary Lagrangian–Euler equations at the cell–ﬂuid boundary were built based on the ﬂuid–solid coupling method to investigate the inﬂuencing factors of the erythrocyte ﬂow path in an I-pillar DLD array and ﬁnd its optimization method. Three different sizes of I-pillars were built and multiple sets of corresponding arrays were constructed, followed by ﬁnite element simulations to separately investigate the effects of these arrays on the induction of erythrocyte motion paths. This work demonstrates the motion paths of erythrocyte models in a series of I-pillar arrays with different design parameters, aiming to summarize the variation modes of erythrocyte motion paths, which in turn provides some reference for designing and optimizing the pillar size and array arrangement methods for I-pillar array DLD chips.


Introduction
Erythrocytes are an essential component of about 45% of human blood. Their deformability directly affects the hemorheology [1][2][3]. Erythrocyte biconcave shape and corresponding deformability are essential features of their biological function, and the deformability of erythrocytes can be seriously affected under the influence of various diseases [4][5][6]. This deformable change has been identified as an intrinsic indicator of some blood disorders. Therefore, sorting and enriching erythrocytes is necessary in researching and diagnosing various blood diseases. As a new technology, in recent years, the microfluidic method can potentially apply tiny biological particles. Instead of the existing bulky and complex blood-cell processing equipment, microfluidic technology is applied to sorting and enriching erythrocytes, which can provide a new, efficient, cheap method for studying blood diseases and erythrocyte properties.
A microfluidic technique called deterministic lateral displacement (DLD) can separate and enrich micro-nano particles of various sizes [7,8]. An array of micro pillars is arranged in a flat microchannel to form a DLD chip, and the array is arranged at an angle to the main channel. The particle suspension flowing through the array is divided into a large number of microfluidic beams, and the particles are affected by the combined action of these beams and the pillars during their motion. The larger particles move in displacement flow mode along the pillar arrangement toward the side of the channel, while the smaller particles move in zigzag mode around the pillars in the same overall direction as the main channel. This eventually causes the two types of particles to be collected at different exits at the end of the chip [9]. movement of red blood cells, and provides new and improved solutions for the design of I-pillar array.

Calculation Method
We adopted version 5.6 of COMSOL ® Multiphysics software as the calculation tool, which is widely used for the computation of particle trajectories [22]. Typically, the I-pillar DLD arrays used for erythrocyte sorting are very narrow and long. In addition, the study of erythrocyte motion requires the longest possible arrays. The complexity of the threedimensional (3D) problem in simulating the motion and deformation of particles in a DLD sorting chip stems from more than just the coupling itself. The interaction between the particles and the fluid and microcolumns can cause the particles to tumble and deform during their motion. This can lead to mesh distortion and overlap in the computational domain, which makes the computation non-convergent. To reduce the computational demand, a 2D approximation is usually used to simplify the 3D problem. Therefore, 2D arrays and erythrocyte models were established to simulate longer arrays. The proposed array height was 15 µm, and the erythrocyte model was driven by the fluid flow through ten rows of pillars leading to the outlet. Figure 1 shows a partial schematic diagram of the I-pillar DLD chip, in which the red blood cells follow the suspension from the left inlet into the array and tumble under the combined effect of the fluid and the pillar. When flowing out from the array, the red cells are displaced in the direction perpendicular to the suspension flow. The main parameters affecting the effect of the array and pillar action are shown in Figure 1, where H is the side length of the I-pillar; r is the radius of the semicircle at the groove of the pillar; G x and G y denote the lateral and downstream pillar distances; and λ is the row shift distance of the array.

Fluid Flow Solver
these arrays were observed. The comparison of these motion paths reflects the effect of parameters such as pillar size, distance between array rows, and fluid flow rate on the movement of red blood cells, and provides new and improved solutions for the design of I-pillar array.

Calculation Method
We adopted version 5.6 of COMSOL ® Multiphysics software as the calculation tool, which is widely used for the computation of particle trajectories [22]. Typically, the I-pillar DLD arrays used for erythrocyte sorting are very narrow and long. In addition, the study of erythrocyte motion requires the longest possible arrays. The complexity of the threedimensional (3D) problem in simulating the motion and deformation of particles in a DLD sorting chip stems from more than just the coupling itself. The interaction between the particles and the fluid and microcolumns can cause the particles to tumble and deform during their motion. This can lead to mesh distortion and overlap in the computational domain, which makes the computation non-convergent. To reduce the computational demand, a 2D approximation is usually used to simplify the 3D problem. Therefore, 2D arrays and erythrocyte models were established to simulate longer arrays. The proposed array height was 15 µm, and the erythrocyte model was driven by the fluid flow through ten rows of pillars leading to the outlet. Figure 1 shows a partial schematic diagram of the I-pillar DLD chip, in which the red blood cells follow the suspension from the left inlet into the array and tumble under the combined effect of the fluid and the pillar. When flowing out from the array, the red cells are displaced in the direction perpendicular to the suspension flow. The main parameters affecting the effect of the array and pillar action are shown in Figure 1, where H is the side length of the I-pillar; r is the radius of the semicircle at the groove of the pillar; Gx and Gy denote the lateral and downstream pillar distances; and λ is the row shift distance of the array. In terms of the selection of chip material parameters, it was assumed that the DLD chip was composed of polydimethylsiloxane (PDMS), which is commonly used in microfluidic experiments, and its Young's modulus was set as 1.5 MPa, Poisson's ratio as 0.49, and dynamic viscosity as 0.001 Pa·s [21,23,24]. The physical parameters of the proposed chip material, as well as the erythrocyte suspension, were all constants.

Fluid Flow Solver
According to the law of conservation of mass, the continuity equation satisfied by the sample solution is: In terms of the selection of chip material parameters, it was assumed that the DLD chip was composed of polydimethylsiloxane (PDMS), which is commonly used in microfluidic experiments, and its Young's modulus was set as 1.5 MPa, Poisson's ratio as 0.49, and dynamic viscosity as 0.001 Pa·s [21,23,24]. The physical parameters of the proposed chip material, as well as the erythrocyte suspension, were all constants.
According to the law of conservation of mass, the continuity equation satisfied by the sample solution is: where ρ f is the density of the fluid and u f is the fluid velocity. The erythrocyte suspension (saline, phosphate buffer, etc.) usually used in the process of DLD sorting experiments has a linear change in velocity gradient in the direction perpendicular to the flow, and can be considered as a Newtonian fluid. At the same time, the flow of fluid in the DLD chip is slow, the Reynolds number is very small (Re << 10) [25], and the flow state belongs to the laminar flow. According to the nature of the actual fluid, the fluid was set as an incompressible steady-state laminar fluid motion during the simulation. At this time, ρ f and u f were constants, and Equation (1) could be simplified as: According to the law of conservation of momentum, the flow of fluid in the DLD sorting chip satisfied the Navier-Stokes (N-S) equation [26]: is the inertial force per unit volume of fluid, −∇p represents the pressure gradient, µ f ∇ 2 u f denotes the fluid's deformation stress, − 2 3 ∇(µ f ∇ · u f ) represents the expansion stress of the viscous body, and ρ f F is the mass force per unit volume of fluid.
The erythrocyte sample solution used for the DLD separation chip was in laminar flow state, where the viscous force had greater influence on the flow incompressible field than the inertial force [27]. This fluid was considered an incompressible flow, and its density and viscosity were uniform throughout the process. Equation (3) can be simplified as: The other feature sizes in the flow field must be significantly less than the depth of the flow field in order for the two-dimensional simulation to be applied. In contrast, the main channel of a DLD chip must be generally flat, and its depth must be significantly smaller than its width. As a result, the simulation approach in this research applied the shallow groove approximation theory [28]. With this approach, the volume force representing the effect of the flat channel inner surface was added to the N-S equation: In the transient flow simulation, the entire flow field changed with time, and the geometric boundaries and the grid in the physical plane changed with time. Considering the change in the grid boundary, the (ξ, η) coordinate system was introduced, and this coordinate system was used as the computational space. The grid in the computational space corresponded to the grid in the physical space one by one. When the grid in the physical space changed with time, the grid in the computational space was invariant, and at this time, the coordinate transformation quantity relationship was: J in Equation (6) is the Jacobian determinant: The transformation relation of the volume when the corresponding grid was solved in two spaces was: dxdy = Jdξdη (8) By bringing this relationship into the N-S equation and continuously calculating the value of J as the simulation proceeded, the physical space could be solved by calculating the differential equation in the space.
The sample solution in the DLD separation chip was driven by the injection pump, and the flow field in the chip reached a stable state after fully developing flow conditions. Therefore, the average flow velocity at the chip inlet was set as u 0 , and the pressure and viscous resistance at the outlet were defined according to Dirichlet conditions: The injection pump injected the sample solution into the DLD chip. The pressure difference was formed between the inlet and outlet of the chip. With the pressure difference, the driving particles arrived at the outlet along with the flow of the fluid. The inner wall of the DLD sorting chip and the wall of the pillar obstacle were uniformly set as no-slip conditions. This paper does not consider the influence of the material properties of the pillar or wall on the motion trajectory of the object. If the acceleration motion of the particle at the entrance of the chip is considered, the initial values of its displacement field u s and velocity field ∂u s /∂t are set to zero. If the particle has moved to a certain part of the chip, the initial value of velocity field ∂u s /∂t can be set to the velocity at the corresponding position under the steady-state flow field.

Erythrocyte Structure Solver
Evans and Feng [29] fitted an expression for the geometry of human erythrocytes based on empirical data: where 0 ≤ ω ≤ π, c 0 = 0.207, c 1 = 2.003, c 2 = −1.123, and R = 4 µm is the radius of erythrocyte. Then the obtained curves were moved and scaled to establish the approximate shape of the erythrocytes with the cell membranes and cytoplasm, and the geometric shape was obtained as shown in Figure 2. Erythrocytes can be thought of as an incompressible fluid encased in an elast brane that has a distinct elastic deformability and is primarily composed of a 5 n phospholipid bilayer and a cytoskeleton embedded in the lipid layer [30]. If the cyte model structure composed of the cell membrane and cytoplasm was used u  Erythrocytes can be thought of as an incompressible fluid encased in an elastic membrane that has a distinct elastic deformability and is primarily composed of a 5 nm thick phospholipid bilayer and a cytoskeleton embedded in the lipid layer [30]. If the erythrocyte model structure composed of the cell membrane and cytoplasm was used under the current continuum 2D model conditions, the degree of folding deformation created during the flow simulation was too great and a significant discrepancy with the actual experiment existed, as shown in Figure 3. This gap was brought about by structural flaws in the 2D model; hence the erythrocyte model had to be modified to account for these flaws. Erythrocytes can be thought of as an incompressible fluid encased in an elastic mem brane that has a distinct elastic deformability and is primarily composed of a 5 nm thick phospholipid bilayer and a cytoskeleton embedded in the lipid layer [30]. If the erythro cyte model structure composed of the cell membrane and cytoplasm was used under the current continuum 2D model conditions, the degree of folding deformation created dur ing the flow simulation was too great and a significant discrepancy with the actual exper iment existed, as shown in Figure 3. This gap was brought about by structural flaws in the 2D model; hence the erythrocyte model had to be modified to account for these flaws. The deformation of red blood cells during movement in the I-pillar array was smal because the force during movement was low. Under a low shear rate and low fluid vis cosity, the cell appeared as a solid [31], and its deformation and displacement satisfied the linear elastic momentum equation. Jiao et al. [21] verified the reliability of linearly elasti solids. Therefore, the whole erythrocyte model was simplified to a homogeneous elasti solid without a cell membrane by combining the Young's modulus, E = 1 kPa, and Pois son's ratio, 0.3 [21,32] of erythrocytes.
The linearly elastic material is characterized by the same deformability in all direc tions, and the controlling equation of its solid mechanics is: Figure 3. Deformation of 2D erythrocyte model with cytoplasmic structure.
The deformation of red blood cells during movement in the I-pillar array was small because the force during movement was low. Under a low shear rate and low fluid viscosity, the cell appeared as a solid [31], and its deformation and displacement satisfied the linear elastic momentum equation. Jiao et al. [21] verified the reliability of linearly elastic solids. Therefore, the whole erythrocyte model was simplified to a homogeneous elastic solid without a cell membrane by combining the Young's modulus, E = 1 kPa, and Poisson's ratio, 0.3 [21,32] of erythrocytes.
The linearly elastic material is characterized by the same deformability in all directions, and the controlling equation of its solid mechanics is: where Equation (13) is the strain-displacement equation, Equation (14) is the Newtonian equation of motion, and Equation (15) is the linear elastic stress-strain equation, and ε s is the Cauchy stress tensor, C is the stiffness matrix, E s is the particle's Young's modulus, and R s is the Poisson's ratio of the particle.

Coupling of Fluid and Erythrocyte Solvers
During the numerical calculation, the bidirectional coupling between the fluid and the solid was achieved through the interaction of the solid-liquid boundary [33,34]. On the boundary of the red cell model, the flow velocity of the fluid was equal to the motion velocity of the solid, and the boundary conditions can be expressed as: where f p is the force of the fluid on the solid particle, u s is the distance the particle travels, and u w is the particle velocity. This calculation process firstly solves the parameters in the flow field to obtain the pressure and stress; secondly, the force is applied to the particle boundary to calculate the boundary motion of the particle; finally, the displacement and deformation of the particle are obtained and fed back to the flow field to repeat the previous process. Therefore, the relationship between the load on the particle boundary and the load applied by the fluid is [35]: where the scale factors of the cell grid in the solid and fluid domains are represented by dv and dV, respectively. In the domains with free displacement, the deforming domain feature and the free features in the moving mesh and deformed geometry interface solve an equation for the mesh displacement. This equation smoothly deforms the mesh, given the constraints placed on the boundaries.
Among the four types of mesh smoothing offered by version 5.6 of COMSOL ® Multiphysics software, Yeoh smoothing is considered to generally produce the best results and allow for maximum boundary displacement before the mesh elements become inverted. It is more suitable for particle flow simulation within complex structural channels.

Model Verification and Flow Simulation
At the beginning of the simulation flow, as shown in Figure 4, the erythrocyte model was positioned close to the first pillar parallel to the flow direction at the start of the simulation flow to examine the rolling action. To match the fluid's velocity, its starting velocity was adjusted to u 0 . For the convenience of grid independence verification later, the point P on the erythrocyte model's edge was designated at its beginning position.
During the establishment of the initial grid, to ensure the accuracy of the flow pos of the erythrocyte model, the erythrocyte model and the rest were divided into diffe free triangular grids. Then the curvature factor, the maximum cell growth rate, and o parameters were adjusted to construct multiple groups of initial grids for compari Skewness (based on the equiangular skew that penalizes elements with large or smal gles as compared to the angles in an ideal element) was used as a reference for grid qua Table 1 shows the partitioning methods of different mesh scales.
During the establishment of the initial grid, to ensure the accuracy of the flow posture of the erythrocyte model, the erythrocyte model and the rest were divided into different free triangular grids. Then the curvature factor, the maximum cell growth rate, and other parameters were adjusted to construct multiple groups of initial grids for comparison. Skewness (based on the equiangular skew that penalizes elements with large or small angles as compared to the angles in an ideal element) was used as a reference for grid quality. Table 1 shows the partitioning methods of different mesh scales. Figure 5 plots the movement trajectory of point P on the erythrocyte model at different grid scales within 5 ms. The trajectory of point P was roughly the same when the number of grids was 164,491 and 255,303. In order to consider the computational efficiency and the convergence of the solver, the grid division method was adopted when the number of grids was 164,491. At this time, the mesh division conformed to the law that the mesh quality is close to one; that is, the mesh quality was good.  Figure 5 plots the movement trajectory of point P on the erythrocyte model a ent grid scales within 5 ms. The trajectory of point P was roughly the same when th ber of grids was 164,491 and 255,303. In order to consider the computational efficie the convergence of the solver, the grid division method was adopted when the nu grids was 164,491. At this time, the mesh division conformed to the law that th quality is close to one; that is, the mesh quality was good.       The above conditions were combined for the simulation calculations, and finally, the cell positions stored in the program were derived by interpolation time and reasonable sparsity to obtain the flow trajectory of the erythrocyte model in this DLD array, as shown in Figure 7. In comparison with the experimental results, the erythrocyte model flowed in the same displacement mode of continuous tumbling.  The square pillar array was also established by imitating the control group in the experiment. According to the pillar shape in the experiment, the section side length of the The square pillar array was also established by imitating the control group in the experiment. According to the pillar shape in the experiment, the section side length of the square pillar was set at 15 µm, and the radius of 2.5 µm was rounded at the square corner. We set array layout parameters of G x = G y = 10 µm, λ = 1.25 µm, and an initial flow rate of u 0 = 1000 µm/s. As shown in Figure 8, The motion path of the erythrocyte model obtained after simulation calculation was also similar to that of the erythrocyte in the experiment. Therefore, the 2D linear elastic model is applicable to simulate the flow of erythrocyte in the DLD array.

2023, 8, x FOR PEER REVIEW
10 of 17 square pillar was set at 15 µm, and the radius of 2.5 µm was rounded at the square corner. We set array layout parameters of Gx = Gy = 10 µm, λ = 1.25 µm, and an initial flow rate of u0 = 1000 µm/s. As shown in Figure 8, The motion path of the erythrocyte model obtained after simulation calculation was also similar to that of the erythrocyte in the experiment. Therefore, the 2D linear elastic model is applicable to simulate the flow of erythrocyte in the DLD array.

Experimental Scheme Design
The pillar side lengths and the groove radius were decreased to form pillars B and C, taking advantage of the simulation model's benefit of being unconstrained by chip processing technology. Table 2 displays the measurements of the three pillar shapes: A, B, and

Experimental Scheme Design
The pillar side lengths and the groove radius were decreased to form pillars B and C, taking advantage of the simulation model's benefit of being unconstrained by chip processing technology. Table 2 displays the measurements of the three pillar shapes: A, B, and C. After the pillar shape was established, respective corresponding DLD arrays were established. Considering that the change in G y might make the period of arrays difficult to calculate, the downstream pillar distance G y of all arrays was set to be 10 µm, and the following scheme was designed and simulated: (1) Set the lateral pillar distance to G x = 10 µm, and set the flow period of the array as 20. According to the pillar side length H and the lateral pillar distance G y , the row shift distance λ of the three arrays was λa = 1.25 µm, λb = 1.1 µm, λc = 1 µm, respectively. Then the u 0 value was increased at a rate of 100 µm/s each time, with 1000 µm/s as the initial value. Multiple experiments were conducted until u 0 = 1500 µm/s; (2) Keep the array period constant and fix u 0 = 1000 µm/s. With G x as the variable and G x = 10 µm as the initial value, the lateral pillar distance G x of each array was continuously increased by 2 µm length each time for the simulation until G x = 30 µm; (3) Fix u 0 = 1000 µm/s, G x = 10 µm, shorten the array period to 10, and find the row shift distance λ after three array changes to λ A = 2.5 µm, λ B = 2.2 µm, and λ C = 2 µm, respectively, and perform the flow simulation.

Effects of Pillar Shape on Pressure Distribution within the Array
With the array arrangement parameters of G x = G y = 10 µm, u 0 = 1000 µm/s, and a period of 20, Figure 9 shows the fluid pressure distribution clouds within the three different I-pillars. The resistance to fluid flow generated by the array decreased with decreasing pillar size for the same array gap. In contrast, the fluid pressure variation in the array with Pillar C was smoother, which facilitated faster sample processing of the chip by increasing the fluid flux, and also allowed a longer array arrangement to improve the accuracy of chip sorting under the same pressure limitation. (3) Fix u0 = 1000 µm/s, Gx = 10 µm, shorten the array period to 10, and find the row shift distance λ after three array changes to λA = 2.5 µm, λB = 2.2 µm, and λC = 2 µm, respectively, and perform the flow simulation.

Effects of Pillar Shape on Pressure Distribution within the Array
With the array arrangement parameters of Gx = Gy = 10 µm, u0 = 1000 µm/s, and a period of 20, Figure 9 shows the fluid pressure distribution clouds within the three different I-pillars. The resistance to fluid flow generated by the array decreased with decreasing pillar size for the same array gap. In contrast, the fluid pressure variation in the array with Pillar C was smoother, which facilitated faster sample processing of the chip by increasing the fluid flux, and also allowed a longer array arrangement to improve the accuracy of chip sorting under the same pressure limitation.

Movement of Erythrocyte Models in Each Array
The number of array rows passed by erythrocytes in displacement mode, the number of times erythrocyte models rolled in the groove of I-pillars, and the number of times erythrocyte models rolled in the lateral channels were all recorded in each simulation. The three different types of data were contrasted to assess the stability of the erythrocyte model in the new array.

Movement of Erythrocyte Models in Each Array
The number of array rows passed by erythrocytes in displacement mode, the number of times erythrocyte models rolled in the groove of I-pillars, and the number of times erythrocyte models rolled in the lateral channels were all recorded in each simulation. The three different types of data were contrasted to assess the stability of the erythrocyte model in the new array.
The statistical situation in the first group of flow simulation is shown in Figure 9. The size of the I-pillar directly affects the movement of erythrocytes in the array. In the array formed by the larger Pillar A, the erythrocyte model flowed in displacement mode within the range of velocity variation and could almost maintain rolling every time it passed through the groove of the I-pillar. However, in the Pillar B array, the position of erythrocyte tumbling was randomly distributed in the gap between the depression of the pillar and the array, indicating that the ability of the I-pillar to induce erythrocyte tumbling decreased after the size reduction, but the Pillar B array could still keep the erythrocyte moving in displacement mode. As the pillar size continued to shrink, the number of erythrocyte rolls in the pillar depression became less in the Pillar C array, indicating that Pillar C almost lost its ability to induce erythrocyte rolling at this size. In addition, with the increase in suspension flow rate, the position of erythrocyte tumbling became more random, and the ability of I-pillar to induce erythrocyte tumbling decreased. However, the increased flow rate made it easier for the red cell model to flow in displacement mode.
(I) in Figure 10a shows the local flow of the erythrocyte model in the Pillar A array at u = 1500 µm/s. (I) in Figure 10b shows the local flow of the erythrocyte model in the Pillar B array when u 0 = 1100 µm/s. (I) in Figure 10c shows the local flow of the erythrocyte model in the Pillar C array at u 0 = 1000 µm/s. In all three cases, the tumbling position of the erythrocyte model was abnormal and did not roll steadily in the depression of the pillar.
x FOR PEER REVIEW 12 of 17 The statistical situation in the first group of flow simulation is shown in Figure 9. The size of the I-pillar directly affects the movement of erythrocytes in the array. In the array formed by the larger Pillar A, the erythrocyte model flowed in displacement mode within the range of velocity variation and could almost maintain rolling every time it passed through the groove of the I-pillar. However, in the Pillar B array, the position of erythrocyte tumbling was randomly distributed in the gap between the depression of the pillar and the array, indicating that the ability of the I-pillar to induce erythrocyte tumbling decreased after the size reduction, but the Pillar B array could still keep the erythrocyte moving in displacement mode. As the pillar size continued to shrink, the number of erythrocyte rolls in the pillar depression became less in the Pillar C array, indicating that Pillar C almost lost its ability to induce erythrocyte rolling at this size. In addition, with the increase in suspension flow rate, the position of erythrocyte tumbling became more random, and the ability of I-pillar to induce erythrocyte tumbling decreased. However, the increased flow rate made it easier for the red cell model to flow in displacement mode.
(I) in Figure 10a shows the local flow of the erythrocyte model in the Pillar A array at u = 1500 µm/s. (I) in Figure 10b shows the local flow of the erythrocyte model in the Pillar B array when u0 = 1100 µm/s. (I) in Figure 10c shows the local flow of the erythrocyte model in the Pillar C array at u0 = 1000 µm/s. In all three cases, the tumbling position of the erythrocyte model was abnormal and did not roll steadily in the depression of the pillar. The statistical situation in the second group of flow simulations is shown in Figure  11. With the increase in array transverse spacing Gx, the red cell model flowed more easily in the word mode. As in the first simulation, the larger Pillar A had more control over the movement of the erythrocytes and caused the erythrocyte model to roll each time it passed through the groove of the pillar. However, with the increase in Gx, the smaller Pillar B and Pillar C rapidly lost the ability to make the erythrocyte model move in displacement mode, and the erythrocyte roll was very random. The statistical situation in the second group of flow simulations is shown in Figure 11. With the increase in array transverse spacing G x , the red cell model flowed more easily in the word mode. As in the first simulation, the larger Pillar A had more control over the movement of the erythrocytes and caused the erythrocyte model to roll each time it passed through the groove of the pillar. However, with the increase in G x , the smaller Pillar B and Pillar C rapidly lost the ability to make the erythrocyte model move in displacement mode, and the erythrocyte roll was very random.
(I) in Figure 11a shows the local flow of the erythrocyte model in the Pillar A array at G x = 30 µm. (I) in Figure 11b shows the local flow of the erythrocyte model in the Pillar B array at G x = 18 µm. These two conditions show the movement of the red cell model after the array was unable to make the model flow in displacement mode as the value of G x increased.
When the array period was shortened and the λ value was increased, the erythrocyte model in the array composed of three pillars flowed in the zigzag mode more easily. As shown in Figure 12, the transfer of erythrocytes from the Pillar A array to another lateral channel occurred behind the second row. In the array of the other two pillars, the shift occurs behind the third row.
11. With the increase in array transverse spacing Gx, the red cell model flowed more easily in the word mode. As in the first simulation, the larger Pillar A had more control over the movement of the erythrocytes and caused the erythrocyte model to roll each time it passed through the groove of the pillar. However, with the increase in Gx, the smaller Pillar B and Pillar C rapidly lost the ability to make the erythrocyte model move in displacement mode, and the erythrocyte roll was very random. (I) in Figure 11a shows the local flow of the erythrocyte model in the Pillar A array at Gx = 30 µm. (I) in Figure 11b shows the local flow of the erythrocyte model in the Pillar B array at Gx = 18 µm. These two conditions show the movement of the red cell model after the array was unable to make the model flow in displacement mode as the value of Gx  Figure 13 shows the surface stress change in the erythrocyte model in the second group of experiments when it moved in the Pillar B array with Gx = 18 µm. Figure 13a shows the surface stress change in the erythrocytes in the rolling process. During this process, the surface stress of erythrocytes increased and reached its maximum value. During the subsequent movement, the erythrocyte model was continuously affected by fluid, and the stress curve was in the shape of wavy lines, indicating that the fluid's effect on the erythrocytes was constantly changing but relatively stable. Figure 13b shows the surface stress changes in the erythrocyte model when it entered the zigzag mode. In this process, the collision between the erythrocyte model and the pillar was weak, resulting in a small thrust of the erythrocyte model as it was returning to the flow channel, and the change in its surface stress was relatively more stable.

Conclusions
In comparison with existing experimental results, it was found that under the condition of a low degree of deformation, the linear elastic 2D red cell model based on the finite element method could be used to simulate and analyze the motion path of red cells in Ishaped DLD arrays with fairly accurate results. After a comprehensive analysis of the research results in this paper, the following conclusions are drawn: (1) The increase in row shift distance λ and lateral pillar distance Gx will affect the effect of I-pillar arrays, inducing erythrocytes to move in displacement mode. If the increase is too high, the array will lose the function of sorting and enriching erythrocytes.
(2) The decrease in pillar size and the increase in fluid velocity will make the effect of  Figure 13 shows the surface stress change in the erythrocyte model in the second group of experiments when it moved in the Pillar B array with G x = 18 µm. Figure 13a shows the surface stress change in the erythrocytes in the rolling process. During this process, the surface stress of erythrocytes increased and reached its maximum value. During the subsequent movement, the erythrocyte model was continuously affected by fluid, and the stress curve was in the shape of wavy lines, indicating that the fluid's effect on the erythrocytes was constantly changing but relatively stable. Figure 13b shows the surface stress changes in the erythrocyte model when it entered the zigzag mode. In this process, the collision between the erythrocyte model and the pillar was weak, resulting in a small thrust of the erythrocyte model as it was returning to the flow channel, and the change in its surface stress was relatively more stable.  Figure 13 shows the surface stress change in the erythrocyte model in the second group of experiments when it moved in the Pillar B array with Gx = 18 µm. Figure 13a shows the surface stress change in the erythrocytes in the rolling process. During this process, the surface stress of erythrocytes increased and reached its maximum value. During the subsequent movement, the erythrocyte model was continuously affected by fluid, and the stress curve was in the shape of wavy lines, indicating that the fluid's effect on the erythrocytes was constantly changing but relatively stable. Figure 13b shows the surface stress changes in the erythrocyte model when it entered the zigzag mode. In this process, the collision between the erythrocyte model and the pillar was weak, resulting in a small thrust of the erythrocyte model as it was returning to the flow channel, and the change in its surface stress was relatively more stable.

Conclusions
In comparison with existing experimental results, it was found that under the condition of a low degree of deformation, the linear elastic 2D red cell model based on the finite element method could be used to simulate and analyze the motion path of red cells in Ishaped DLD arrays with fairly accurate results. After a comprehensive analysis of the research results in this paper, the following conclusions are drawn: (1) The increase in row shift distance λ and lateral pillar distance Gx will affect the effect of I-pillar arrays, inducing erythrocytes to move in displacement mode. If the increase is too high, the array will lose the function of sorting and enriching erythrocytes.
(2) The decrease in pillar size and the increase in fluid velocity will make the effect of I-pillar-induced erythrocyte tumbling worse. When the pillar side length is reduced to 10

Conclusions
In comparison with existing experimental results, it was found that under the condition of a low degree of deformation, the linear elastic 2D red cell model based on the finite element method could be used to simulate and analyze the motion path of red cells in I-shaped DLD arrays with fairly accurate results. After a comprehensive analysis of the research results in this paper, the following conclusions are drawn: (1) The increase in row shift distance λ and lateral pillar distance G x will affect the effect of I-pillar arrays, inducing erythrocytes to move in displacement mode. If the increase is too high, the array will lose the function of sorting and enriching erythrocytes.
(2) The decrease in pillar size and the increase in fluid velocity will make the effect of I-pillar-induced erythrocyte tumbling worse. When the pillar side length is reduced to 10 µm and the semicircle diameter at the pillar groove is reduced to 3 µm, it is difficult for the I-pillar to induce erythrocytes to roll stably. However, if appropriate array arrangement parameters are selected, smaller I-pillars can be used for sorting and enrichment of erythrocytes.
(3) The collision degree of erythrocytes in contact with the pillar affects the movement mode of the following erythrocytes, so that the erythrocytes carry out stable collisions at the pillar, which can be used as a method to optimize the performance of the I-pillar DLD chip.
In this paper, three I-pillars were designed, and the corresponding arrays were established. Through analytical simulations, it can be seen that the performance of DLD arrays composed of smaller I-beam columns can be further improved by using methods such as adjusting the period and offset distance and increasing the flow rate of the suspension, which provides a new design direction for the optimization of I-beam DLD chips.