Computational inertial microfluidics: a review.

Since the discovery of inertial focusing in 1961, numerous theories have been put forward to explain the migration of particles in inertial flows, but a complete understanding is still lacking. Recently, computational approaches have been utilized to obtain better insights into the underlying physics. In particular, fundamental aspects of particle focusing inside straight and curved microchannels have been explored in detail to determine the dependence of focusing behavior on particle size, channel shape, and flow Reynolds number. In this review, we differentiate between the models developed for inertial particle motion on the basis of whether they are semi-analytical, Navier-Stokes-based, or built on the lattice Boltzmann method. This review provides a blueprint for the consideration of numerical solutions for modeling of inertial particle motion, whether deformable or rigid, spherical or non-spherical, and whether suspended in Newtonian or non-Newtonian fluids. In each section, we provide the general equations used to solve particle motion, followed by a tutorial appendix and specified sections to engage the reader with details of the numerical studies. Finally, we address the challenges ahead in the modeling of inertial particle microfluidics for future investigators.


Introduction
Microfluidics, a technology characterized by the engineered manipulation of fluids at the microscale, has shown considerable promise in point-of-care diagnostics and clinical studies [1]. Since its birth in the 1990s, this technology has been matured into a complex field impacting many commercial applications [2]. Isolation, fractionation, and purification of cells using microfluidic platforms have been a flourishing area of development in recent years with many successful translations into commercial products. Several review articles are available on particle/cell separation using microfluidics systems, with some of them being presented in [3][4][5].
Microfluidic systems normally leverage on the disparities in the intrinsic properties of different particle/cell populations (i.e., size, deformability, surface charge, and density) to achieve separations. These systems can be broadly classified as active and passive separation techniques. Active techniques rely on external force fields such as acoustic or magnetic for operation, while passive techniques (e.g., pinched flow fractionation, deterministic lateral displacement (DLD), and inertial microfluidics) rely only on the channel geometry and inherent hydrodynamic forces for functionality [6]. Among all existing microfluidic systems, inertial microfluidics, which takes the advantage of size-dependent hydrodynamic effects in microchannels, has become a promising approach for particle and cell separation due to its capacity in high-volume and high-throughput sample processing [7][8][9][10][11]. Inertial microfluidics has been employed for various applications, influencing a wide range of industries such as microbiology, biochemistry, and biotechnology, most of which can be found in [9,10,12]. The most prevalent structures used in inertial microfluidics are demonstrated in Fig. 1A.
Inertial microfluidics is defined as the migration of randomly dispersed particles toward specific equilibrium positions inside a microchannel. This phenomenon was first reported by Segre and Silberberg in 1961 [13]. Later, with the advent of microfluidics, it was extensively explored, resulting in numerous articles evaluating the underlying physics of this phenomenon.
To date, numerous research groups have attempted to provide a numerical solution to better understand fundamental aspects of particle focusing inside inertial microfluidic systems. These numerical solutions can be divided into three distinct classifications. First, asymptotic analysis by simplification of fluid equations predicts inertial forces acting on a particle within a channel.
It has been reported that the phenomena pictured by this solution are far from the real scenario, due to its oversimplifying assumptions and limitations. In the second classification, scientists utilize Navier-Stokes-based solutions for inertial microfluidics, where all prevalent issues of particle size or particle effects on fluid streamline in asymptotic solutions are addressed.
However, tracking the solid-fluid interface or the demanding calculation time for solving equations are challenges for Navier-Stokes-based approaches. As an alternative, researchers use the lattice Boltzmann method (LBM) for inertial flow modeling. The relative simplicity of algorithm parallelization, adding new physics, and the strength of LBM in the intermediate Re regime explain its rapidly increasing use for inertial microfluidics over the last decade.
Although numerous review papers exist for inertial microfluidics, there is a lack of a comprehensive review of computational inertial microfluidics. Accordingly, the primary aim of this review paper is to provide researchers with the latest updates regarding computational solutions of inertial particle motion within a microfluidic device. Here, we have reviewed all computational attempts for inertial particle motion, covering asymptotic calculations, Navier-Stokes-based approaches, and LBM.

