Three-Dimensional Numerical Simulation of Particle Focusing and Separation in Viscoelastic Fluids

Particle focusing and separation using viscoelastic microfluidic technology have attracted lots of attention in many applications. In this paper, a three-dimensional lattice Boltzmann method (LBM) coupled with the immersed boundary method (IBM) is employed to study the focusing and separation of particles in viscoelastic fluid. In this method, the viscoelastic fluid is simulated by the LBM with two sets of distribution functions and the fluid–particle interaction is calculated by the IBM. The performance of particle focusing under different microchannel aspect ratios (AR) is explored and the focusing equilibrium positions of the particles with various elasticity numbers and particle diameters are compared to illustrate the mechanism of particle focusing and separation in viscoelastic fluids. The results indicate that, for particle focusing in the square channel (AR = 1), the centerline single focusing becomes a bistable focusing at the centerline and corners as El increases. In the rectangular channels (AR < 1), particles with different diameters have different equilibrium positions. The equilibrium position of large particles is closer to the wall, and large particles have a faster lateral migration speed and few large particles migrate towards the channel center. Compared with the square channel, the rectangular channel is a better design for particle separation.


Introduction
Particle focusing is essential for particle counting and detection [1]. In addition, particle separation from mixed samples is a key step in applications such as medical diagnosis and chemical analysis [2]. In recent years, microfluidic technology has become an important means of particle focusing and separation, because of the advantages of small sample volume, high throughput, and simple control [3][4][5][6]. Among them, the viscoelastic microfluidics using viscoelastic fluids can achieve a single equilibrium position focusing of particles in a simple straight channel, and particles in viscoelastic fluids can also migrate to the walls or corners by overcoming the repulsion between the wall and the particles due to inertia. These interesting phenomena can be utilized and applied to particle focusing and separation through parameter control. A variety of experiments have been conducted for the further study of particle focusing and separation in viscoelastic fluids. The concept of elasto-inertial particle focusing was first proposed by Yang et al. [7], who realized the single-line focusing of particles in the center of the square channel by combining the elastic force and inertial lift. The two most important parameters for the migration of particles in the viscoelastic fluids are the Reynolds number Re (Re = ρUL/η t , ρ is the fluid density; U is the characteristic velocity; L = 2hw/(h + w) is the characteristic length, h and w are the height and width of channel, respectively; and η t is the total viscosity of the viscoelastic fluid) and Weissenberg number Wi (Wi = λ p U/L, λ p is the relaxation time of the polymer solution), and they represent the inertial and elastic effect, respectively. The ratio of the two parameters