Asymptotic solutions
In 1961, Segre and Silberberg demonstrated that a suspension of neutrally buoyant particles flowing in a tube with radius R migrate to positions 0.6 R from the centerline when inertial The asymptotic solutions are among the early attempts made for predicting particle migration in confined and non-confined flows. Despite their weakness to capture specific details, such as the particle-fluid interaction and disturbed fluid velocity profiles around the particle, and to apply to a complex 3D domain, these solutions are convincing enough to explain the phenomenon in the most time-efficient manner.
In 1961, Brenner investigated the motion of a particle toward a semi-infinite plane surface [16].
The effect of the plane surface is considered using two distinct cases: a rigid plane surface and a free plane surface. They found that resistance against a particle path increases in a wallbounded case (e.g., particle motion within a confined channel) compared to an unbounded case (e.g., particle movement in an infinite fluid domain) when a particle moves under otherwise identical conditions. Nowadays, we know that a particle experiences more resistance along its way near the walls. The results obtained for a solid surface or a free plane surface resulted in a drag force ( ) on the sphere which is shown by Eqs. (1) and (2) ] } (2) where U is the fluid velocity, and is a term that relates the ratio of sphere radius (b) to the distance of its surface from the plane (h), as shown by Eq. (3).
In another work, the force exerted on a spinning particle moving in a viscous fluid was investigated using a combination of Stokes and Oseen expansions [17]. Since both expansions are for the same solution and each is valid over its own range, the idea was that these two expansions should be matched at the interface. Through the mathematical derivations, lift and drag forces over the particle were extracted as by Eqs. (4) and (5), respectively.
where is the fluid density, is the angular velocity of the particle, is the particle diameter, Ω and V is the relative velocity between fluid and particle. These two formulas are acceptable for a moving particle in an initially quiescent flow. The transverse force on a particle in a 3D Poiseuille flow was calculated based on Eq. (6).
where is the radius of the tube. Equation (6) shows a force in the direction of the central 0 axis where the force vanishes; however, since the Segre and Silberberg observations [13,18] show that particles tend to concentrate on an annulus of a tube, this formula requires further development and modification.
These studies were further continued by Saffman [19] who showed that if a spherical particle travels through a highly viscous liquid with the velocity ( ) relative to a uniform and simple shear flow, it will encounter a lift force perpendicular to both rotational and translational velocity vectors. If the particle lags behind the fluid with a relative velocity, fluid would push the particle toward the center, and vice versa. This force can be calculated by Eq. (7).
here, is the magnitude of the velocity gradient. This model was valid only when there was a difference between fluid and particle velocity. Besides, the effect of the wall on the velocity profile around the particle was omitted by the assumptions of this model. As a result, this model was unable to resolve the motion of particles near the wall.
Lateral force on a deformable sphere in an inertia-less steady shear flow was observed by Tam and Hyman [20]. A rigid particle does not experience any lift force in a creeping flow; however, a deformable particle is affected even in a Stokes flow regime. To characterize the deformation of a particle, it was assumed that the elastic displacement was sufficiently small; hence, the linear elastic theory was used, and the total transverse force was obtained as Eq. (8).
and are constant parameters that are available in [20], and and are Lamé constants. Through particle tracking analysis, it was shown that particles reached an equilibrium position for both simple shear and Poiseuille flow, independent of the initial particle position, where particles reside on the centerline for simple shear flow and 0.3 times half the channel height for Poiseuille flow for . Moreover, a general lateral force (Eq. (9)) was proposed, in ≪ 1 which is the ratio of the particle diameter to the channel height, is for 4 (1 -2 ) Poiseuille flow and for simple shear flow, is the shear rate, and G 1 as well as G 2 are numerically predefined functions in terms of s, which is illustrated in Fig. 2 [21]. Fig. 2 The schematic view of a cylindrical particle inside the confine channel flow Also, it has been illustrated that rigid particles ( in a second-order fluid (a well-known ≪ ) viscoelastic fluid model) migrated toward the centerline due to the induced force from normal stress differences in Poiseuille flow, and toward the outer cylinder in the case of circular Couette flow [22]. When the ratio of Weissenberg number (Wi) to Re is much more than unity , the effect of inertial terms in the Navier-Stokes equations are negligible ( ≫ 1) where is the blockage ratio, and are the first and second normal stress differences, 1 2 respectively, and s is normalized vertical location of the particle inside the domain (Fig. 2).
In 1989, Schonberg et al. presented that in a 2D Poiseuille flow, by increasing Re, particles tend to migrate toward the wall [23]. However, the particle was assumed to be small enough so that the velocity profile inside the channel remains undisturbed, indicating that the particle . The singular perturbation technique was employed to solve coupled flow and particle ≪ 1 equations in this article. It is demonstrated that the particle equilibrium position becomes closer to the wall as Re increases, which was validated with experimental data of Segre and Silberberg [13].
The majority of studies so far have dealt with particle migration in the small and intermediate Re flow regime (Re<3000). However, in 1999, Asmolov investigated the lateral migration of buoyant and non-buoyant particles in 2D confined flow over a wide range of Re and presented some results for unbounded and semi wall-bounded cases [24]. It was shown that two equilibrium positions exist for neutrally buoyant particles. It was mentioned that net inertial lift force exerted on the neutrally buoyant particle in 2D Poiseuille flow could be determined using Eq. (11), where is the lift coefficient that is a function of particle lateral position and Re, is the fluid density, is the local shear rate, and is the particle diameter.
While many studies have addressed the migration of particles in planar fluid flow, it is crucial to investigate the effect of a realistic 3D fluid flow on particle migration in confined flow. To address this issue, Matas and his colleagues were one of the pioneers who analyzed 3D particle inertial migration over a wide range of Re  in a pipe flow [25]. Accordingly, the matched asymptotic expansion was used for cylindrical pipe flow, which was previously used for the planar case [24]. Based on their results, lift force in a pipe flow was qualitatively similar but quantitatively different from the one in a planar channel flow. Furthermore, the equilibrium position of particles in the planar flow was closer to the wall, and the magnitude of the lift force in planar flow was considerably larger than that of pipe flow. For Re>700 and in the case of pipe flow, one additional equilibrium position exists [26].
Previous matched asymptotic solutions and point particle methods assumed that the particle is small enough to not disturb the fluid flow. However, particle size has a vital role in particle migration. Therefore, these methods cannot fully capture the particle trajectory. To meet this demand, Hood et al. developed a new theory for the exerted force on finite-size particles within a 3D microchannel flow with a square-shaped cross-section [27]. They employed a numerical method to analyze the dominant balancing forces within the Navier-Stokes equations and then derived an asymptotic model to calculate the lift force on spherical particles as a function of their particle size. Their results were valid up to Re of 80, with maximum particle size being limited only by the proximity of the particles to the walls. Furthermore, they illustrated that the lift force tended to zero at eight symmetrically placed points around the microchannel section, while four points were stable, but the others were unstable. This observation was in good agreement with the experimental data [27]. Earlier in 2018, Asmolov et al. developed a new theory for the lift force on finite-size particles in a confined microchannel flow for Re < 20 [28]. They showed that this theory might be applicable for the use beyond Re of 20 and other confined flows such as pipe flow. The predicted results obtained from their new model were verified by LBM. The proposed lift force is defined by Eq. (12). = 4 2 ( 0 ( ) + 1 ( ) + 2 ( ) 2 ) (12) where is the particle radius, is the maximum shear rate at the wall, is slip velocity, and values are listed as follows.
This study proposed one of the first models for inertial migration of finite-sized particles in confined microchannel fluid flow.
In particular, this model shows a remarkable increase in the lift force near the walls.

Navier-Stokes-Based Solutions
Semi-analytical and asymptotic solutions based on perturbation theories can explain the physics of particle motion since they provide an explicit formula for the forces acting on the particle [29][30][31]. However, these methods have some restrictions that limit their application to a variety of scenarios. These restrictions include particle size, Re, and the distance between the particle and the wall [32][33][34][35]. There are many cases of particle motion that do not comply with these restrictions [36][37][38]. For these cases, numerical methods are powerful alternatives that can determine the motion of particles for a wide range of sizes and Re [39]. Most existing commercial CFD packages employ the point particle assumption when simulating particle motion. Since there is not an exact formula for calculating inertial lift forces acting on a particle, these packages cannot predict the particle trajectory in inertial flows precisely. DNS can be considered as an alternative and robust method, in which forces acting on a particle are calculated from the interaction between the fluid and the particle. The most frequently used methods, which adopt DNS to simulate particulate motion, are briefly explained in the following.

Flow at Specific Particle Position (FSPP)
In some studies, the steady-state flow around the particle and inertial forces acting on the particle at a specific point have been studied, obviating the need for calculating the whole particle trajectory. Accordingly, in 2009, Di Carlo [40] proposed a DNS method for calculating the steady-state flow fields around a single particle at specific positions of the channel crosssection to obtain the net inertial force acting on a particle. This method is highly efficient as it avoids problems encountered in the calculation of transient particle motion such as remeshing at each time step. To calculate the steady-state flow field around a single particle in the position (y,z) of the channel cross-section (Fig. 3A), a part of the channel with sufficient length should be considered (usually 20a is adequate for neglecting the effect of the particle-particle interaction, where a is the diameter of the particle). Then, the particle is located at the center of this part at position (y,z). The boundary conditions shown in Table 1 are implemented in the channel. In the beginning, the particle is considered to be at rest. Using continuity and Navier-Stokes equations (Eqs. (16) and (17)), pressure and velocity fields around the particle are calculated. Next, the values of the forces and torques exerted on the particle are obtained. Linear and angular velocities of the particle are calculated via Newton's second law (Eqs. (18) and (19)).
In these equations, and are the particle surface and outward unit normal vector, respectively. ∂ At the next time step, values of the boundary conditions are updated by velocities obtained from the previous time step, and the flow field is solved again. This iterative algorithm is terminated when the momenta in y and directions and the force per surface area of the particle become smaller than a specified value. Finally, components of the inertial lift force acting on the particle in the y-z plane, i.e., and at point (y 1 , z 1 ) are calculated. By repeating this process for different positions of the channel cross-section, the distribution of inertial forces in the cross-section of the channel can be obtained.

Arbitrary Lagrangian-Eulerian (ALE) method
In the ALE method proposed by Hughes et al. [41], the system of suspended rigid particles in an incompressible fluid flow is considered. The surface velocity of each particle is = + Ω . Here, , and are the particle's center velocity, angular velocity, and center × ( -) , Ω position, respectively. Fluid flow is solved using incompressible continuity and Navier-Stokes equations. In most studies, a finite element scheme is used for solving the flow. At the next time step, particle positions, orientations, and their velocities are updated by integrating the total stress on the surface of each particle and Newton's second law (Eq. (18)and (19)). A new mesh should be generated at this time step to recalculate the flow field. Using this iterative solution, one can track the particle trajectories at each time step. The difficulties of this method are updating the mesh at each time step, which is a time-consuming process. In this method, if the density of a particle is less than that of the fluid, the above explicit time-stepping method leads to an unstable solution. This instability can be overcome by separating part of the drag force, which contributes to the acceleration of fluid around the particles from the total drag force [42]. For more information about the ALE method and methods for stabilizing its solution, readers can refer to [43].

Fictitious Domain Methods
As stated in the ALE method, the computational domain is re-meshed at each time step when particles move along the domain. In fictitious domain methods, a fixed mesh is used for the simulation of particulate flows. Unlike the ALE method, this mesh includes the domain inside the particles, and the flow field is solved for the entire computational domain using the Navier-Stokes equations. Moreover, to ensure that the fluid inside each particle obeys the rigid body motion, a force is applied to its domain. Fictitious methods are classified according to this applied force. If this force is applied to the surface of the particle, it is called the immersed boundary method (IBM), while if it is applied to the body of the particle, it is called the immersed body method. Here, we only describe the Distributed Lagrange Multiplier (DLM) method as a subset of the immersed body method that is frequently used for simulation of particle motion in inertial microfluidics.

Distributed Lagrange Multiplier (DLM)
The DLM method is an immersed body type of the fictitious domain method introduced by Glowinski et al. [44] in the framework of the finite element method. In summary, to solve the flow, one mesh for the pressure field and one finer mesh for the fluid velocity field are required. Particle domains are discretized with additional meshes, which move with their corresponding particles. The size of the fluid mesh should be smaller than that of the particle mesh to avoid the overestimation of the system. The flow field is solved over the whole computational domain, including the particles' interior. A distribution of Lagrange multipliers as body force densities is used to force the flow inside the particles to show a rigid body motion. The value of the Lagrange multipliers is calculated using the constraint at each point of the mesh inside the particle. For more = + Ω × ( -) details about this method, please refer to [43].

Immersed Boundary Method (IBM)
In IBM, there are two computational meshes [45][46][47][48][49]. The first one, which is for the fluid, is a fix staggered Cartesian grid which is referred to as the Eulerian grid. The second mesh is represented by a set of markers, , which distribute evenly over the surface of the ∀ 1 ≤ ≤ particle. This set is referred to as the Lagrangian mesh (Fig. 3B). A Lagrangian mesh has an element size of and the makers are located at the center of each element. Since in the IBM, ∆ the no-slip boundary condition is not imposed at the surface of the particle explicitly, a force distribution should be defined on the particle's surface and added to the Navier-Stokes equation to guarantee the constrained rigid motion of the particle. This constraint should be equal to the local particle velocity in Lagrangian coordinates (Eq. (20)). Moreover, because there are two different meshes, the value of force and velocity between Eulerian and Lagrangian meshes should be transferred using interpolation based on the regularized Dirac delta function δℎ ( ) [45]. To calculate the force distribution, first, the velocity predicted on the Eulerian mesh is interpolated at Lagrangian markers (Eq. (20)). Second, the forces on the Lagrangian mesh are computed based on the difference between the interpolated velocity and the particle velocity (Eq. (21)). Finally, the force is added to the Navier-Stokes equation to calculate the flow field and the motion of particles using Newton's second law [45,46].

Fig. 3
A) Schematic illustration of particle motion in a part of the channel. The surface of particle rotates with an angular velocity of Ω, and walls have the backward velocity of . Reprinted from [50] with the permission of AIP Publishing. B) Representation of the Eulerian and distribution of Lagrangian grids. Eulerian grid is shown with fix staggered Cartesian grid with gray color. Lagrangian gird is represented by a set of green markers, which are distributed evenly over the surface of the particle with an element size of ∆ .

Straight channels Particle motion in a tube and Couette flow
The first attempts to numerically investigate particle motion in a shear flow was made by Feng et al. [36]. They considered neutrally and non-neutrally buoyant particle motion in 2D horizontal Couette and Poiseuille flows, which was followed by similar studies for 3D cases [51]. These studies used the ALE method for the simulation of particle motion and illustrate that the density difference between particle and fluid plays a crucial role in the lateral equilibrium position of particles.
Similarities in the particle motion between confined and unconfined flows motivate several groups to propose lift formulas analogous to the lift force formula acting on a particle in classical aerodynamics using ALE [39,52] and DLM methods [39]. In these studies, since the motion of a single particle was investigated, ALE method was a better choice to obtain more precise results than DLM.
Shao et al. [37] investigated particle motion in a tube at different Re by introducing some modifications to the DLM method. They found that by an increase in Re, the equilibrium position shifts from the Segre and Silberberg equilibrium position to the inner radius equilibrium position (Fig. 4A). The critical Re in which the inner equilibrium position becomes stable is a function of the particle size and the distance between each particle in the flow direction.