3D Lattice Boltzmann Method (LBM) of Viscoelastic Fluid
In the simulations, the dimensionless Navier-Stokes equations of incompressible flow are given by [26]: ρ[ ∂u ∂t + (u·∇)u] = −∇p + η s ∇ 2 u + ∇·τ + F D , (2) where ρ, u, p, η s , τ, and F D are the density, velocity, pressure, Newtonian solvent viscosity, viscoelastic stress tensor, and external force, respectively. Viscoelastic solvent can be regarded as a mixture of polymer and Newtonian solvent. The total viscosity of viscoelastic solvent η t is the sum of polymer viscosity η p and η s , and the parameter β = η s /η t is used to describe the viscosity ratio. The external force F D is imposed to drive the flow. The Oldroyd-B constitutive equation [22] that describes viscoelastic fluids is given by: where Dτ Dt = ∂τ ∂t + (u·∇)τ is the material derivative, λ P is the polymer relaxation time, ν is the polymer diffusion parameter, d = (∇u + (∇u) T )/2 is the rate of strains.
The viscoelastic flow field in microchannel is calculated by three-dimensional nineteen velocity (D3Q19) single relaxation time lattice Boltzmann Bhatnagar-Gross-Krook (LBGK) model [27], as follows: where α is the 19 discrete directions in the D3Q19 model, f α (x, t) is the velocity distribution function at node x and time t, λ is the relaxation time, and c α is the discrete lattice velocity. The 19 discrete lattice velocity can be given by: where c= ∆x/∆t, ∆x and ∆t are the lattice spacing and the time step. The f eq α (ρ, u) in Equation (4) is equilibrium distribution function and can be calculated from the macroscopic fluid density ρ and velocity u as: (c α ·u) 2 c 4 − 1.5 where weight coefficients ω 0 =1/3, ω 1-6 =1/18, ω 7-18 =1/36. The F α (x, t) in Equation (4) is external force term [28] and can be calculated as where F = F D + F A + F V is the total external force acting on the fluid. F D , F A , and F V are the driving force, the acting force from particles, and the elastic force, respectively. The macroscopic fluid density ρ, velocity u, and kinematic viscosity υ of the fluid can be given as: The elastic force F V = ∇·τ is the calculation result from the stress tensor τ of the viscoelastic flow. The stress tensor τ can be computed from the tensor distribution function G ijα (x, t) [22] which is on the same lattice node of f α (x, t). The LBGK model of G ijα (x, t) is also based on D3Q19 as: where G ijα (x, t) is the distribution function at node x and time t, ∆ t, and c α are the time step and discrete lattice velocity of G ijα (x, t), λ is the relaxation parameter, G eq ijα (x, t) is the equilibrium tensor distribution function and can be calculated as: where c = ∆x/∆ t, τ ij is the component of τ.
The source term χ ijα in Equation (9) can be given as: where χ ij is the component of tensor χ which can be given as: The stress tensor τ of the viscoelastic flow can be obtained through: and then the elastic force F V can be calculated. In this study, periodic boundary condition is applied to simulate an infinite length channel and non-slip boundary condition is applied for walls. The distribution function f α (x, t) and G ijα (x, t) on the walls are calculated by non-equilibrium extrapolation method [29], as follows: where u w = 0 is the velocity of walls, ρ f and τ ijf are the density and stress tensor component of adjacent ijα (x f , t) are the distribution functions and equilibrium distribution functions of adjacent fluid nodes.

IBM for the Interaction between Particles and Fluid
In the simulations, particle can be built through an 3D finite element membrane model which divides the surface of the particle into N triangular surface elements. The vertices of all triangular elements constitute the nodes of the particle surface boundary. The information of particle surface nodes and fluid nodes can be transferred to each other through IBM [30] to realize the interaction between particles and fluid. The total energy E of the particle deformation produced by fluid shear is calculated from E = E S + E B + E A + E V , where E S , E B , E A , and E V denote strain energy, bending energy, area energy, and volume energy, respectively, and can be calculated as: Micromachines 2020, 11, 908 5 of 17 where I 1 = λ 2 1 + λ 2 2 − 2 and I 2 = λ 2 1 λ 2 2 − 1 are strain invariants (λ 1 , λ 2 are eigenvalues). κ S , κ α , κ B , κ A , and κ V are the surface elastic shear modulus, area dilation modulus, binding modulus, surface modulus, and volume modulus, respectively. When the elastic moduli are large enough, the simulated particle can be regarded as a rigid particle. φ n and φ n,0 , A and A 0 , V and V 0 are the angle between normal vectors of adjacent triangular faces, particle surface area, particle volume after, and before the particle deformation, respectively.
Then the restoring force on particle nodes can be calculated from f (x n ) = −∂E/∂x n , where x n is the coordinate of the particle node. Similarly, the coordinate of the fluid node is expressed by x f . The acting force F A from particle on x f can be calculated by IBM as: where At this time, F D , F A , and F V in the force term F have all been obtained. The velocity distribution u f x f of the flow field can be obtained by calculating Equations (4)- (8). The velocity distribution u p (x n ) of the particle node can be updated by , which completes the interaction between particle and fluid for one time step.

Channel Model and Validation
The schematic diagram of the particles in the straight channel with viscoelastic fluid is shown in Figure 1. For a rectangular channel, the height and width are h = 40 µm, w = h/AR, and for the square channel, w = h = 50 µm. The trajectories of particles in the viscoelastic fluid can be simulated by the method given above. In the simulation, the fluid can be driven by driving force F D . Thus, the acceleration a of the fluid can be obtained from the Poiseuille equation as a = 32υ 2 Re/L 3 . The channel can be simulated with infinite length in the x-axis direction benefiting from the periodic boundary conditions of the inlet and outlet. The number of membranes N of 6 µm and 12 µm diameter particles in the simulation are 120 and 480, respectively. The elastic moduli for a rigid particle are set as κ S = 3.2 × 10 −1 N/m, κ α = 3.2 × 10 −1 N/m, κ B = 3.2 × 10 −13 Nm, κ A = 3.2 × 10 −2 N/m, and κ V = 3.2 × 10 4 N/m 2 . The other parameters in simulation can be given as ∆x = 1.0 × 10 −6 m, ∆t = ∆ t = 1.0 × 10 −7 s, υ = 1.0 × 10 −6 m 2 /s, To validate the numerical simulation model, the particle-free viscoelastic fluid evolution and the particle rotation in viscoelastic shear fluid are conducted respectively. For the particle-free fluid model, a microtube with a circle cross-section is employed and the simulation parameters are set as Re = 1.0, Wi= 0.3, and β = 0.3. The numerical and analytical results of the velocity distribution, normal-stress components τ xx , and extra-stress components τ xy at the channel center plane are shown in Figure 2. The simulation results are in good agreement with the analytical solutions. The components of stress tensor τ can be calculated as: To validate the numerical simulation model, the particle-free viscoelastic fluid evolution and the particle rotation in viscoelastic shear fluid are conducted respectively. For the particle-free fluid model, a microtube with a circle cross-section is employed and the simulation parameters are set as Re = 1.0, Wi= 0.3, and β = 0.3. The numerical and analytical results of the velocity distribution, normalstress components τxx, and extra-stress components τxy at the channel center plane are shown in Figure 2. The simulation results are in good agreement with the analytical solutions. The components of stress tensor τ can be calculated as: For fluid-particle interaction, the simulation of a spherical particle rotation in an Oldroyd-B viscoelastic shear flow is employed to validate the present method. The shear flow model is shown in Figure 3a. In the shear flow model, the velocities on the upper and lower board nodes are set as the same constant values with opposite directions. Periodic boundary conditions are applied around the fluid field, and the particle is placed at the fluid filed center. The angular velocities of the particle are calculated under different Wi, as shown in Figure 3b. The angular velocity of particle decreases monotonously as Wi increases, and the present results agree well with the previous results.  For fluid-particle interaction, the simulation of a spherical particle rotation in an Oldroyd-B viscoelastic shear flow is employed to validate the present method. The shear flow model is shown in Figure 3a. In the shear flow model, the velocities on the upper and lower board nodes are set as the same constant values with opposite directions. Periodic boundary conditions are applied around the fluid field, and the particle is placed at the fluid filed center. The angular velocities of the particle are calculated under different Wi, as shown in Figure 3b. The angular velocity of particle decreases monotonously as Wi increases, and the present results agree well with the previous results.

Results and Discussion
The shape of cross-section, velocity distribution, and stress tensor distribution on the cross-section are all symmetrical. Therefore, the square cross-section can be divided into 8 sections according to its four symmetry axes (the rectangular cross-section is divided into 4 sections), as shown in Figure 4. The trajectories of the particle on the entire cross-section can be obtained by only simulating the migration of particle on Section 1.

Results and Discussion
The shape of cross-section, velocity distribution, and stress tensor distribution on the crosssection are all symmetrical. Therefore, the square cross-section can be divided into 8 sections according to its four symmetry axes (the rectangular cross-section is divided into 4 sections), as shown in Figure 4. The trajectories of the particle on the entire cross-section can be obtained by only simulating the migration of particle on Section 1.