Particle motion in channels with a rectangular cross-section
The physics of inertial microfluidics in channels with rectangular cross-sections is more complicated compared to tubular flows. Liu et al. [38] and Mashhadian and Shamloo [53] investigated the migration of particles in a straight channel with a rectangular cross-section using the FSPP method. Particle migration in this channel can be divided into two steps. First, particles enter the channel at random positions (Fig. 4BI) and focus in a line around the walls Experimental results of different studies illustrate various focusing patterns at the outlet of the channel [54,55]. To address these diversities, several researchers numerically investigated the migration of particles in straight rectangular channels [38,40,53,56,57]. Hence, the value of inertial lift force ( ) at different positions within the cross-section should be investigated.
Here, acting on the particle in a rectangular channel is a function of particle size ( ), the aspect ratio of the cross-section ( ), Re, and the position of the particle in the channel crosssection. Di Carlo and his colleagues [40] were pioneers in investigating the distribution of the inertial lift force in the cross-section of square channels using the FSPP method (Fig. 5AI). The which is calculated using the FSPP method for three different sizes of the particle at a constant Re. By increasing the particle size, the equilibrium position approaches the channel center [57], while by increasing Re, this equilibrium position moves toward the wall [40]. Also, in rectangular channels, depending on Re, there are several equilibrium positions near the channel walls [38,53]. At low Re, the equilibrium positions at the center of long walls are stable (Fig.   5DI). By increasing Re, the equilibrium positions near the center of short walls also become stable (Fig. 5DII). Moreover, if Re increases to high values, new equilibrium positions emerge near the long walls ( Fig. 5DIII and Fig. 5BII). The focusing pattern of particles in a channel with a square cross-section at low (I) and high (II) Re [56]. C) By increasing the particle size (blockage ratio), the equilibrium position of particle moves toward the center [57]. D) Distributions of lift force acting on the particle in a rectangular channel with an aspect ratio of two under various Re and particle sizes Reproduced from Ref. [38] with permission from The Royal Society of Chemistry. I) Only the centers of the long walls are stable equilibrium positions (the blockage ratio and Re=100). ( / ) = 0.3 II) Stable equilibrium particle position at Re=200 and blockage ratio of 0.3 III) By increasing Re to high values (Re=200) for the blockage ratio , new equilibrium positions emerged near the ( / ) = 0.1 centers of long walls.

Effect of the particle on fluid flow
In the majority of studies on particulate flows, the effect of fluid on the particle has been investigated to trace the particle migration. However, there are limited works investigating the effect of particles on the carrier fluid, which is of utmost importance to fully understand the exact mechanism of particle focusing within inertial microfluidic systems [58][59][60].  [58]. On the other hand, outside the closed streamline area, the open streamline area can be observed where streamlines after passing the obstacle return to their previous directions. Moving to an inertial regime causes the closed streamlines around the particle to collapse. This creates two reversed streamline areas (Fig. 6B). Both the spiraling streamline and reversed streamline areas meet each other at a stagnation point [61]. Fig. 6A and B are obtained via numerical simulations using the front-tracking finite difference method [58].
Unlike open boundary shear flows, reversed streamlines are observed in confined particulate flows in the Stokes regime (Fig. 6C) [59]. The distribution of streamlines around a particle in A secondary flow in the channel cross-section is generated once there is a particle in the channel, even the particle has no rotation and is located at the channel center [60]. Fig. 6E, obtained through the FSPP method, represents the value and direction of secondary flows created due to the presence of a particle located at different distances from the channel crosssection center. When the particle is closest to the channel center, the weakest secondary flow is observed.
The pressure distribution across the surface of a particle plays a significant role in the lift force acting on the particle. Numerical studies show that there are four distinct regions on the surface; two minimum and two maximum pressure areas (Fig. 6F) [38,50]. Since the magnitude of the pressure distribution across the areas near the channel wall is larger, the particle is repelled from the wall toward the center of the channel. Closed streamline areas near the particle in Stokes flow and B) reversed streamline areas around the particle in the inertial flow [58]. Distribution of streamlines around a rotating particle in confined flow. Reversed streamline observed in both C) inertia-less and D) inertial flow regimes [59]. E) The lateral position of the particle in the cross-section of the channel has a significant impact on the strength and direction of induced secondary flows [60]. F) I. Isometric view of the velocity distribution. II. Top view of pressure distributions on the surface of the particle near the wall. Reprinted from [50] with the permission of AIP Publishing.

Deformable and non-spherical particles
Rigid, spherical particles are commonly used in numerical and analytical solutions of particulate flows. However, the non-spherical shape and elasticity of a particle play a vital role in the inertial particle motion. Using FSPP, Masaeli and her co-workers [62] studied the motion of rigid rod-like particles with different aspect ratios in rectangular channels. They showed that the rotation of particles is a combination of in-plane and out-of-plane rotations (Fig. 7A). By increasing Re, particles mainly exhibit in-plane rotation. Also, by increasing the aspect ratio of a particle, its equilibrium position moves toward the center of cross-section (Fig. 7B). This phenomenon has been used to separate particles with different rotational diameters [63] (Fig.   7C). Also, Lashgari et al. [64] investigated the migration of a single oblate particle in a duct using the IBM method from moderate to high Re (Fig. 7D). To reduce the numerical calculation time, particles are released from an initial position far from the center as the particle migration velocity at the center vicinity is relatively low.
Deformability is another property of real particles (e.g., cells and microgels). Deformable particle models can be divided into three categories: 1) droplet [65,66], 2) deformable capsule [67,68], and 3) elastic particles [69]. The deformability leads to approaching the particle equilibrium position to the channel center (Fig. 7E). When deformable particles are released in the channel, after a while, they reach a steady-state shape and migrate toward their equilibrium position. This feature of the deformable particle is due to the fact that they exhibit tank-treat motion (for periodical deformation of different points of a deformable particle see Fig. 7H) while the particle is rotating rather than tumbling (Fig. 7I), which can be observed in rigid non- Hadikhani et al. [66] using the ALE method observed that by increasing the rectangular channel aspect ratio, bubbles focus at the center of the short walls, which is not similar to the behavior of rigid particles with the same size (Fig. 7G).
There are several constitutive laws for the simulation of the deformation of capsule membranes [70,71]. However, in most cases, the neo-Hookean law is used because of its compatibility with experimental results [67]. The front-tracking method [72], as an accurate method to track the interface of a deformable object, has a formulation similar to that of IBM. However, it requires a dynamic remeshing in each iteration and thus is more time-consuming than IBM.
Using a front-tracking method, Doddi and Bagchi [67] investigated the migration of deformable capsules in Poiseuille flow (Fig. 7H). Near the wall, the deformation rate is high, and after a quick deformation, the particle shape becomes approximately steady. In another study, Raffiee et al. [73] investigated the equilibrium position of cells with a blockage ratio of 0.2 in a square channel using the front-tracking method. According to their results, unlike rigid particles that reach equilibrium position located on the center of the channel walls, the focusing position of the cell is on the channel diagonal at the same range of Re. Villone et al. [69] investigated the motion of elastic spheroid particles under unbound and confined shear flows using the ALE simulation method. The results show that the wall-induced force causes more particle deformation compared to that in an unbound shear flow. toward the center of the channel. Reprinted from [67], with permission from Elsevier. I) Schematic of definitions of tumbling and tank-treat motions.