Particle Focusing
In the square channel, the particle with a diameter of 6 μm (d/h = 0.12) is released at different initial positions on Section 1 under different elasticity numbers. The Reynolds number is Re = 1.0 and the Weissenberg number increases from 0.05 to 0.6, which makes the El range from 0.05 to 0.6. The migration trajectories of particle at different initial positions of the cross-section are shown in Figure   Figure 4. Schematic diagram of the selected simulated section on the cross-section.

Particle Focusing
In the square channel, the particle with a diameter of 6 µm (d/h = 0.12) is released at different initial positions on Section 1 under different elasticity numbers. The Reynolds number is Re = 1.0 and the Weissenberg number increases from 0.05 to 0.6, which makes the El range from 0.05 to 0.6. The migration trajectories of particle at different initial positions of the cross-section are shown in Figure 5. In the particle trajectory figures, the dots denote the initial positions of the particles, the lines of different colors correspond to the migration trajectories of the particles at different initial positions of the cross-section, and the arrow is the moving direction of the particle on the cross-section. For El = 0.05, particles at different initial positions all migrate towards the centerline of the channel and the particles migrate slowly due to the small elastic force. As El increases to 0.1, the efficiency of particle focusing increases, and the particle close to the corner migrates towards the corner instead of the centerline of the channel. What is more, there seems to be a separatrix in the channel, the particle within the separatrix migrates towards the channel centerline and the particle outside the separatrix migrates towards the wall or corner. This phenomenon is consistent with previous experimental [8,34,35] and numerical [20,36] results. For El = 0.3, the separatrix shrinks, and the particles close to the wall migrate towards the corner. As El further increases to 0.6, this phenomenon becomes more obvious, and the separatrix shrinks further. These phenomena indicate that there always has a centerline focusing equilibrium position under different El. For small El, the particle focusing efficiency is slow, but for large El, some particles migrate to corners and single-line focusing cannot be achieved. Therefore, El is the key parameter to obtain different focusing patterns. In addition, the trajectories of the particles at the symmetry axes are relatively straight due to the symmetry of the fluid field, while the particles on the other positions deflect towards the diagonal during the migration process. Yu et al. [20] also discovered this phenomenon, and they realized the particle focusing on the diagonal by increasing Re. In their research, the diagonal equilibrium position exists under the conditions of (Re = 10, Wi = 0.05, El = 0.005), (Re = 50, Wi = 0.25, El= 0.005), and (Re= 100, Wi = 0.5, El = 0.005), but as El increases, the particles will migrate to the center of the channel.
For the rectangular channel, the channel with AR = 1/2 is studied firstly. Similarly, the performance of particle focusing under different El is compared. The Reynolds number is Re = 2.5 and the Weissenberg number increases from 0.025 to 1.2, which makes the El range from 0.01 to 0.48. The migration trajectories of particle at different initial positions of the cross-section are shown in Figure 6a-d. For El = 0.01, the particles migrate very slowly, and the elastic forces acting on the particles are very small and not enough to make the particles fully focused to the channel center. For El = 0.08, the efficiency of particle focusing increases, and the particles close to the corner migrate towards the corners, which is consistent with the square channel. However, the longitudinal migrations of particles are significantly faster than the lateral migrations, and the particles are firstly focused at the long midline plane (z = 0) and then migrate towards the equilibrium position at the centerline. For the particles at the long midline plane, due to the smaller inertial lift and elastic force, particles close to the channel center migrate more slowly. Meanwhile, particles close to the channel center also have slower rotation speed as shown in Figure 6e. As El increases to 0.24, particle focusing becomes more interesting, the off-centerline focusing equilibrium position appears, which means that the particles can be focused into double lines or even multiple lines. Compared with El = 0.08, since the large El may attenuate the particle rotation, and the particle rotation speed is further reduced as shown in Figure 6f, which weakens the fluid stretching [37] and makes the elastic force directing to the channel center smaller, so that the particle can be focused at the off-centerline equilibrium position. As Wi further increases to 0.48, the equilibrium positions of the particles change again, almost all the particles migrate towards the walls and corners.
trajectories of the particles at the symmetry axes are relatively straight due to the symmetry of the fluid field, while the particles on the other positions deflect towards the diagonal during the migration process. Yu et al. [20] also discovered this phenomenon, and they realized the particle focusing on the diagonal by increasing Re. In their research, the diagonal equilibrium position exists under the conditions of (Re = 10, Wi = 0.05, El = 0.005), (Re = 50, Wi = 0.25, El= 0.005), and (Re= 100, Wi = 0.5, El = 0.005), but as El increases, the particles will migrate to the center of the channel. For the rectangular channel, the channel with AR = 1/2 is studied firstly. Similarly, the performance of particle focusing under different El is compared. The Reynolds number is Re = 2.5 and the Weissenberg number increases from 0.025 to 1.2, which makes the El range from 0.01 to 0.48. The migration trajectories of particle at different initial positions of the cross-section are shown in Figure 6a-d. For El = 0.01, the particles migrate very slowly, and the elastic forces acting on the particles are very small and not enough to make the particles fully focused to the channel center. For El = 0.08, the efficiency of particle focusing increases, and the particles close to the corner migrate Channels with lower aspect ratios are also studied. Figure 7 shows the migration trajectories of particles with the diameter of 6 µm (d/h=0.15) in the channels with the aspect ratio of 1/3 and 1/4 under Re = 2.5 and Wi = 0.2. Compared with Figure 6b, with the aspect ratio decreasing, the focusing equilibrium position is farther away from the centerline, and the lateral migration speed of particles at the long midline plane further slows down. Within the same time (1,200,000 dimensionless time units), particles in the channel with the aspect ratio of 1/2 can be focused near the equilibrium position, while particles in the channel with the aspect ratio of 1/4 only migrate to the long midline plane. We can predict that when the aspect ratio is further reduced, the phenomenon that particles are focused at the entire long midline plane will be formed. Seo et al. [11] confirmed this prediction and realized the focusing of the long midline plane of the red blood cells (RBCs) by using a low aspect ratio channel with the height of 50 µm and the width of 500 µm, and used a holographic microscope to perform quantitative phase imaging on the midline plane to achieve the monitoring and counting of RBCs. This plane focusing can make the height of particles is highly uniform, which reduces the out-of-focus blurring and improves the detection sensitivity. In addition, for particles whose initial positions are close to the center of the channel, it is difficult to migrate laterally after they migrate to the long midline plane. Therefore, particles are actually focused at the region between the two equilibrium positions on the long midline plane, and the smaller the aspect ratio, the larger the region is. towards the corners, which is consistent with the square channel. However, the longitudinal migrations of particles are significantly faster than the lateral migrations, and the particles are firstly focused at the long midline plane (z = 0) and then migrate towards the equilibrium position at the centerline. For the particles at the long midline plane, due to the smaller inertial lift and elastic force, particles close to the channel center migrate more slowly. Meanwhile, particles close to the channel center also have slower rotation speed as shown in Figure 6e. As El increases to 0.24, particle focusing becomes more interesting, the off-centerline focusing equilibrium position appears, which means that the particles can be focused into double lines or even multiple lines. Compared with El = 0.08, since the large El may attenuate the particle rotation, and the particle rotation speed is further reduced as shown in Figure 6f, which weakens the fluid stretching [37] and makes the elastic force directing to the channel center smaller, so that the particle can be focused at the off-centerline equilibrium position. As Wi further increases to 0.48, the equilibrium positions of the particles change again, almost all the particles migrate towards the walls and corners. Channels with lower aspect ratios are also studied. Figure 7 shows the migration trajectories of particles with the diameter of 6 μm (d/h=0.15) in the channels with the aspect ratio of 1/3 and 1/4 under Re = 2.5 and Wi = 0.2. Compared with Figure 6b, with the aspect ratio decreasing, the focusing equilibrium position is farther away from the centerline, and the lateral migration speed of particles at the long midline plane further slows down. Within the same time (1,200,000 dimensionless time units), particles in the channel with the aspect ratio of 1/2 can be focused near the equilibrium position, while particles in the channel with the aspect ratio of 1/4 only migrate to the long midline plane. We can predict that when the aspect ratio is further reduced, the phenomenon that particles are focused at the entire long midline plane will be formed. Seo et al. [11] confirmed this prediction and realized the focusing of the long midline plane of the red blood cells (RBCs) by using a low aspect ratio channel with the height of 50 μm and the width of 500 μm, and used a holographic microscope to perform quantitative phase imaging on the midline plane to achieve the monitoring and counting of RBCs. This plane focusing can make the height of particles is highly uniform, which reduces the out-offocus blurring and improves the detection sensitivity. In addition, for particles whose initial positions are close to the center of the channel, it is difficult to migrate laterally after they migrate to the long midline plane. Therefore, particles are actually focused at the region between the two equilibrium positions on the long midline plane, and the smaller the aspect ratio, the larger the region is. To further explore the particle migration in the channels with AR < 1, the stress distribution and particle rotation on the cross-section are studied. In the rectangular channel, the velocity distribution and stress distribution on the long and short symmetry axis of the cross-section are different. The first normal stress difference N1 on the long symmetry axis is significantly smaller than that on the short symmetry axis, as shown in Figure 8a, and there is a larger "zero N1" region (N1 is close to zero in this region) along the long symmetry axis, and this region increases as the aspect ratio decreases. Since the elastic force acting on the particles is proportional to N1 [38], the elastic force acting on the particle in the "zero N1" region is small and close to zero. Moreover, for particles near the center on the long symmetry axis, especially for AR = 1/4, the velocity difference around the particles is also very small, as shown in Figure 8b. It means that the shear-induced inertial lift and wall-induced inertial lift (the particles are far from the wall) acting on the particles are very small. For particle rotation, the rotation process of the particle is described by three angles θx, θy, and θz, which are the angles between the vector P and the x-axis, y-axis, and z-axis, respectively, as shown in Figure 8c. The vector P points from the particle center to a fixed node on the surface of the particle. The initial vector P is set along the x-axis, and the changes of θx, θy, and θz with particle migration are shown in Figure 8d,e. The θy remains unchanged at the beginning, which indicates that the particle rotates around the y-axis. After the particle is focused at the long midline plane, for (0, −0.3), the particle hardly rotates, and for (−0.4, −0.3), the particle turns to rotate around the z-axis, but the To further explore the particle migration in the channels with AR < 1, the stress distribution and particle rotation on the cross-section are studied. In the rectangular channel, the velocity distribution and stress distribution on the long and short symmetry axis of the cross-section are different. The first normal stress difference N 1 on the long symmetry axis is significantly smaller than that on the short symmetry axis, as shown in Figure 8a, and there is a larger "zero N 1 " region (N 1 is close to zero in this region) along the long symmetry axis, and this region increases as the aspect ratio decreases. Since the elastic force acting on the particles is proportional to N 1 [38], the elastic force acting on the particle in the "zero N 1 " region is small and close to zero. Moreover, for particles near the center on the long symmetry axis, especially for AR = 1/4, the velocity difference around the particles is also very small, as shown in Figure 8b. It means that the shear-induced inertial lift and wall-induced inertial lift (the particles are far from the wall) acting on the particles are very small.   For particle rotation, the rotation process of the particle is described by three angles θ x , θ y , and θ z , which are the angles between the vector P and the x-axis, y-axis, and z-axis, respectively, as shown in Figure 8c. The vector P points from the particle center to a fixed node on the surface of the particle. The initial vector P is set along the x-axis, and the changes of θ x , θ y , and θ z with particle migration are shown in Figure 8d,e. The θ y remains unchanged at the beginning, which indicates that the particle rotates around the y-axis. After the particle is focused at the long midline plane, for (0, −0.3), the particle hardly rotates, and for (−0.4, −0.3), the particle turns to rotate around the z-axis, but the rotation speed is very slow. This proves that the velocity shear around particles is small. Figure 8f more clearly shows the migration trajectories of the two particles on (0.4, −0.3) and (0, −0.3). The particles migrate to the long midline plane after the time t1 (t1 = 300,000 dimensionless time units), but hardly migrate laterally during the longer time t2 (t2 = 900,000 dimensionless time units). Therefore, in the "zero N 1 " region, both the elastic and inertial effects acting on particles are small. The particles are relatively stable and difficult to migrate laterally. The above results are obtained under Re = 2.5, Wi = 0.2, El = 0.08.

Particle Separation
In the viscoelastic fluid of straight microchannel, equilibrium position of particles may differ with various diameters, which can be utilized to particle separation [15,16]. To study the different size particle performance in various AR microchannels, the trajectories of particle with diameters of 6 and 12 µm are calculated and systematically analyzed.
In the square channel, the migration trajectories of particles with diameters of 6 and 12 µm (d/h = 0.12 and 0.24) under Re = 1.0 and Wi = 0.3 are shown in Figure 9a,b. The small particle has a larger region within the separatrix and more of them can be focused to the centerline of the channel. The large particle has a faster lateral migration speed as shown in Figure 9c,d. It is difficult to realize particle separation in the square channel due to the same focusing equilibrium positions (the center and corner of the channel) of the large and small particles. Even if the initial position of the mixed particles can be set near the wall and corner with the help of sheath flow, so that some small particles can migrate to the center of the channel while large particles migrate to the corner, the separation efficiency is extremely low. It is why few researchers use square channel for particle separation. However, the characteristic that large particles have a faster lateral migration speed to the centerline can be utilized in the sudden expansion channel. Nam et al. [39,40] proposed a two-stage flow channel structure for particle separation. In the first stage, particles with different sizes are focused to the center of the channel, and then in the second stage, large particles can be separated from small particles in the sudden expansion channel due to the faster migration. In addition, large particles are easier to deflect towards diagonal line during the migration, which may be caused by the larger velocity shear around the large particles.
In the rectangular channel, particles of various diameters have different focusing equilibrium positions for separation and the problems in the square channel disappear. Figure 10 shows the migration trajectories of particles with d/h = 0.15 and 0.30 in the channel with the aspect ratio of 1/2 under Re = 2.5, Wi = 0.2, and El = 0.08. The focusing equilibrium position of small particles is at the center of the channel, while large particles migrate towards off-center equilibrium positions. This is consistent with the results of Liu et al. [15] and Nam et al. [16]. The channel with smaller aspect ratios is also studied. Figure 11 shows the migration trajectories of particles with the diameters of d/h = 0.15 and 0.30 in the channel with the aspect ratio of 1/3 under different El. For El = 0.08, the migration performance of large particles is similar to small ones, but the focusing equilibrium position of large particles is farther away from the centerline of the channel than small particles. As El increases to 0.24, the elastic effect acting on the particles is enhanced. For small particles, the equilibrium position is farther away from the centerline of the channel, while for large particles, almost all of them migrate towards the walls and corners with faster lateral migration speed. Comparing the equilibrium positions of different size particles in El = 0.08 and 0.24 rectangular channels, El has a greater influence on the migration of large particles. When El = 0.24, the large particle can be fully focused near the channel walls and the separation efficiency for different size particles is higher than that of the El = 0.08 channel.
can migrate to the center of the channel while large particles migrate to the corner, the separation efficiency is extremely low. It is why few researchers use square channel for particle separation. However, the characteristic that large particles have a faster lateral migration speed to the centerline can be utilized in the sudden expansion channel. Nam et al. [39,40] proposed a two-stage flow channel structure for particle separation. In the first stage, particles with different sizes are focused to the center of the channel, and then in the second stage, large particles can be separated from small particles in the sudden expansion channel due to the faster migration. In addition, large particles are easier to deflect towards diagonal line during the migration, which may be caused by the larger velocity shear around the large particles.  In the rectangular channel, particles of various diameters have different focusing equilibrium positions for separation and the problems in the square channel disappear. Figure 10 shows the migration trajectories of particles with d/h = 0.15 and 0.30 in the channel with the aspect ratio of 1/2 under Re = 2.5, Wi = 0.2, and El = 0.08. The focusing equilibrium position of small particles is at the center of the channel, while large particles migrate towards off-center equilibrium positions. This is consistent with the results of Liu et al. [15] and Nam et al. [16]. The channel with smaller aspect ratios is also studied. Figure 11 shows the migration trajectories of particles with the diameters of d/h = 0.15 and 0.30 in the channel with the aspect ratio of 1/3 under different El. For El = 0.08, the migration performance of large particles is similar to small ones, but the focusing equilibrium position of large particles is farther away from the centerline of the channel than small particles. As El increases to 0.24, the elastic effect acting on the particles is enhanced. For small particles, the equilibrium position is farther away from the centerline of the channel, while for large particles, almost all of them migrate towards the walls and corners with faster lateral migration speed. Comparing the equilibrium positions of different size particles in El = 0.08 and 0.24 rectangular channels, El has a greater influence on the migration of large particles. When El = 0.24, the large particle can be fully focused near the channel walls and the separation efficiency for different size particles is higher than that of the El = 0.08 channel.
Overall, rectangular channel is a better design for particle separation than the square channel. What is more, the 1/2 AR rectangular channel can realize center focusing for small particles, while the 1/3 AR channel is able to focus the large particles fully near the sidewall. In addition, both of them can be applied to particle separation device according to the specific requirement.