Inertial microfluidics in non-straight and non-planar microchannels
In the previous section, all reported studies are conducted using straight channels to investigate the physics underpinning inertial microfluidics. However, most channels used for particle manipulations are non-straight channels. In the existing numerical studies on particles' motion in these channels, the distribution of the inertial lift force in the cross-section of the channel has been calculated through the FSPP method. Then, the particle trajectories in the channel can be obtained through the point particle model by combining this force with other existing forces that have explicit formula such as drag force. These trajectories have an explicit formula [38] shown in Eq. (22).  Fig. 8A presents the focusing mechanism in serpentine channels [50]. In this channel, drag and inertial lift forces play a vital role in particle focusing. If the velocity of a particle exceeds a threshold, it will be focused at the center of the long walls of the serpentine channel [74].
Similar to serpentine channels, spiral microchannels have been frequently used for the separation of bio-particles and cells [75][76][77][78][79][80]. If the size of particles in spiral channels is less than a particular threshold ( ), the effect of the drag forces overcomes the inertial < 0.07 forces and traps particles in their vortex streamlines (Fig. 8B).
By changing the shape of the channel cross-section, we can control the focusing positions of the particles in a straight channel (Fig. 8C) [53,81,82]. The equilibrium positions of particles can be obtained through the vector plot of the total force in any arbitrary channel cross-section by combining the calculated inertial forces from the FSPP method with the drag forces ( Fig.   8D) [83]. Therefore, it can be stated that particle focusing is highly related to the channel shape and cross-section, and a deep understanding of the underlying physics helps the investigator for better manipulation of the particles. IV. Particle migration along the length of the channel [77]. C) Particle focusing in a hybrid straight channel with rectangular, triangular, and semicircular cross sections. Reprinted from [53], with permission from Elsevier. D) I. Force-maps for triangular microchannel over a wide range of Re. The magnitude of lift forces for II. Re=20 and III. Re=150 [83].

Inertial microfluidics in non-Newtonian fluids
Particles in a non-Newtonian fluid experience several lateral forces in addition to the inertial forces. The magnitude and direction of these lateral forces depend on the non-Newtonian fluid's rheological properties. Therefore, the analysis of particle motion in a non-Newtonian fluid becomes more complicated compared to the Newtonian fluid. Seeing as this field is a subject of ongoing investigations, only recent major studies on particle motion in viscoelastic fluids have been reviewed.
Similar to Newtonian fluids, to obtain the flow field in a viscoelastic fluid, continuity and momentum equations need to be solved (Eqs. (23) and (24)). However, an extra stress tensor ( ) must be added to the total stress tensor ( ) in order to consider the effect of viscoelasticity τ σ (Eq. (25)) [84].
σ = -pI + 2 ( ) + τ (25) where and are the Newtonian solvent viscosity and strain rate tensor, respectively. To ( ) solve this set of equations, the relation between and other parameters, such as fluid velocity, τ should be obtained by a constitutive equation [85,86]. Two constitutive laws that are frequently used are Oldroyd-B and Giesekus equations [87]. Equation (26) shows the Giesekus constitutive law [84]. In this equation, , , and are the mobility factor, polymeric viscosity, and relaxation time, respectively. When the mobility factor ( ) is equal to zero, the Giesekus equations is reduced to the Oldroyd-B equation.
The migration and focusing pattern of the particles in the channel cross-section mainly depend on the competition between elastic and inertial forces acting on particles (Fig. 9A) [88][89][90][91]. Li and his colleagues [92] investigated the effects of fluid elasticity, fluid inertia, and shearthinning viscosity using Oldroyd-B (Fig. 9BI) and Giesekus equations (Fig. 9BII). They used the DLM method in the framework of the finite volume method on a staggered grid. Particles in Oldroyd-B fluid focus at the center of the channel, while in the Giesekus fluid particles reach their equilibrium at a position away from the center (Fig. 9B).
More recently, Raffiee and colleagues using the DLM method [89] studied the focusing pattern of particles in the cross-section of a square channel with Giesekus fluids for different values of Wi and Re (Fig. 9C). In another study, with a modified version of the DLM method, Yu et al. [90] investigated focusing patterns of particles in Oldroyd-B fluids. Fig. 9D shows the result of their study for a particle with a blockage ratio of 0.15 in a square channel under various conditions. Using FSPP, Raoufi et al. [88] investigated the effect of channel cross-section on focusing efficiency of elasto-inertial flows in straight channels. They have shown that by increasing the channel corner angle, the elastic force pushes the particles toward the center of the cross-section more efficiently.
To simulate the motion of cells in a square microchannel in Newtonian and viscoelastic fluid ( Fig. 9E), Raffiee et al. [93] used the front-tracking method. The results revealed that an increase in the solution concentration led to significant enhancement in the volumetric flow rate, which can help to boost the total throughput of a microfluidic device.  [88], with the permission of AIP Publishing. B) Investigation of fluid behavior and particle equilibrium position in I. Oldroyd-B II. and Giesekus fluids [92]. C) the distance of the offcenter lateral position of particles from the center of the channel along the I. main axis and II. diagonal of the channel for all ranges of Wi and Re. Stable equilibrium positions are identified by filled symbols, whereas unstable ones are represented by unfilled symbols [89]. D) focusing pattern of a particle with a blockage ratio of 0.15 under different conditions of EI and Re [90]. E) Increasing in Wi leads to relocation of the cells toward the centerline [93].

Lattice Boltzmann Method for Inertial Particle Migration
The LBM, as an explicit alternative method for fluid dynamics problems, was introduced in 1988 by McNamara and Zanetti [94]. In this method, the fluid is approximated by discrete particle distributions moving and colliding on a regular lattice. LBM can be easily adapted for parallel processing. The numerical investigation of particle dynamics in fluids requires tracking of a solid-fluid interface where LBM is particularly powerful [95][96][97][98]. As providing the details of LBM for inertial particle microfluidics is beyond the scope of this review, please refer to textbooks, such as [99,100].

Inertial particle motion in 2D/3D straight channels
Ladd and his team were the first to investigate inertial phenomena in microfluidic devices [101][102][103], where they investigated the equilibrium positions of particles at Re of 100 to 1000 in a square duct [104]. At low Re, the numerical results were consistent with experimental ones, while at high Re, some discrepancies are raised (Fig. 10A-D). Afterward, Nakagawa and his colleagues [105] addressed this problem in their numerical study. They illustrated that for Re less than 260, the equilibrium positions were at channel faces, while for higher Re, channel corners were added to the equilibrium positions. Also, they showed an outward shift of the equilibrium positions at channel faces when Re was increased. For larger Re, the positions moved back inward (Fig. 10E-H). By combining IBM and LBM, the physical mechanism for this behavior was identified. The two vortices generated next to the particles grew upon an increase of Re, which pushes away particles from the wall [106,107].  [104], at Re = 100, particles migrated to the eight equilibrium positions. For Re = 500, particles migrated to one of four stable locations near the corner of square duct, while for Re =1000, the particle equilibrium configuration changed. Reprinted from [104], with the permission of AIP Publishing. Map of lateral forces for Re of E) 260 and F) 514 calculated by Nakagawa and co-workers [105] who showed that the particle behavior and equilibrium positions change by increasing Re. Equilibrium position in G) channel face and H) channel corner. Circles are the results of their study, while experimental results are indicated by other shapes.
In a simple-shear flow profile, Mao and Alexeev evaluated particle inertia, fluid inertia, and their combination for spheroid particles with different aspect ratios at low and moderate Re [108]. The authors pointed out the possibility of superposition principle for the effects of fluid inertia and particle inertia at sufficiently low Re. They showed that Stokes number, aspect ratio, Re, and initial orientation of a single spheroid affect its trajectory. Stokes number (Eq. (27)) defines the ratio of particle response time to fluid response time This number is essential for the separation of non-biological particles that are denser than water. At low Stokes numbers, particle behavior is similar to a neutrally buoyant particle [109], while an increase in Stokes number leads to the oscillatory behavior of particles [110]. Mao and Alexeev proposed a microchannel decorated with diagonal symmetric aligned ridges, resulting in the generation of secondary flows, for the hydrodynamic sorting of microbeads ( Fig. 11A) [111]. The induced secondary flows and inertial migration led to the separation of microparticles based on their sizes (Fig. 11B). However, the values of Re investigated are relatively small, and it is unclear how relevant inertial forces are. Besides ridges, wall roughness also affects the particle trajectory [112].
With the combination of IBM, LBM, and FEM, the hydrodynamic focusing of rigid particles ( in a straight channel was investigated [113]. The particle trajectory (Fig.  = 3, 6, and 12 µm) 11C) for two particles with a radius of 6 and 12 µm at Re = 83 reveals that while these two different particles show small oscillations in the interacting stage, both followed two different paths at the separation stage. Besides, particle equilibrium position ( demonstrated = 6 µm) that for high enough Re, inertial effects became dominant and pushed particles closer to the channel wall. However, an even larger Re might result in unfavorable focusing and unstable fluid dynamics (Fig. 11D). In another study, Krüger and his co-worker investigated the interplay of inertial and deformability over a wide range of Re (3-417) and Ca (0.003-0.3) for the volume fraction of [114]. They showed that for a fixed Re, softer particles focused = 0.1 near the center-line of the microchannel (Fig. 11EI). According to their results, the higher the rigidity of the particles, the thinner the depletion layer (d), which is defined as the minimum distance of the particle surface to the wall (Fig. 11EII). For Re > 45, as identified with the dotted line in Fig. 11EII, the growth of the depletion layer is significantly increased.  Fig. 11 A) Schematic illustration of a microfluidic device decorated with diagonal ridge developing secondary flows. B) The particle trajectory of various particles from the starting point of y/w=0.5 was evaluated. As can be seen, the smaller particle migrated toward the positive y while larger ones migrated in the opposite direction. Different lines belong to different particles [111]. C) particle trajectory for two particles with a radius of 6 and 12 µm. The whole process can be divided into three separate sections of initial, interacting, and separation stages. In the separation zone, particles with different size migrated By increasing Re, equilibrium positions relocate from the center of the longer wall to the center of the smaller wall, and it pushes particles more toward the walls. Reprinted from [113], with permission from Elsevier. E) I. for a fixed value of Re, softer particles are closer to the centerline II. an increase in stiffness leads to an increase in the depletion layer. Re > 45 is identified by a dotted line [114].
The entropy of particle focusing in microfluidic devices (Eq. (28)) demonstrates the ( ) connection between each particle to the cumulative performance of several particles, where probability function was identified as and the total number of states in the system (either ( ) fluid or solid state) was recognized as [115]. A higher ordering degree means lower focusing entropy, indicating a better focusing behavior.
Focusing entropy of rigid particles was lower than soft ones (Fig. 12A), signifying better ordering behavior of rigid particles, and the rectangular cross-section showed more viability regarding hydrodynamic focusing than the square and circular ones. With the coupling of LBM and LSM, Kilimnik et al. investigated the migration of a deformable capsule [116]. The larger and the softer the particles are, the closer to the center they become, indicating that deformability leads to a center-facing lift force, which is a well-known effect. In addition, increasing the viscosity of the encapsulated fluid resulted in equilibrium positions closer to the wall. It has also been illustrated that in a Poiseuille flow, soft particles moved away from walls, whereas the position of hard ones is independent of the Re. Besides, the migration velocity in Poiseuille flow was 3 to 4 times higher compared to the simple shear flow due to the higher variation of velocity across the particle [117]. Prohm and Stark showed that the lateral position of particles could be manipulated by axial control forces ( Fig. 12BI and BII) [118]. In addition, they used the feedback control approach to apply a non-constant axial force in order to efficiently control the position of particles and increase the particle throughput (Fig. 12BIII). Based on their method, this group [119] extended their theory by investigating of deformable capsules in a microchannel by evaluating on Laplace number (Tutorial Box 1). While deformability entirely is depended on Re, the equilibrium position approximately is independent of Re and falls into a master curve line (for each particle diameter) when plotted versus Laplace number. Moreover, they used external forces to control the equilibrium position of particles. Although rigid particles move away from the centerline in high Re, very soft capsules behave oppositely. Moreover, it was shown that by choosing the proper amount of flow rate, separation of a deformable capsule with different aspect ratio and membrane shear elasticity through bifurcation is achievable [120].
The interparticle spacing in inertial microfluidics is of significant importance for such applications as imaging or flow cytometry [121]. The favored spacing in low particle Re was measured to be 5D while at high particle Re, it went to 2.5D, where D equals to particle diameter. However, an increase in concentration led to a decrease in particle spacing such that a single train of the particle could not be identified [122]. To shed more light on this matter, Liu and Wu [123], in a 2D lattice domain, proposed a dimensionless focusing number and represented that resulted in complete inertial ( = , = 0.36, = 2.33) > + migration while shows unfocused particles and led to partial particle < --< < + focusing ( and are upper and lower limits of particle focusing, respectively). Besides + the train of particles, Haddadi and Morris developed a study regarding dynamics and trajectory of isolated and suspended (solid volume fraction was less than 0.3) pair of particles [124]. They showed that pair trajectories and streamline around an isolated particle have similarities, including reversing, in-lane spiraling, off-plane spiraling (Fig. 12C), and open but fore-aft asymmetric streamlines. More recently, with more emphasis on inertial microfluidics, Schaaf and his colleagues investigated on flowing pair of particles [125]. They believed that the inertial lift profile is strongly dependent on the position of the particle in the fluid, whether it is leading or lagging (Fig. 12D). Nonetheless, by increasing the axial distance or Re, the profiles become similar to each other while having a constant shift to each other. ) rectangular microchannel. The average values of entropy after initial falling were highlighted by dash lines. Reprinted from [115], with permission from Elsevier. B) Schematic illustration of the feedback control. I. when the particle left the desired area (i.e., ), the axial force turned on, and it was [ -, ] turned off when the particle went to the centerline (i.e., ). II. the velocity of a particle with and = 0 without axial force control in the desired portion (i.e., III. An example of particle trajectory in [ -, ] the desired interval. Reproduced from Ref. [118] with permission from The Royal Society of Chemistry. C) off-plane spiraling streamlines for an isolated pair of particles for I. Re = 0.05 and II. Re = 0.6. An increase in Re leads to decrease in the off-plane spiraling zone [124]. D) representation of color-coded lift profile for a pair of leading and lagging particles [125]. Published by The Royal Society of Chemistry.