Conclusions
In summary, we systematically explore the elasto-inertial focusing and separation of particles in viscoelastic fluids with different aspect ratios through the numerical simulation. For particle focusing, a single line focusing at the channel center can be achieved in a square channel. In the rectangular channel, particle focusing becomes more complicated and interesting. With the aspect ratio decreasing, the double-line focusing, multi-line focusing, and even plane focusing of particles will appear. The equilibrium position of particle also can be controlled by parameter El. More particles can be focused at the corners under a larger El in the square channel, while in the rectangular channel, the equilibrium positions of particles can be moved closer to the side walls with El increasing. Overall, rectangular channel is a better design for particle separation than the square channel. What is more, the 1/2 AR rectangular channel can realize center focusing for small particles, while the 1/3 AR channel is able to focus the large particles fully near the sidewall. In addition, both of them can be applied to particle separation device according to the specific requirement.

Conclusions
In summary, we systematically explore the elasto-inertial focusing and separation of particles in viscoelastic fluids with different aspect ratios through the numerical simulation. For particle focusing, a single line focusing at the channel center can be achieved in a square channel. In the rectangular channel, particle focusing becomes more complicated and interesting. With the aspect ratio decreasing, the double-line focusing, multi-line focusing, and even plane focusing of particles will appear. The equilibrium position of particle also can be controlled by parameter El. More particles can be focused at the corners under a larger El in the square channel, while in the rectangular channel, the equilibrium positions of particles can be moved closer to the side walls with El increasing.
For particle separation, large particles have a faster lateral migration speed and a smaller region within the separatrix, but cannot be effectively separated from the small particles in the square channel, because of their same focusing equilibrium position at the centerline. In the rectangular channel, large and small particles have different focusing equilibrium positions, and the large particles with equilibrium positions closer to the wall can be separated from small ones. Moreover, El has a greater influence on the migration of large particles, and higher El can help large particles to be focused closer to the wall, while small particles stay near the center of the channel at the long midline plane. We expect the results will be helpful for the design of microfluidic chips for particle/cell focusing and separation.