Inertial particle motion in non-straight microchannels
Generally, straight microchannels enjoy the features of a simple fabrication process and easy operation, and the inertial focusing mechanism is almost clear for certain specific crosssections. However, these channels suffer from the long footprint, which impedes their further applications and makes limitations on their commercialization aspect. As an alternative, scientists set to induce secondary flows within the channel by means of obstacles, changing channel geometry, or using curved channels to assist inertial particle motion. Accordingly, the principal mechanism of particle migration becomes more intricate, requiring a robust, solid description. So far, the number of studies using LBM for the investigation of inertial particle migration in a non-straight channel has not been significant. Due to the intricate nature of flow in curvature, it is challenging to extract the governing equations of fluid with the corresponding dominant forces acting on particles. Moreover, most of results focus on experimental studies, resulting in unclear particle-particle and particle-fluids interaction during the focusing process and unspecified particle migration in cross-section of the channel. Hence, it is anticipated that in the near future, more numerical studies focus on inertial particle migration within nonstraight microchannels. In the following, we review studies conducted using LBM in a nonstraight microchannel.

Serpentine microchannels
Serpentine microchannels have a smaller footprint compared to straight channels and bear the feature of massive parallelization. Gaining the efficiency of IB-LBM, Jiang and colleagues investigated numerically and experimentally particle focusing with a diameter of 5 and 10 µm inside a symmetrical serpentine microchannel [126]. Calculating fluid flow by LBM and particle structure by FEM, the authors integrated these two parts using IBM. In low Re, inertial lift force increased, and particles were pushed toward the sidewalls where smaller particles were closer to the walls. In high-enough Re, the dominant force turned to be drag force, Dean flow forced particles to swing out of inertial lift force trap, and particles were focused at the vicinity near the center of Dean flow vortex ( Fig. 13AI and Fig. 13AII). This focusing pattern is expressed by Eq. (29) where is curvature ratio ( , is hydraulic diameter and = ℎ 2 ℎ is the channel radius) and is particle diameter [127,128].
Based on Eq. (29), an increase in Re leads to make the Dean drag force dominant.

Cavities and Contraction-expansion arrays
Cavities, contraction-expansion arrays, and constrictions provide additional secondary flows, assisting in inertial particle migration [129]. Wu et al. proposed a contraction-expansion microchannel to induce additional secondary flows for particle separation (Fig. 13BI and BII) [130]. They observed that particles with the diameter of 9.9 µm moved more toward the centerline than those with a diameter of 5.5 µm. However, since the concentration of the used microparticles for numerical simulation was not high enough, particle-particle, particle-wall, and particle-fluid interaction were not thoroughly investigated. Also, although contractionexpansion channels with higher expansion to contraction ratio required higher flow rates to focus size-based particles, their focusing performance is better [131].
In another study, the vortex entrapment of particles inside a microchannel with flat and curved edges was evaluated [132]. The authors suggested that the combination of repulsive forces and inertial effects could potentially lead to liberating particles trapped in a vortex zone at high Re.
It was shown later that the entrapment inside a microcavity is related to both particle dynamics and flow morphology [133] where the feasibility of these microcavities was showcased by cancer cell separation [134]. Using IB-LBM, it was revealed that the dynamics of a particle within a cavity is affected by two competitive outward centrifugal and inward inertial forces [135]. There are three particle entrapping phases: no trapping (Re < 50, due to the lack of enough inertial forces), stable trapping (50 ≤ Re ≤ 200), and unstable trapping (Re > 200, the existence of strong fluid inertia) (Fig. 13C). Also, four trapping modes of outer to inner (50 ≤ Re < 100), invariable (100 ≤ Re < 150), inner to outer (150 ≤ Re ≤ 200), and inner to escape (Re > 200) within a microcavity occur. Within a cavity, rotating and orbiting velocity of a particle is not constant while both motions are counter-clockwise. Fig. 13 A) The comparison of the lateral position of particles with different sizes in a channel. Data shown in the figure is obtained experimentally and numerically. Reproduced from Ref. [126] with permission from The Royal Society of Chemistry. B) I. Schematic illustration and the design principle of the proposed contraction expansion microchannel by Wu and co-workers. II. randomly dispersed particles were first focused due to the inertial effect. Next, at the presence of contraction-expansion arrays, secondary flows were induced resulted in the migration of small particles to the side of channel walls while maintaining large particles near the centerline of the channel. Afterward, each stream was collected from separate outlets. Reproduced from Ref. [130] with permission from The Royal Society of Chemistry. C) Three particle entrapping phases I. no trapping, II. stable trapping, and III. Unstable trapping.

Concluding remarks and outlook
In this review, we have summarized all computational techniques for inertial microfluidic modeling and categorized them into three subsections of semi-analytical solution, direct numerical simulation, and Lattice Boltzmann method. In the first, all relevant articles utilizing semi-analytical methods to simulate inertial particle focusing were reviewed. In these methods, the main important parameters are Re and confinement ratio, allowing an analytic treatment of inertia-induced migration for the calculation of lift force profiles. As these methods require the particle radius to be much smaller than the channel diameter, it is hardly applicable to microfluidic particle flow, where this assumption is often violated. Also, they cannot be implemented for coupled inertia-viscoelastic problems where complex constitutive equations exist. There are some crucial issues that need to be addressed in order to achieve more accurate and realistic solutions. The first one is overcoming the inherent complexity of solving 3D partial differential equations (PDEs) using proper analytical or semi-analytical approaches.
Analytically solving a set of PDEs for the fluid (three momentum equations) coupled with a set of PDEs for particles (three linear momentum plus three angular momentum equations) and also two coupling equations at the interface of fluid and solid domains using non-homogeneous boundary conditions is nearly an impossible task. Therefore, all covered articles tried to simplify this complex problem. The second issue stems from convective terms in the momentum equations. These terms alter prior first-order PDEs to more complex second-order PDEs, implying more difficulty in solving the equations by analytical methods. Last but not least, the effect of particles on fluid flow through disturbances near the particles is another intricacy of inertial microfluidics which should be considered.
In the second section, all numerical studies on inertial particulate flows based on the Navier-Stokes equations are reviewed. These include methods such as ALE, DLM, IBM, and FSPP which are mainly used for assessing the effects of different parameters, such as particle shape, particle deformability, channel geometry, and the type of fluid on particle migration within a microchannel.
Although there has been significant progress in inertial modeling of particle motion, there are still various cases which have not been thoroughly investigated. As can be seen in Table S1,   the most uninvestigated cases are non-straight channels, non-spherical/deformable  interpolation functions [136]. Although there are some efforts (e.g., picturing the steep boundary condition using a logarithmic representation) to address this problem to some extent, this field is actively being investigated. The need for further investigation is more evident when one considers a combination of the above-mentioned conditions, namely deformable and nonspherical particles in non-straight channels in viscoelastic fluids.
It would be significant to assess the capability of existing commercial CFD packages for simulation of inertial microfluidics and give an insight to those who want to use them for their studies. There exist various commercial CFD software such as COMSOL Multiphysics, ANSYS Fluent, OpenFOAM, or Flow-3D. Using these software packages dramatically reduces the efforts required to write bespoke CFD codes. Nevertheless, using these packages for the simulation of particle motion in inertial and viscoelastic flows brings about several problems.
Most of these packages can deal with particle motion, such as the particle tracing module in COMSOL. As mentioned earlier, these modules treat particles as point particles, indicating that forces acting on particles are calculated through a series of equations (e.g., Stokes drag). The problem with this approach is that the exact formula for the inertial forces acting on the particle in terms of flow parameters and particle properties is unknown. Although COMSOL, as a pioneer in this field, provides an approximate formula for inertial forces, the numerical model in some cases gives incorrect and inaccurate results for 3D flows. Therefore, it is necessary to calculate the interaction between the finite-size particles and the fluid. This requires fluidstructure interaction (FSI) which is in principle available in most commercial software packages. However, since current FSI modules do not have periodic boundary condition for moving particles at the inlet and outlet of channel unit cell, they cannot be implemented to the particle motion simulation. For simulating viscoelastic fluids, it is appropriate to use constitutive equations such as Giesekus or Oldroyd-B that are consistent with experimental results; but these models are not defined by default in current software packages. Altogether, it is not possible or practical to accurately simulate particle motion in inertial and viscoelastic flows using commercial software packages in their default mode, and additional scripts must be written. Using these scripts, it is possible to use periodic boundary conditions in a part of the domain or to add constitutive equations for viscoelastic fluids to commercial packages.
In the third section, all numerical studies based on LBM were reviewed. Initially, numerical results obtained using LBM were not entirely in line with those obtained from experiments.
More recent results, however, illustrate that LBM can now accurately predict the particle behavior in inertial microfluidics. LBM has a strong potential to be applied to inertial microfluidics and can be considered as a robust, efficient, and powerful alternative to conventional Navier-Stokes solvers in many fluid flow problems. LBM is mostly limited to viscous fluids (both Newtonian and non-Newtonian), and more research is needed to extend the method to viscoelastic fluids. Since the LBM in its original form depends on a cubic lattice, it is challenging to apply the method to geometries such as spirals due to the large "dead" volume. This explains why the LBM has mostly been used to study inertial flows in straight channels. Sparse geometry LBM codes can be applied to geometries such as spirals more easily, and more development work in this field is needed.
Due to the large number of different lattice-Boltzmann boundary condition schemes available for moving and stationary objects, it is yet to be decided which methods are most suitable and reliable for inertial particle microfluidic applications. Since the LBM is a weakly compressible scheme, it is not very accurate for the pressure field, which may negatively affect the lift and drag forces acting on moving particles. Although recent LBM studies have shown considerable progress in enhanced computational performance, non-spherical particles (using DSP-LBM method), or non-Newtonian fluids [137][138][139][140], more research is needed to investigate the accuracy and reliability of the LBM for inertial particle microfluidics.
In conclusion, computational inertial microfluidics is a nascent field, requiring more devoted studies to develop and describe the underlying physics. Please refer to Table S1 in electronic supplementary information for the current contribution of computational methods in inertial microfluidics. Most of the studies published so far are not comprehensive since the developed codes and models are limited to specific channel geometries or particle properties, mostly rectangular straight microchannels and rigid spherical particles. New computational packages, commercial and scientific, need to be developed that can predict the migration of realistic particles in complex flow geometries, depending on initial conditions, flow parameters, and particle concentration.

Dimensionless parameters in inertial microfluidics Reynolds number (Re):
Re is one of the most important dimensionless numbers in fluid mechanics, representing fluid behavior in various situations. Re is the ratio of inertial to viscous effects, as shown in Eq. (S1). In this text, we refer to channel Reynolds number as Re. Here, is the fluid density, is the fluid dynamic viscosity, is the mean velocity of the fluid, and is the hydraulic diameter of the channel. where is particle diameter and is channel dimension.

Blockage Ratio ( ):
The blockage ratio is defined as the ratio between particle diameter to the characteristic length of channel (height of the channel cross-section). If the blockage ratio is more than 1, the channel gets blocked.  where is the dynamic viscosity of the continuous phase fluid, is the maximum velocity, and is the channel width. For a deformable capsule or elastic particle, Ca is considered as the ratio between the viscous force and elastic force, while the elastic modulus takes different forms depends on capsules or bulk elastic particles [67,142].

Laplace number (La):
Ca for a deformable capsule depends on the flow characteristics (i.e., flow speed). La is a dimensionless number that characterizes the rigidity of a capsule based on the elastic shear force ( , is shear modulus) and intrinsic inertial force scale ( ), without depending on the explicit flow speed [119]. 2

Using
, La can be written as Eq. (S6).