Numerical simulation of inertial microfluidics: a review

Abstract

Since the proposal of inertial microfluidics in 2007, it has been widely used for particle focusing and separation with superior performances. The particle migration behaviour in microchannel is extremely complicated involving various parameters. Recently, computational approaches have been utilized to obtain better insights into the underlying physics, which facilitates a comprehensive and intuitive understanding of particle migration in inertial flow. The numerical methods are ideal ways for assessing the effects of various parameters on particle focusing. In this review, numerous theories and models proposed to systematically explore the effects on migration of particles are summed concisely, especially numerical methods. Besides, latest research advances in direct numerical simulation of diverse particle motion and migration, whether rigid or deformable, spherical or non-spherical, suspended in simple or complex channel, are introduced comprehensively. Then, numerical simulations for microfluidic devices design are provided and discussed to facilitate the optimization of microchannel and the design of coupling microfluidic device for real- particles separation. Finally, we further discuss the remaining challenges in the modeling of inertial particle microfluidics and provide several suggestions for future investigators.

KEYWORDS:

• Inertial microfluidics
• direct numerical simulation
• cell separation
• CTCs

1. Introduction

In the 1990s, Manz et al. (Manz et al., 1990) proposed the concept of a micro-total analysis system and produced a microfluidic chip, which was gradually developed into a microfluidic technology for controlling fluids in micro-nano-scale structures (Kulrattanarak et al., 2008). Due to the advantages of integration, automation and high throughput, microfluidic technology has shown considerable promise in point-of-care diagnostics and clinical studies. Isolation, enrichment, and purification of cells using microfluidic platforms have been a flourishing area of development in recent years with many successful translations into commercial products. Microfluidic technology enriches circulating tumour cells (CTCs) based on biochemical or physical properties, requiring small sample volumes, controllable flow rates, and the ability to capture live cells. Microfluidic technology based on biochemical properties can effectively sort different types of cells with similar shape or size in a short period of time. The phenotypes of CTCs are diverse, including epithelial, mesenchymal and other cell phenotypes, and the phenotypes are prone to dynamic changes with the microenvironment (Cheng et al., 2018 Miao & Tang, 2019;). At present, in most commercial technology platforms, epithelial cell adhesion molecule (EpCAM) is used as the only surface-specific antigen, ignoring the extremely aggressive mesenchymal phenotype CTC that has undergone epithelial–mesenchymal transition (EMT). The label-free CTC microfluidic separation techniques based on physical properties have developed rapidly. According to the technical mechanism, it can be broadly classified as active and passive separation techniques. Active techniques rely on external force fields (e.g. acoustic (Collins et al., 2016), optical (Hu et al., 2019), and dielectrophoresis (Cheng et al., 2015)), while passive techniques (e.g. microfluidic filters (Liu et al., 2018), deterministic lateral displacement (Loutherback et al., 2012), and inertial microfluidics (Kwak et al., 2018)) only rely on the channel geometry and inherent hydrodynamic forces for separating. Among all existing microfluidic systems, inertial microfluidics has emerged as a promising way for particle and cell separation due to its capabilities in high-volume and high-throughput sample processing (Cruz et al., 2019; Gou et al., 2018; Shen et al., 2017; Stoecklein & Di Carlo, 2018; Zhang et al., 2020). The flow channel structure of the inertial microfluidic system can be divided into straight and curved forms. In addition to the inertial forces, the Dean forces in transverse Dean flow arose from curved channel structure not only provide a compact design, but also allow the exquisite control of particle migration (Warkiani et al., 2016). In 2015, Clearbridge BioMedics developed the ClearCell® FX system based on a microfluidic chip with a single curved channel, which has successfully achieved high-activity separation of heterogeneous CTCs using inertial focusing technology. The system was approved as a Class I medical device by the US Food and Drug Administration in 2018, which provides reliable liquid biopsy solutions for clinical research used in vitro diagnostics.

In the last decade, many research groups focused on the application of inertial microfluidics system with curved channels for the separation of CTCs (Dhar et al., 2015; Liu et al., 2015; Ozkumur et al., 2013; Warkiani et al., 2013) or bacteria (Lee et al., 2019; Mashhadian & Shamloo, 2019). However, inertial microfluidics may have a low recovery rate or purity due to the selectivity only based on size. Moreover, the multi-variable characteristics such as microchannel structure, fluid velocity and particle dynamics size, also pose realistic challenges to the design and fabrication of the inertial focusing system. Numerical modeling is a crucial tool to predict and assess the underlying physics and performance of inertial microfluidic devices, which will eventually facilitate their optimum design for a high throughput, recovery rate and purity. Currently, a few review papers have been reported to introduce inertial microfluidics (Chung, 2019; Tabatabaei et al., 2022; WenlaiTang et al., 2020). However, some of them focus on the computational inertial microfluidics with numerous theories (Bazaz et al., 2020), and there is a lack of a comprehensive review of numerical simulations for the design of inertial microfluidic devices. Accordingly, the primary aim of this review article is to provide researchers with the latest updates regarding numerical simulations of inertial particle motion within microchannel. First, we summarize the models developed for inertial particle motion, covering asymptotic calculations, Navier-Stokes-based approaches, lattice Boltzmann method (LBM). Next, we introduced numerical studies of inertial particle motion within microchannel based on those models, highlighting the effects of different parameters on particle migration and focusing on the significance of numerical studies for optimizing inertial microfluidic devices.

2. Numerical models

In 1961, Segre and Silberberg reported the phenomenon that the randomly dispersed particles migrated towards specific equilibrium positions inside a microchannel, which was defined as inertial microfluidics (Segr & Silberberg, 1961). Subsequently, it was extensively explored and numerous articles were reported to evaluate the underlying physics mechanism. The rapid development of computer technology has facilitated the promising prospects of numerical methods in the field of inertial focusing. To date, numerous researchers have provided a multitude of numerical solutions, which can be divided into asymptotic solution, Navier-Stokes-based approaches and LBM. First, the asymptotic solution predicts the inertial forces of a particle within microchannel by simplifying fluid equations. However, the oversimplification of assumptions and constraints, resulting in the phenomenon described by this solution is far from the reality. The Navier-Stokes-based approaches address all prevalent issues in asymptotic solutions including particle size or particle effects on fluid streamline. Nevertheless, the challenge still exists in the solid-fluid interface tracking or the demanding calculation time for solving equations. Over the last decade, increasing researchers have used the LBM for inertial flow modeling as an alternative of others. Due to the relative simplicity of algorithm parallelization, the addition of new physics and the strength in the intermediate Reynolds number (Re) regime, it is extensively explored and applied. For obtaining particle trajectories in complex microchannel with a low computing load, a relatively simple approach was proposed by incorporating the results of asymptotic solutions into the Lagrangian particle tracking (LPT) scheme. This method has the ability to obtain realistic picture of real applications, which facilitates the design of inertial microfluidic devices effectively.

2.1. Asymptotic solutions

Since the inertial microfluidic was reported, many asymptotic solutions have been established to predict this phenomenon. In 1961, Brenner found that resistance against a particle path increased in a wall-bounded case compared to an unbounded case when a particle moved under otherwise identical conditions. This phenomenon is caused by drag force (F_d) on the sphere induced by a solid surface or a free plane surface. In addition to F_d, suspended particles within a microchannel also experience appreciable inertial lift forces (F_L), including shear-gradient lift force (F_LS) and wall-induced lift force (F_LW). The former is formed due to the curvature of the velocity profile and tends to direct the particle away from the channel centre, whereas the latter repels the particle away from the wall. The lateral force field formed by the combination of these two forces ultimately determines the focusing and equilibration of particle within microchannel. In the recent half century, increasing asymptotic solutions have been proposed to solve the force acting on various particles within diverse flow, which are summarized in Table 1. Initially, the studies of asymptotic solutions focused on addressing the migration of spherical particles in planar fluid flow (Ho & Leal, 1974; Ho & Leal, 1976; Rubinow & Keller, 1961; Saffman, 1965; Tam & Hyman, 1973). With the increasing diversity and complexity of microfluidic applications in particles sorting and focusing, it is crucial to investigate the real particle migration in confined microchannel flow. Consequently, a novel asymptotic equation was proposed for investigating the lateral migration of buoyant particles in 2D confined flow over a wide range of Re (Asmolov, 1999). To predict how particle size affects the focusing position, an asymptotic model valid for a wide array of particle sizes and Re was derived. The asymptotic approach allows to compute the focusing forces for particles placed at arbitrary positions in the three-dimensional (3D) channel by developing a perturbation series expansion for the lift force (Hood et al., 2015). Subsequently, the asymptotic theory was evidenced by the direct measurement of particle positions and inertial migration velocities via hybridizing particle image velocimetry (PIV) and particle tracking (Hood et al., 2016). Moreover, the previous asymptotic expansion used for the planar case has been matched accordingly to suit cylindrical pipe flow (Asmolov et al., 2018).

Table 1. Asymptotic solutions.

Years	The equations of asymptotic solutions	Applications	References
1961	$F_d=6\pi bU\left\{{\displaystyle \frac{4}{3}}\sinh \left(\alpha \right)\sum _{n=1}^{\mathrm{\infty }}{\displaystyle \frac{n\left(n+1\right)}{\left(2n-1\right)\left(2n+3\right)}}\left[{\displaystyle \frac{2\text{sinh}\left(2n+1\right)\alpha +\left(2n+1\right)\sinh 2\alpha }{4{\sinh }^2\left(n+\frac{1}{2}\right)\alpha -{\left(2n+1\right)}^2{\sinh }^2\alpha }}-1\right]\right\}$	Fd on the sphere resulted by a solid surface.	Rubinow & Keller, 1961
	$F_d=6\pi bU\left\{{\displaystyle \frac{4}{3}}\sinh \left(\alpha \right)\sum _{n=1}^{\mathrm{\infty }}{\displaystyle \frac{n\left(n+1\right)}{\left(2n-1\right)\left(2n+3\right)}}\left[{\displaystyle \frac{4{\cosh }^2\left(n+{\displaystyle \frac{1}{2}}\right)\alpha +{\left(2n+1\right)}^2{\sinh }^2\alpha }{2\sinh \left(2n+1\right)\alpha -\left(2n+1\right)\sinh \left(2\alpha \right)}}-1\right]\right\}$ $\alpha ={\cosh }^{-1}\left({\displaystyle \frac{h}{b}}\right)=\ln \left\{{\displaystyle \frac{h}{b}}+\sqrt{\left[{\left({\displaystyle \frac{h}{b}}\right)}^2-1\right]}\right\}$	Fd on the sphere resulted by a free plane surface.
1965	$F={\displaystyle \frac{81.2\mu V{a}^2{\zeta }^{0.5}}{v^{0.5}}}$	A spherical particle travels through a highly viscous liquid.	Saffman, 1965
1973	$D_y=4\pi \left[{\displaystyle \frac{9\mu a}{20}}{U}_{\mathrm{\infty }}{C}_{21}^1-5\mu a^2\omega \left({\displaystyle \frac{\frac{27}{20}\mu U_{\mathrm{\infty }}}{aG\left(2+\frac{\lambda }{G}\right)}}+a{B}_{10}^1\right)\right]$	Lateral force on a deformable sphere in an inertia-less steady shear flow.	Tam et al., 1973
1974	$F_L=k^{2}Re\left({\text{β}}^2{G}_1\left(s\right)+\text{β}\gamma G_2\left(s\right)\right)$	Particles travel through simple shear and Poiseuille flow.	Ho et al., 1974
1976	$F_L=-{\displaystyle \frac{80}{3}}\pi k^3{V^2}_{Max}\left(1-2s\right)\left(N_1-2N_2\right)$	Particles travel through Poiseuille flow.	Ho et al., 1976
1999	$F_L=c_l{\rho }_f{G}^2{\alpha }^4$	The neutrally buoyant particle in 2D Poiseuille flow.	Asmolov 1999
2015	${F_L}^{\left(3D\right)}={\displaystyle \frac{\rho U_m^2{c}_4^{\left(3D\right)}{a}^4}{l^2}}+{\displaystyle \frac{\rho U_m^2{c}_5^{\left(3D\right)}{a}^5}{l^2}}+O\left(a^6\right)$	Inertial migration of a rigid sphere in 3D Poiseuille flow.	Hood et al., 2015
2018	$F_L=\rho \alpha G_m^2\left(C_0\left({\displaystyle \frac{z}{H}}\right)+\gamma C_1\left({\displaystyle \frac{z}{a}}\right)V_s+C_2\left({\displaystyle \frac{z}{a}}\right)V_s^2\right)$	Finite-size particles in a confined microchannel flow for Re < 20.	Asmolov et al., 2018

The asymptotic solutions based on perturbation theories are convincing enough to explain the inertial microfluidic phenomenon in the most time-efficient manner. However, they are too weak to capture specific details, such as the particle-fluid interaction and disturbed fluid velocity profiles around the particle, and unable to apply in a complex 3D domain. What’s more, these methods have some restrictions including particle size, Re, and the distance between the particle and the wall, which limit their application in a variety of scenarios.

2.2. Direct numerical simulation methods

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 (Auton, 1987; Yang et al., 2005). However, their application to various scenarios have been limited by some restrictions and deficiencies, which are not comply with many complex cases of particle motion. As a powerful alternative, numerical method based on Navier-Stokes equations can determine the motion of particles for a wide range of sizes and Re. A robust method is direct numerical simulation (DNS), in which the forces acting on a particle are calculated from the interaction between the fluid and the particle. The most frequently methods adopt DNS for simulating particle motion can be divided into three significant classifications, which are briefly explained in the Table 2.

Table 2. Most common numerical methods for inertial particulate flows.

Methods		Applications	The advantages	The shortcomings	References
FSPP		Calculate the steady-state flow fields around a single particle at specific positions of the channel cross-section.	Highly efficient.	The solution of particle trajectories requires explicit formulas.	Carlo et al., 2009
ALE		Rigid particles in incompressible fluid flow.	The particle trajectories can be tracked at each time step.	Time-consuming.	Hughes et al., 1981
FDM	IBM	Apply to irregular and complex boundary conditions, and has greater geometric flexibility.	Simplicity and efficiency.	The accuracy of the solution near the immersed boundary is lower.	Glowinski et al., 1998
	DLM	Calculate the movement of various particles in both Newtonian and viscoelastic fluids.	The task of mesh generation is significantly simplified.	The size of the fluid mesh should be smaller than the particle mesh resulting in a large number of grids.	Xu et al., 2019

The first method to obtain the net inertial force acting on a particle at specific positions of the channel cross-section is Flow at Specific Particle Position (FSPP) (Carlo et al., 2009). The solution process consists two steps: calculating the pressure field and velocity field around the particle according to the continuity and Navier-Stokes equations, and then obtaining the values of the force and moment exerted on the particle. FSPP is highly efficient as it avoids problems encountered in the calculation of transient particle motion such as remeshing at each time step. But there are still some deficiencies in calculating particle trajectories, the extra explicit formulas (e.g. Newton’s second law) are required to obtain the linear and angular velocities of the particle.

In 1981, the Arbitrary Lagrangian-Eulerian (ALE) method was proposed for the system of suspended rigid particles in an incompressible fluid flow (Hughes et al., 1981). This method can track the particle trajectories at each time step by an iterative way. At the first time step, using incompressible continuity and Navier-Stokes equations to solve the fluid flow. At the next time step, integrating the total stress on the surface of each particle and Newton’s second law to update the particle positions, orientations, and their velocities. Meanwhile, owing to the mesh excludes the domain inside particles, it needs to be updated for recalculating the flow field at each time step, which is an inconvenient and time-consuming progress.

Compared to the ALE method, a fixed mesh used in Fictitious Domain Methods (FDM) (Glowinski et al., 1998; Prosperetti & Tryggvason, 2009; Xu et al., 2019) for the simulation of particulate flows includes the domain inside the particles. Therefore, a force is exerted on the domain to ensure that the fluid within each particle obeys the rigid body motion. According to the type of the applied force, FDM can be divided into the Immersed Boundary Method (IBM) and the Distributed Lagrange Multiplier (DLM) method. In IBM, the force is applied to the surface of the particle, whereas it’s applied to the body of the particle in DLM. The fixed mesh used in IBM consists of two different computational meshes, an Eulerian grid for the fluid and a Lagrangian mesh for particles. Consequently, the value of force and velocity between two meshes should be transferred using interpolation. By introducing the framework of the finite element method (FEM) into IBM, DLM was formed to deal with the movement of the particles. The flow field is solved over the whole computational domain, including the particles’ interior, where the motion of particle is enforced as rigid-body motion by using a distributed Lagrange multiplier. The method has been widely used in both Newtonian and viscoelastic fluids in 2D and 3D flow geometries.

All of the above DNS methods have been developed and applied to deal with the movement of the particles in various fluids. But the usability and validity of them differ in various microchannel structure, fluids and particles. Thanks to the calculation steps of flow field and particle trajectories are independent, FSPP is highly efficient and adaptable for complex microchannels, such as serpentine and spiral microchannel. However, the method is barely used for investigating the particles motion in non- Newtonian fluid due to the flow boundary condition is normally set to laminar flow. Besides, the particle size is constant in proposed formula for the calculation of inertial lift, which restricts its application in deformable particles. As for ALE method, it’s a technique based on a combined formulation of the fluid and particle momentum equations, together with an arbitrary Lagrangian–Eulerian moving unstructured finite element mesh technique. The method with a high applicability and accuracy has been used to solve particle motions in both Newtonian and viscoelastic fluids. It also handles particles of different sizes, shapes, and materials. To track the free movement of particles in the fluid, the domain occupied by the fluid is irregular. And the mesh needs to be updated for recalculating the flow field at each time step to track the particles trajectories with time. The update of mesh at each time step leads to an inconvenient and time-consuming progress. Therefore, ALE method is not suitable for particle motions in complex microchannel, e.g. serpentine and spiral microchannel. As stated in the ALE method, a time-dependent geometrically complex domain is used and re-meshed at each time. In FDM, it is replaced by a stationary fictitious domain, which is larger but simpler. Compared with ALE method, FDM provides mathematically simple and computationally efficient algorithms at the expense of accuracy and effectiveness. It is also a convenient approach to solve a great variety of particles motion in various fluids. IBM, as the one of the most widely used FDM, is well-suited for elastic particles and non-spherical particles. And the DLM, which is initially used to simulate the fluid–particle systems, has been extended to handle non-Newtonian fluids and massive particles problems, showing excellent performance.

2.3. Lagrangian particle tracking

In spite of being the only means to accurately obtain the lift forces, DNS often consumes high computational cost and even becomes impractical when applied to complex geometries (e.g. long serpentine and spiral channels) encountered in real-world microfluidic devices. To reduce the computing load, the results of asymptotic solutions are implemented in the LPT method to realize fast prediction of particle trajectories in complex microchannels for inertial microfluidic applications (Liu et al., 2016; Rasooli & Çetin, 2018; Zhou et al., 2017). This method can obtain a more realistic picture of real applications according to the assignment of statical distributions for the starting point, size and particle properties, which facilitates the design of inertial microfluidic devices. However, the calculation of asymptotic solutions requires the researchers with strong knowledge of computational fluid dynamics. Moreover, if one term of the formula generalized by DNS data, the implementation of these results for different problems requires a look-up table generated from the computationally heavy DNS (Liu et al., 2016).

3. Applications of numerical methods

3.1. Particle motion in Newtonian fluid flow

Viscosity of Newtonian fluid is independent of the shear rate, e.g. water, urine and cerebrospinal fluid. The particle in a Newtonian fluid experiences the inertial forces, which decide the particle motion. Feng et al. (Feng et al., 1994; Patankar et al., 2001) made the first attempts to investigate particle motion in Newtonian fluid flow numerically by using ALE method. These studies illustrated that the density difference between particle and fluid played a crucial role in the lateral equilibrium position of particles. Besides ALE method, Yang et al. (Yang et al., 2005) proposed lift law for circular and spherical particles by using DLM method, which was analogous to the aerodynamics lift law in specific case. However, ALE method was a better choice to obtain more precise results than DLM due to the investigation was about the motion of a single particle. Using the DLM method with some modifications, Shao et al. (Shao et al., 2008) found that the equilibrium position shifted from the Segre and Silberberg equilibrium position to the inner radius equilibrium position by increasing Re, as shown in Figure 1A. Compared with above tubular flows, the physics of inertial microfluidics in channels with rectangular cross-sections is more complicated, and the equilibrium position of particles is related to the cross-sectional shape. In this channels, the most important force acting on the particles is inertial lift force ($F_L$), which is a function of particle size ($a$), the aspect ratio of the cross-section ($AR$), Re, and the position of the particle in the channel cross-section. For investigating the migration of particles in straight rectangular (Liu et al., 2015; Mashhadian & Shamloo, 2019) and square (Carlo et al., 2009) channels, the FSPP method was used to calculate the value of $F_L$ at different positions within the cross-section (Figure 1B). By using numerical simulation with a combination of the Discrete Element Method (DEM) and the IBM, the migration of particles is obtained as Figure 1C (Carlo et al., 2009; Udono & Sakai, 2017; Wang et al., 2017). It is obviously that the equilibrium position of a constant-size particle moves towards the wall with the increase of Re. At a constant Re, conversely, the equilibrium position of particles approaches the channel centre with increasing the particle size (Figure 1D). Similarly, the equilibrium positions of particles in rectangular channels depend on Re and blockage ratio (k) (Liu et al., 2015; Mashhadian & Shamloo, 2019). By increasing Re, the equilibrium positions near the centre of short walls become stable in addition to the centre of long walls (Figure 1E). Moreover, if k decreases to lower value at a constant Re value, new equilibrium positions will emerge near the short walls.

Figure 1. Inertial microfluidics in Newtonian fluid flow. (A) The equilibrium positions in tube flow, square cross-section channel, and rectangular cross-section channel. (B) Inertial lift forces in a square cross-section channel. Reproduced from Ref. (Carlo et al., 2009) with permission from The Royal Society of Chemistry. (C) I. Snapshots of spatial distributions of particles in microchannels for AR = 1 at different Re. II. Lateral positions of focused particles as a function of Re. Reproduced from Ref. (Udono & Sakai, 2017) with permission from Springer. (D) Transverse distributions of the inertial lift force coefficients (C_FL) of particles with different diameters, $\lambda$ is ratio of particle diameter and channel width. Reproduced from Ref. (Wang et al., 2017) with permission from MDPI. E) The inertial lift distribution (upper) and fluorescence images of particle distribution (bottom) in rectangular microchannels at I. Re = 100, k = 0.3, II. Re = 200, k = 0.3 and III. Re = 100, k = 0.1. The nondimensional coordinates 2z/H and 2y/H are twice the ratio of short axes coordinate (z) to the channel height (H) and twice the ratio of long main axis coordinate (y) to the channel height (H), respectively. Reproduced from Ref. (Liu et al., 2015) with permission from The Royal Society of Chemistry.

To improve the particle focusing, the channels with curved structure or non-rectangular cross-section have been proposed. As the complexity of channels used for particle manipulation increases, using straight channels to investigate the physics of inertial microfluidics is no longer satisfactory, because particle focusing is highly related to the channel shape and cross-section. In the existing numerical studies, the calculation of particles’ motion in non-straight channels divided into two steps. First, using FSPP method to calculate the distribution of inertial lift force in the cross-section of the channel. Second, by combining inertial lift force with other existing forces, the particle trajectories can be obtained using the explicit formula (Eq. (1)) (Liu et al., 2015). (1) $mp\frac{d^2{x}_p}{d{x}^2}=F_{Inertialforce}+F_{Drag}+F_{Virtualmass}$(1) Apart from inertial lift force, particles in serpentine channels experience extra dean drag force. If the velocity of a particle exceeds the threshold, it will focus at the centre of the long walls of the serpentine channel, which is shown clearly in Figure 2A (Liu, 2016; Shamloo & Mashhadian, 2018). By changing the channel rectangle cross-section to triangle, the focusing positions of the particles in a straight channel can be controlled and strengthened (Figure 2B) (Kim et al., 2016, Kim et al., 2018; Mashhadian & Shamloo, 2019;). Similarly, the equilibrium positions of particles in any arbitrary channel cross-section can be obtained through the vector plot of the total force by combining the inertial force with the drag force (Figure 2C) (Kommajosula et al., 2019). Owing to the additional dean drag force, spiral channels have been frequently used for the separation of bio-particles and cells (Condina et al., 2019; Garcia & Pennathur, 2019; Ghadami et al., 2017a; Kwon et al., 2017; Rasooli & Çetin, 2018; Rzhevskiy et al., 2020). The effect of the drag force is dominant for particles with a size comparable to channels dimensions $a_c/h<0.07$, and traps particles in vortex streamlines (Figure 2D)

Figure 2. Inertial microfluidics in non-straight channels. (A) Particles move towards the equilibrium position after passing each U-turn. Reproduced from Ref. (Shamloo & Mashhadian, 2018) with the permission of AIP Publishing. (B) The magnitude of lift forces for triangular microchannel. Reproduced from Ref. (Kommajosula et al., 2019) with permission from Springer. (C) Statistics of particle positions with R_P (total N = 3000) in I. half-circular channel, II. wide triangular channel and III. narrow triangular channel. $R_{\text{p}}=\text{Re}{\left(\frac{a}{H}\right)}^2$ is particle Reynolds number. The normalized count (P) is shown for the top view and the side view. The y/W and z/H are the ratio of y-coordinates to the channel weight (W) and the ratio of z-coordinates to the channel height (H), respectively. Reproduced from (Kim et al., 2016) with the permission of The Royal Society of Chemistry. (D) I. Schematic illustration of the spiral channel in research. II. The cross-sectional trajectories the three particles in this spiral channel at Re = 100. The y/W and z/W are the ratio of y- and z-coordinates to the channel weight (W). The a/D_h is the ratio of particle diameter (a) to hydraulic diameter of the channel (D_h). Reproduced from Ref. (Garcia & Pennathur, 2019) with the permission of Springer. (E) The cell trajectories in the microchannel at Re = 100. The colour contours indicate the axial flow velocity at different cross sections. Reproduced from Ref. (Shiriny & Bayareh, 2021) with permission from Elsevier. F) The schematic representation of the hybrid spiral microfluidic platform and the numerical results of the cell separation study. Reproduced from Ref. (Altay et al., 2022) with permission from MDPI. G) The design of microchips with I. corrugated dragonfly wing patterns and II. with cilia walls and the simulation of particle trajectory in these microchips. Reproduced from Ref. (Yang et al., 2022) with permission from MDPI.

In conclusion, the particle motion in Newtonian fluid flow can be divided into two categories according to the channel structure, and the numerical methods for investigating each case are clearly shown in Table 3. The numerical researches have mainly focused on the particle motion within straight channel over the last two decades, presenting obvious diversity and maturity. Compared to straight channels, non-straight channels have a shorter development time in particle separation, resulting in a lack of numerical investigations, both in quantity and category. Nevertheless, the non-straight channels not only provide a compact design, but also allow for more exquisite control of particle migration. Therefore, the increasing separation devices with non-straight channels for real particles have been produced in recent 5 years. Correspondingly, numerical investigations for particle migration within non-straight channels attract various researchers, and will usher in an era of explosive growth.

Table 3. Numerical methods for investigating particle motion in different channels.

Channel structure		Numerical methods	References
straight	tube	ALE	Feng et al., 1994; Patankar et al., 2001; Joseph et al., 2002
		DLM	Shao et al., 2008; Yang et al., 2005
	rectangular cross-section	FSPP	Liu et al., 2015; Mashhadian & Shamloo, 2019; Carlo et al., 2009
		DEM + IBM	Udono & Sakai, 2017; Wang et al., 2017
	non-rectangular cross-section	FEM	Kim et al., 2016; Kommajosula et al., 2019
non-straight	serpentine channels	FSPP	Shamloo & Mashhadian, 2018; Liu et al., 2016
	spiral channels	FSPP	Garcia & Pennathur, 2019

Real-time and fast trapping and tagging of microfeatures, such as microparticles and cells, are of great significance for biomedical researches. Inertial technique for isolating CTC from whole blood serves as a valuable assistance to accelerate cancer diagnosis and guide the selection of treatment medications. In recent two years, increasing numerical methods have been used for investigating the inertial separation of real particles (e.g. CTCs and bacteria). Using FEM, Afshin et al. (Shiriny & Bayareh, 2021) investigated the inertial separation of two kinds of CTCs from blood cells in a spiral microchannel consisting of only one loop with a rectangular cross-section (Figure 2E). The results showed that the number and distribution of cells at three outlets were related to Re, and the separation efficiency reached 100% over a range of Re from 90 to 110. Dariush et al. (Bahrami & Bayareh, 2022) designed a novel spiral micromixer with sinusoidal channel walls to enhance the mixing index. To quantify the effects of Re, the amplitude, the wavelength, and the number of loops on the performance of the micromixer, a set of numerical simulations were performed by using the finite-difference method. For particle/CTC separation, Rana et al. (Altay et al., 2022) proposed a multi-stage hybrid microfluidics platform utilizing inertial and acoustic radiation forces. Based on the FEM, the simulation of passive separation of particle in the spiral microchannel was performed using the Laminar Flow and Particle Tracing Module of the software COMSOL Multiphysics 5.5 (Figure 2F). Moreover, the active separation induced by pressure distribution due to the applied SAWs, was also simulated using the Pressure Acoustics, Frequency Domain, and Particle Tracing Module. To improve the classical straight-line channel, Yang et al. (Yang et al., 2022) set the dragonfly corrugated wing structure and the cilia wall structure inside the two stages microchannels (Figure 2G). Using Fluid-Particle Interaction Module of COMSOL Multiphysics, the vortex regimes, boundary layers, and particle migration in microchannels were comprehensively revealed at a varying Re flow regimes. The computational fluid dynamics (CFD) of particle flow through the corrugated dragonfly wing microchip and the cilia wall microchip are shown in Figure 2G. Compared with the corrugated dragonfly wing microchip, more micro vortices were induced in microchannel by the irregular cilia, which also swept fluids and foreign particles from one side to the other very easily. The research found that the cilia wall microchip with very small dimensions could do the same job rapidly and efficiently with a lower flow rate of the fluid. Besides dispersed particles, the self-ordering and organization of an in-line particle chain flowing through a square microchannel were also studied by Liu et al. (Liu & Pan, 2022). The results obtained by the Immersed Boundary-Lattice Boltzmann method implied that the migration of small-particle chains was relatively complex considering the effects of Re, length fraction, characterizing particle concentration, and particle size. Owing to the subtle interaction among vortex-induced repulsive force, wall-induced lift force and shear gradient lift force, migratory patterns are categorized as Stable Pattern, Spring Pattern, and Chaotic Pattern.

3.2. Particle motion in non-Newtonian fluids

In contrast to the particle motion in Newtonian fluid (Joseph & Ocando, 2002), the particles in a non-Newtonian fluid experience several lateral forces in addition to the inertial forces, depending on the non-Newtonian fluid’s rheological properties. Therefore, the analysis of particle motion in a non-Newtonian fluid is more complicated. As a new field, the recent major studies have preferentially focused on particle motion in viscoelastic fluids. In order to obtain the flow field in a viscoelastic fluid, continuity and momentum equations need to be solved (Eqs. (2) and (3)). Besides, considering the effect of viscoelasticity, an extra stress tensor ($\tau$) must be added to the total stress tensor (Eq. (4)) (Beris et al., 1992). (2) $\begin{array}{rl}& \mathrm{\nabla }.u=0\end{array}$(2) (3) $\begin{array}{rl}& \rho \left(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+\left(u.\mathrm{\nabla }\right)u\right)=\mathrm{\nabla }.\sigma \end{array}$(3) (4) $\begin{array}{rl}& \sigma =\text{ - pl}+2{\eta }_{s}e\left(u\right)+\tau \end{array}$(4) where ${\eta }_s$ and $e\left(u\right)$ are the Newtonian solvent viscosity and strain rate tensor, respectively. To solve this set of equations, the relation between $\tau$ and other parameters should be obtained firstly by a constitutive equation (Keunings, 2003; Larson, 2013). Two frequently used constitutive laws are Giesekus (Eq. (5)) and Oldroyd-B equations (Eq. (6)) (Avino & Maffettone, 2015). Where $\alpha$, ${\eta }_p$ and $\lambda$ are the mobility factor, polymeric viscosity, and relaxation time, respectively. (5) $\begin{array}{rl}& \tau +\lambda \left(\frac{\delta \tau }{\delta t}+u.\mathrm{\nabla }\tau -\mathrm{\nabla }{u}^T.\tau -\tau .\mathrm{\nabla }u+\frac{\alpha }{{\eta }_p}\left(\tau ,\tau \right)\right)\\ & ={\eta }_p\left(\mathrm{\nabla }u+\mathrm{\nabla }{u}^T\right)\end{array}$(5) (6) $\begin{array}{rl}& \tau +\lambda \left(\frac{\delta \tau }{\delta t}+u.\mathrm{\nabla }\tau -\mathrm{\nabla }{u}^T.\tau -\tau .\mathrm{\nabla }u\right)\\ & ={\eta }_p\left(\mathrm{\nabla }u+\mathrm{\nabla }{u}^T\right)\end{array}$(6) In viscoelastic fluids, the particles experience elastic force and inertial force, and the competition of them decides the migration and focusing pattern of the particles in the channel cross-section, as shown in Figure 3A (Raffiee et al., 2019; Raoufi et al., 2019; Trofa et al., 2015; Yu et al., 2019). Unlike Newtonian fluids, the centre of channel is the unique equilibrium position of particle in viscoelastic fluids. Based on the DLM method in the framework of the finite volume method on a staggered grid, Li et al. (Li et al., 2015) investigated the effects of fluid elasticity, fluid inertia, and shear-thinning viscosity. Research revealed that particles in Oldroyd-B fluid (Figure 3B I) focused at the centre of the channel, while reached a position away from the centre in Giesekus fluid (Figure 3B II).

Figure 3. Inertial microfluidics in non-Newtonian fluid flow. A) Migration of particles in a straight channel filled with the viscoelastic fluid is different from the Newtonian fluid flow. Reprinted from (Raoufi et al., 2019), with the permission of AIP Publishing. (B) Steady flow field around the particle in a channel filled with I. Oldroyd-B and II. Giesekus fluids. Reprinted from (Li et al., 2015), with the permission of Cambridge University Press. (C) At various Wi and Re, distance of the off-centre equilibrium points of particles from the channel centre along the I. main axis and II. diagonal axis. The filled and open symbols indicate stable and unstable equilibrium points, respectively. The d/w is the ratio of the distance to the half the width of channel (w). Reprinted from (Raffiee et al. 2019) with the permission of Elsevier. D) Phase diagram for the equilibrium positions as a function of El and Re at a blockage ratio of 0.15. Reprinted from (Yu et al., 2019), with the permission of Cambridge University Press. (E) Temporal evolution of cell distance from the centreline and the distribution of cells in channel. Reprinted from (Raffiee et al., 2019), with the permission of Springer. (F) Particle trajectories for different Wi and the existence of possible multiple particle equilibrium states in very dilute pressure-driven flows. Reprinted from (Banerjee et al., 2019), with the permission of Springer.

Recently, a great quantity of investigators focused on the particles motion in elasto-inertial flows. Raffiee et al. (Yu et al., 2019) performed 3D computational investigation of a particle in viscoelastic channel flow considering inertial and elastic effects. It was found that the fluid elasticity affected the particle dynamics indirectly by changing the velocity field. And the equilibrium points of particle for the entire range of Re and Wi numbers are illustrated in Figure 3C. According to a modified version of the DLM method, Yu et al. (Raoufi et al., 2019) investigated focusing patterns of particles in Oldroyd-B fluids with k of 0.15 under various conditions. The Figure 3D clearly shows that the equilibrium position (EP) depends strongly on the elasticity number (El, El = Wi/Re) and weakly on the Re. When El = 0.001–0.15, the Diagonal EP occurs largely at Re = 10–100. When El > 0.015, the equilibrium positions are always located on the channel centreline. Raoufi et al. (Raffiee et al., 2019) used FSPP to investigate the effect of channel cross-section on the focusing efficiency of elasto-inertial flow. Research found that the elastic force pushed the particles towards the centre of the cross-section more efficiently by increasing the channel corner angle. Using the front-tracking method, Raffiee et al. (Raffiee et al., 2019b) simulated the motion of cells within a square microchannel in viscoelastic fluid, and the evolution of the average cell distance from the centreline (r*) was shown in Figure 3E. It can be seen that r* initially decreases and reaches a quasi-steady state, and it decreases with increasing Weissenberg (W_i/We) number. The results also revealed that increasing the solution concentration significantly led to the growth of volumetric flow rate, which was beneficial for boosting the total throughput of a microfluidic device. Using direct numerical simulations, Bannerjee et al. (Banerjee et al., 2021) studied the phenomenon of particle elasto-inertial focusing in presence of moderate inertia and weak fluid elasticity, and found a general tendency of the particle to assume multiple equilibrium positions based on competitive effects of inertia and elasticity. An IBM method was used to account for the presence of the particles and a FENE-P model was used to simulate the presence of polymers in a non-Newtonian fluid. For the first time, the results demonstrated the possibility of the analog tuning of the particle focusing positions in elasto-inertial flows. Using a single-loop microdevice, Afshin et al. (Shiriny et al., 2022) investigated the inertial separation of microparticles suspended in shear-thinning fluids. Numerical simulations were performed for focusing the trajectory and the separation rate of polystyrene microparticles with the impact of parameters, such as diameters, power-law index and relaxation time constant. Owing to the shear-thinning effect, the centre of the Dean vertices deviated from the symmetrical state and moved towards the outer wall of the microchannel. And the separation efficiency decreased by reducing the power-law index and relaxation time constant.

Seeing as this field is a subject of ongoing investigations, major studies focus on the particles motion in viscoelastic fluids with different rheological properties. Compared with channel diversity of Newtonian fluids, recent researches for non-Newtonian fluids almost selected straight channel as object, while ignoring the effect of channels structure on the separation of particles in practical applications. Moreover, the real particles are rarely used in above numerical studies for particles motion in non-Newtonian fluids, which isn’t content with requirements of leading-edge application in fast real-time trapping and tagging of cells. Therefore, there are still a vast gap between numerical researches and practical applications, especially in terms of complex channels structure and real particles aspects.

3.3. Effect of the particle on fluid flow

In order to trace the particle migration, the majority of studies focus on the effect of fluid ignoring 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 (Collins et al., 2016; Kulrattanarak et al., 2008; Manz et al., 1990). In general, particles significantly perturb the fluid around them, causing reversed streamlines or secondary flows in different boundary shear flow. In an inertia-less open boundary shear flow (i.e. Stokes flow), two kinds of streamlines could be obtained via numerical simulations using the front-tracking finite difference method, namely closed and open streamlines (Singh & Sarkar, 2011). They are depicted in Figure 4A, and the closed streamlines spiral around the particles near the vorticity axis. By adding fluid inertia to the system, two reversed streamlines are created owing to the collapse of the closed streamlines around the particles, as clearly shown in Figure 4B. Besides, there are open spiralling streamlines outside the closed streamlines. Both the spiralling streamline and reversed streamline areas meet each other at a stagnation point. In confined inertial flows, the reversed streamline can also be found by using FSPP method (Figure 4C), which is similar to that in open boundary flows (Lee et al., 2010). However, owing to the parabolic nature of confined inertial flows, the symmetry streamlines area in the vertical direction doesn’t exist. Except for reversed streamlines, the net secondary flows were also induced by particles in finite-Reynolds number flow. And the presence of reversed streamlines in confined inertial flows has been experimentally demonstrated for very low Re as well as finite Re flows, as shown in Figure 4D. Besides the patterns of streamlines along the fluid direction, the streamlines around a particle in channel cross-section can be obtained using FSPP method, as shown in Figure 4E. We can see that the secondary flows are generated once there is a particle in the channel, even the particle is located at the channel centre without rotation (Amini et al., 2012). The direction and magnitude of the secondary flows depend heavily upon fluid inertia, the location of the particle in the channel, particle rotation rate, etc.

Figure 4. The effects of the particle on fluid flow. (A) Distribution of streamlines around a particle in an open boundary shear flow. Reprinted from (Singh & Sarkar, 2011), with the permission of Cambridge Core. (B) Reversed streamline areas around the particle in the inertial flow. Reprinted from (Singh & Sarkar, 2011), with the permission of Cambridge Core. (C) Reversed streamline observed in both I. Stokes flow and II. finite-Reynolds number flow. Reprinted from (Lee et al., 2010), with the permission of PNAS. (D) The net secondary flows were induced by particles in confined channels with finite inertia. Reprinted from (Lee et al., 2010), with the permission of PNAS. (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. σ is a nondimensional measure of the net secondary flow strength. Reprinted from (Amini et al., 2012) with the permission of PNAS. F) The evolution of vorticity contours and particle distribution at three different times and conditions. Reprinted from (Lin et al., 2022) with the permission of AIP Publishing.

Recently, the effects of the particle on turbulent flow have also been investigated. The presence of solid particles or water droplets in continuous fluid flow can induce turbulence attenuation or amplification, which depends on the characteristics of the particles (e.g. volume fraction, mean diameter and mass density) and carrier phase flow properties (Gai et al., 2020). Using a localized artificial diffusivity scheme and an Eulerian–Lagrangian approach, Xia et al. (Xia et al., 2016) simulated two-way interactions between heavy point particles and forced compressible homogenous turbulence. The results showed that the turbulence was suppressed more obviously by increasing the density of particles. In order to improve the efficiency of real particles sorting, more and more complex inertial microfluidic structures were proposed, such as combining lateral structural columns, adding grooves and herringbones. Therefore, it is necessary to numerically simulate the flow of fluid containing particles around obstacle. Lin et al.(Lin et al., 2022) developed a two-way fluid-solid coupling model by using the DEM and IBM to numerically study the effect of particles on the dynamics of fluid flows around an obstacle. The results showed that the interaction between fluid and particles caused the original regular vortices to collapse behind the obstacle, as shown in Figure 4F. As the particle mass loading ratio (β) increases, the variation of vortex intensity and shedding frequency become more complicated, and the collapse of the vortices becomes more obvious.

3.4. Non-spherical particles

In majority of the theoretical, numerical, and experimental works conducted on inertial migration, particles are commonly assumed to be rigid and spherical particles. However, particles found in real suspensions are rarely spherical, and the shape of a particle plays a vital role in the inertial particle motion. According to the shape, the non-spherical particles can be classified into three main categories: ellipsoids, non-ellipsoidal axisymmetric and asymmetrical types (Figure 5A). The transport of non-spherical particles in a microchannel is indeed more complicated than the spherical ones, since their translational and rotational motions are strongly coupled. The rotational behaviour of ellipsoidal particle is extremely complicated, the effects including fluid inertia, buoyancy, the aspect ratio and the initial orientation of particle. Jeffery’s theory is the first one that explains the rotational behaviour of an ellipsoid in general shear flow. In a simple shear flow, the ellipsoid rotates around the vorticity axis along the Jeffery orbits, as shown in Figure 5B (Tohme et al., 2021). The rotation modes can be divided into a conventional rotation mode and two extreme motion modes, namely ‘kayaking’, ‘tumbling’ and ‘log-rolling’, which are indicated by red dotted line, blue dotted line and green curve with arrow respectively. Since the proposal of this theory, more researchers have investigated the motion of non-spherical particles in fluid using numerical method, such as IBM and LBM. Lashgari et al. (Lashgari et al., 2017) were the first to investigate the entire migration dynamics of oblate particles in a duct from moderate to high Re. The diagrams of oblate particles migration obtained by IBM method are clearly shown in Figure 5C. The figure consists of three charts, the left chart is the lateral trajectories of a rigid oblate particle in a square duct at Re_b = 100 (blue line), Re_b = 150 (green line), Re_b = 200 (red line), the chart on the top right is zoom of a particle trajectory at Reb = 200, and the chart on the bottom right is the evolution of the x-component of the orientation vector (n_x) with time at various Re (same colour coding as in left chart). For Ri_b = 100 and Ri_b = 200, the oblate spheroids move laterally towards the equilibrium position, and their final orientation is stable. And the focusing length of the oblate spheroid at Ri_b = 150 is shorter than that at Ri_b = 100, which indicates that the reduction in the focusing length by increasing the Reynolds number is also evident for the oblate particles. However, when the Ri_b = 200, the oblate particle exhibits time-dependent rolling and tumbling motions around an equilibrium position at the face centre. Therefore, there isn’t a steady orientation and rotation-rate at Ri_b  = 200, and the zoom on the particle trajectory shown in top right proves its chaotic dynamic. Using LBM, Nizkaya et al. (Nizkaya et al., 2020) studied the inertial focusing of oblate spheroidal particles in channel flow at moderate Re. The research found that the oblate spheroids eventually reoriented to the stable log-rolling motion around the axis of symmetry, as shown in Figure 5D I. The Figure 5D II show that the equilibrium positions for all spheroids with different b/a values are very close. But the speed of migration to the equilibrium position is related to b/a, which is decreased with the decrease of the b/a. Rosén et al. (Rosén et al., 2015) used the LBM with an external boundary force to describe the effect of inertia on the rotational behaviour of an oblate spheroidal particle, which was confined between two parallel opposite moving walls. Using the same method, Miao et al. (Miao & Xiao, 2021) studied the behaviour of a neutrally buoyant ellipsoidal particle in vortical flow confined by a microcavity. To summarize, the conclusions converge on the fact that prolate spheroid tends to enter a tumbling mode whereas an oblate spheroid is apt to achieve a rolling mode. However, at higher value of Re, the results are divergent, and different rotational modes are observed.

Figure 5. Migration of non-spherical particles in confined flows. (A) Three kinds of non-spherical particles. (B) Components of the rotation axis of a rod-like particle in a rectangular channel. (C) The orientation of an oblate particle at the equilibrium position for different Re. Open circles and triangles show the initial and final equilibrium position, respectively. And the solid (for stable cases) and dashed (for unstable cases) lines show the final orientation of the oblate spheroid. Reproduced from Ref. (Lashgari et al., 2017) with permission from ResearchGate. (D) I. An oblate spheroid orienting in a pressure-driven flow to perform a stable log-rolling state. II. trajectories for spheroids with a/H = 0.15 and b/a = 0.5(solid curve), 0.8 (dashed curve), and 1 (dotted curve). Reproduced from Ref. (Nizkaya et al., 2020) with permission from AIP Publishing. (E) The motion of rigid rod-like particles with different aspect ratios in rectangular channels. I. The schematic picture and II. the image of multiple overlays of frames captured at the channel. Reproduced from Ref. (Masaeli et al., 2012) with permission from The Royal Society of Chemistry. (F) Lateral position of a biconcave particle plotted versus time when moving in the midplane of the channel. La is deformability of particle. Reproduced from Ref. (Patel & Stark, 2021) with permission from The Royal Society of Chemistry. (G) Instantaneous distribution of Str = 5 fibres in turbulent channel flow. Reproduced from Ref. (Dotto et al., 2020) with permission from The Royal Society of Chemistry.

Besides ellipsoids, the rotation of axisymmetric non-ellipsoidal particle in flow has received growing attention in the last decade. Using FSPP, Masaeli et al. (Masaeli et al., 2012) studied the motion and dynamic equilibrium position of rigid rod-like particles with different aspect ratios in rectangular channels. The rotation modes of particles are classified as in-plane rotation, out-of-plane rotation and no rotation. By increasing Re and the aspect ratio of particles, particles mainly exhibit in-plane rotation and the equilibrium position approach to the centre of cross-section, which can be seen in Figure 5E. This phenomenon has been used to separate particles with different rotational diameters. Kuntal et al. (Patel & Stark, 2021) performed 3D lattice Boltzmann simulations combined with the IBM to unravel the dynamics behaviour of two-particle (e.g. mono- and bi-dispersed pairs) in inertial microfluidics. As shown in Figure 5F, the stability of biconcave particle depends on the particle’s ability and the initial position. The mono-dispersed pair with soft capsules is always stable irrespective of the starting position, whereas the bi-dispersed pair performs bounded oscillation.

On a higher level of complexity, few works investigated the effect of the particle asymmetry on its rotational behaviour. The motion of asymmetric particles has no fixed period of rotation around the vorticity axis of fluid, which is qualitatively different from axisymmetric particle. Recently, the slender and asymmetric fibres in shear flow were numerically investigated (Thorp & Lister, 2019). Using an Eulerian-Lagrangian approach, Dotto et al. (Dotto et al., 2020 ; Dotto & Marchioli, 2019) developed a comprehensive database to investigate the orientation, distribution and deformation of inertial flexible fibres within turbulent channel flow. The results revealed that flexible fibres exhibited a stronger tendency to accumulate in the very-near-wall region owing to inertia-driven mechanisms namely turbophoresis, as shown in Figure 5G. And the bending of flexible fibers with small inertia is enhanced by mean shear and turbulent Reynolds stresses, while being reduced for fibres with large inertia, indicating that the fibre deformation is affected by flow anisotropy and fibre inertia.

3.5. Deformable particles

Except for non-spherical, real particles (e.g. cells and microgels) have another significance property is deformability, which leads to the equilibrium position of particle approaches the channel centre (Figure 6A) (Hadikhani et al., 2018). Due to deformable particles exhibit tank-treat motion while the particle is rotating, they will reach a steady-state shape and migrate towards their equilibrium position (Figure 6B) (Doddi & Bagchi, 2008). Deformable particles can be divided into three categories, droplet (Chen et al., 2014; Hadikhani et al., 2018), deformable capsule (Ebrahimi et al., 2021) and elastic particle (Villone & Maffettone, 2019). Droplet and bubble are common real particles with deformability, which have been extensively researched (Cannon et al., 2021). The Figure 6C shows the evolution of the y position of a droplet released from different positions of rectangular channels at various situations (Chen et al., 2014). For every Re, droplets migrate to two equilibrium y positions between the centreline and the wall, namely inner equilibrium position (y_eⁱ) and outer equilibrium position (y_e^o). And the equilibrium positions are decreased with increasing the Re and We. Unlike the behaviour of rigid particles, bubbles focus at the centre of the short walls by increasing the rectangular channel aspect ratio (Figure 6D) (Hadikhani et al., 2018). As for capsule membranes, several constitutive equations were employed to model the mechanical behaviour, namely the Hookean, the neo-Hookean law, and Gent models (Lac et al., 2004). Owing to the compatibility with experimental results of neo-Hookean law, it has been used in most cases. Using a front-tracking method, Raffiee et al. (Raffiee et al., 2017) investigated the role of polymeric fluids on the deformable cell migration behaviour. The results showed that the equilibrium position of the cell was on the channel diagonal, in contrast to that of rigid particles. Saman et al. (Ebrahimi et al., 2021; Ebrahimi & Bagchi, 2021) studied the cross-streamline migration and the focusing of deformable capsules in 3D curved vessels using a numerical model based on IBM. The results suggested that the equilibrium location of capsule within finite inertia flow depended on both capsule deformability and Re. As shown in Figure 6E I, decreasing centre radius (R_C) of a torus vessel results in the equilibrium radial location (R’_eq) being increasingly closer to the inner edge of the torus. But Capsules migrate to the same radial location that is independent of capillary number (Ca), as shown in Figure 6E II. Using an immersed-boundary lattice Boltzmann method, Feng et al. (Feng et al., 2021) investigated the effects of capsule volume fraction and bending stiffness (Eb) on the rheology of the suspension in two-dimensional confined Poiseuille flow, which decided the overall lateral position of the capsules. The equilibrium lateral location and Taylor deformation of a single capsule are functions of Re and bending stiffness (E_b), as shown in Figure 6F. For a given Re, the equilibrium position of the capsule with softer capability (smaller E_b) is closer to the centreline. However, the capsules with larger E_b are less likely to deform, as shown in Figure 6F II.

Figure 6. Migration of deformable particles in confined flows. (A) Schematic of the important forces (e.g. wall-induced lift force, the shear gradient lift force and the deformation-induced lift force) acting on a deformable particle. Reproduced from Ref. (Hadikhani et al., 2018) with permission from The Royal Society of Chemistry. (B) Migration of a capsule in a pressure-driven channel flow. Reproduced from Ref. (Doddi & Bagchi, 2008) with permission from Elsevier. (C) Evolution of the y position of a droplet released from different positions in a microchannel at I. Re = 6.12, We = 1.5 and II. Re = 19.3, We = 15. Reproduced from Ref. (Chen et al., 2014) with the permission of AIP Publishing. (D) Effect of the channel geometry on the lateral position of the bubble while the Re and Ca are constant. Reproduced from Ref. (Hadikhani et al., 2018) with permission from The Royal Society of Chemistry. (E) The effect of I. various the torus geometry (Rc) and II. Ca on the equilibrium radial location ($R^{\mathrm{\prime }}$) of the capsule. Reproduced from Ref. (Ebrahimi et al., 2021) with permission from Cambridge Core. (F) I. Equilibrium lateral location and II. Taylor deformation of a single capsule as a function of Re and Eb. Reproduced from Ref. (Feng et al., 2021) with the permission of AIP Publishing.

As for elastic particles, they are good models of biological cells for exploiting particle mechanical-property-dependent flow behaviour to diagnose illnesses. Villone et al. (Villone & Maffettone, 2019) reviewed papers that deal with the flow dynamics of individual elastic particles, as well as the applications involving soft particles. 1950, Roscoe predicted the dynamic behaviour of an elastic particle in unbounded shear flow, and proposed that the trends of its steady-state deformation parameter and orientation angle were related to elastic capillary number. Subsequently, this opinions have been validated and extended through several numerical techniques, namely the FEM (Villone et al., 2015), the polarization method (Gao et al., 2012), the Lattice-Boltzmann method (Wu & Aidun, 2010), the finite difference method (Sugiyama et al., 2011), and the boundary element method (Subramaniam & Gee, 2016). Recently, more researches have focused on the motion of elastic spheroid particles in confined shear flows. Using the ALE method, Villone et al. (Villone et al., 2015) investigated the motion of elastic spheroid particles in confined shear flows. The wall-induced force causes more particle deformation compared to that in an unbound shear flow (Figure 7A). Alghalibi et al. (Alghalibi et al., 2019) performed fully Eulerian numerical simulations of a hyperelastic particle suspended in a Newtonian pressure-driven flow to disentangle the interplay of inertia and elasticity on the particle focusing. The simulations verified that the migration dynamics and the final equilibrium position were strongly dependent on the particle elasticity. In particular, as the elasticity increases, the particles migrate faster and accordingly reach the final equilibrium position at the centreline in a shorter time (Figure 7B). To simulate the deformation of biological cells in arbitrary 3D flows (Figure 7C I), a hyperelastic Mooney–Rivlin model was incorporated into lattice Boltzmann simulations via an Immersed-Boundary algorithm (Müller et al., 2021). The Taylor deformation parameter D of cell is increased with the increase of capillary number Ca, as shown in Figure 7C II. The elastic behaviour of cells simulated by the numerical model is in excellent agreement with the measurement, which provides an improved set of tools to predict cell deformation.

Figure 7. Migration of elastic particles in flows. (A) Shape evolution of a neo-Hookean elastic prolate spheroid in Newtonian fluid. Reproduced from Ref. (Villone and Maffettone, 2019) with permission from APS Physics. (B) The particle radial position over time at various We. The radial position r is normalized by the radius of the pipe R whereas time is scaled by the ratio of pipe of diameter (D) to the bulk velocity across the domain (U_b) (D /U_b). Reproduced from Ref. (Alghalibi et al., 2019) with permission from APS Physics. (C) I. Simulation snapshots of the converged shapes of a single cell with different capillary number Ca in shear flow. II. The comparison of Taylor deformation parameter D at different capillary number Ca in different researches. Reproduced from Ref. (Müller et al., 2021) with permission from Springer.

4. Numerical simulation for microfluidic devices design

Based on numerical methods, more details of particles separation process in microchannels were presented by simulation data. A deeper insight into the particle migration characteristics, focusing mechanism, even channel deformation analysis can be promoted by comprehensive analysis of lateral force field maps and Dean vortex configurations, which is beneficial for the design and optimization of cell separation microfluidic devices. It is widely known that the geometric parameters, Re and cell characteristics are key factors of inertial particle separation in microfluid. Consequently, optimizing the structure of microfluid as a universal way to improve the performance of microfluidic devices has attracted extensive researchers. And various numerical methods were applied as an auxiliary means to evaluate the availability of innovative structure recently (Do et al., 2020; Ghadami et al., 2017b; Hafemann et al., 2020; Jiang et al., 2019; Palumbo et al., 2020).

4.1. Optimization of microchannel structure

Palumbo et al. (Palumbo et al., 2020) examined the effects of various geometric parameters (e.g. channel pitch, diameter, taper angle and depth) on inertial particle separation in a helical channel (Figure 8A) with trapezoidal cross sections using a calibrated numerical model. The simulation results were in good agreement with the experimental data under various Re and AR. This study is beneficial to the design and assessment of microfluidic devices for optimal size-based inertial particle separation in a quick way. To elucidate the mechanisms about the uneven distribution of particles over the cross section, Hafemann et al. (Hafemann et al., 2020) conducted a series simulations to geometrically resolve the particles in spiral ducts and couple them to the fluid using IBM. More new rich information on the details of separation process in continuous phase was provided by simulation data, including migration time, particle positioning in the cross section, streamwise particle spacing and velocity field. The focusing mode of particle is related to particle concentrations and the mean flow of the continuous phase (Figure 8B). Compare Figure 8B II with Figure 8B I, it can be seen that the particles interact vigorously via the perturbation they introduce in the flow field when the concentration is 15 times higher than before. As a result, unlike the propagation of individual particle along straight trajectories (the bottom left figure), the particles intermingle erratically in a crowd (the bottom right figure). Using a IBM numerical solver developed in the framework of OpenFOAM open-source software, Do et al. (Do et al., 2020) investigated the inertial focusing of a microparticle in spiral channels of rectangular and trapezoidal cross sections numerically for the first time (Figure 8C). In the channel of rectangular cross sections, there are two stable equilibrium position, which are near the inner channel. On the contrary, the stable equilibrium position of particle in trapezoidal channel is near the outer channel, and its number is only one. The comprehensive analysis of simulation data provides investigators a deeper insight into the focusing mechanism of particle in irregular channel. For improving the performance of spiral microchannel to meet potential clinical applications, Ghadami et al. (Ghadami et al., 2017b) introduced a novel spiral microchannel with stair-like cross section. The design fundamental concepts, design criteria and efficacy were investigated thoroughly using a robust numerical model and tested by experiments on the fabricated spiral microchannel. Based on LBM-IBM model, Jiang et al. (Jiang et al., 2019) studied the particle-focusing mechanics in contraction-expansion microfluidic channels through numerical simulation and validation experiments (Figure 8D). According to the simulation results, it could be found that a large γ continuous contraction-expansion microchannel required higher flow rate to ensure the effective separation of particles. The particle migration characteristics in channels with different contraction–expansion structure can provide a helpful insight for particle/cell detection chip design in the future.

Figure 8. Numerical method for the optimization of microchannel structure. (A) Flow contours and Dean flow vectors of the left reference channel and schematic of forces acting on the particles. Reproduced from Ref. (Palumbo et al., 2020) with the permission of AIP Publishing. (B) Particle positions and instantaneous flow field at conditions of I. MT100-80-2P and II. HC 100-80-22P. Reproduced from Ref. (Hafemann et al., 2020) with the permission of AIP Publishing. (C) Lateral force field in I. rectangular and II. trapezoidal cross-sectional spiral channels. The unstable (triangular makers) and stable (circular markers) equilibrium positions are determined from the intersection of the zero lines of radial (green) and vertical (red) force components. Reproduced from Ref. (Do et al., 2020) with the permission of AIP Publishing. (D) Simulation results of different diameter particle trajectories in the contraction-expansion microchannel. Reproduced from Ref. (Jiang et al., 2019) with permission from Springer. (E) A 3D profile of the overall channel deformation and the channel deformation at cross section of each loop. Reproduced from Ref. (Jeon et al., 2022) with permission from The Royal Society of Chemistry.

Except for the whole microchannels structure, the bifurcating (e.g. Y- and T-shaped junction) microchannels to trap the cells of a specific type were numerically studied by Hymel et al. (Hymel et al., 2019) Using a custom computational algorithm, the migration of separated deformable cell in Y-junction microchannel with different bifurcation angles, was systematic numerically analyzed to identify the effects (e.g. cell size, cytoplasmic viscoelasticity, flow rate and bifurcation angle, etc.) on the critical separation distance. The results showed that the Y-junction microchannel with a bifurcation angle of 120^o was most efficient and least dependent on the flow rate for trapping white blood cells and CTCs. This numerical study provides important design criteria to optimize microfluidic technology for deformability-based cell sorting and isolation. To overcome the limitations of conventional PDMS devices, Jeon et al. (Jeon et al., 2022) designed a deformation-free and mass-producible plastic spiral inertial microfluidic system. By visualizing particle trajectories and deformation profiles according to numerical analysis, the significant impacts of channel deformation on the device performance in PDMS devices were demonstrated. Based on channel deformation analysis aided by numerical simulation and confocal imaging (Figure 8E), the empirically optimized design of PDMS devices was performed for blood separation. The study also provides a systematic framework to transfer the accumulated library of inertial microfluidics design towards real-world industrial applications.

4.2. Assessment of complex microchannel

Since the widely applying of LPT, more realistic pictures of particle motion in complex microchannel in real applications were obtained. To understand and predict the separation mechanism of inertial focused diamagnetic particle (Figure 9A I), Zhou et al. (Zhou et al., 2017) developed a 3D numerical model based on LPT to simulate particle motions in ferrofluid flow through the entire microchannel. This study investigates the parameters (e.g. flow rate of ferrofluid and ferrofluid concentration) on the magnetic separation of a mixture of inertially focused polystyrene particles. As the magnetic force is a linear function of ferrofluid concentration, the magnetic deflection of each type of particles increases at a higher ferrofluid concentration as viewed from the superimposed images in Figure 6 (left and middle columns for 10 and 20 µm particles, respectively). And the numerical predictions shown in Figure 9A II (right column) are approximately consistent with the experimental observations. Owing to DNS often consumes high computational cost, Su et al. (Su et al., 2021) proposed a fast numerical algorithm in conjunction with machine learning techniques for the analysis and design of inertial microfluidic devices. The inertial lift distribution of microchannels was obtained by the machine-learning model. Then, it was integrated into the LTP method to quickly predict the particle trajectories in multistage spiral channel with rectangular cross section (Figure 9B I). The simulation results of particles trajectories with different diameters at four specific locations are shown in Figure 9B II. The particles were effectively separated based on size after through the channel, and the simulation results are good agreement with experimental observations. The database and the associated codes allow researchers to expedite the development of the inertial microfluidic devices for particle manipulation. For cancer cell separation, Nasiri et al. (Nasiri et al., 2022) carried out numerical modeling to study the fluid flow, particles’ trajectories in the serpentine inertial device before its fabrication. At the flow rate of 1000 µl/min, the best separation performance is obtained for the proposed inertial biochip for the separation of two particles with sizes of 10 and 15 µm. It is shown in Figure 9C, the streak width of the particles is decreasing by moving from the inlet of the channel to the outlet, and the CTCs are concentrated at the centre of channel, significantly away from other cells. Compared to the serpentine microchannels, a new zigzag microfluidics (Z-RISE) proposed by Bazaz et al. (Bazaz et al., 2022), demonstrates superior separation and purity efficiency due to the sudden channel cross-section expansion at the corners. Its separation performance is studied by experiments and numerical simulation based LPT method. The provided Z-RISE device and the numerical modeling allow the generation of a new class of solid microfluidic devices for efficient particle manipulation.

Figure 9. (A) I. Schematic diagram illustrating the mechanism of particle focusing and separation for the proposed ferrofluid-based microfluidic device. II. Inertially focused diamagnetic particle separation in EMG 408 ferrofluid of different concentrations at a constant flow rate of 1 ml/h. The black zones at the lower corner of the expansion are the buildup of ferrofluid. Reproduced from Ref. (Zhou et al., 2017) with the permission of Springer. (B) I. Schematic diagram of a multistage rectangular channel. II. The simulation results of particles trajectories with diameter of 15 µm (green), 10 µm (red) and 5 µm (cyan) at different locations. Reproduced from Ref. (Su et al., 2021) with the permission of the Royal Society of Chemistry. (C) The Simulation results for trajectories of three particles with diameters including 6, 10 and 15 µm corresponding to RBCs, WBCs and CTCs, respectively at flow rates of 1000 µl/min. Reproduced from Ref. (Nasiri et al., 2022) with the permission of Elsevier.

4.3. Design of coupling microfluidic device

Except for general particles, the numerical simulation also be used for the proposal of coupling microfluidic device for separation of magnetic microparticles in recent years (Alnaimat & Mathew, 2020; Ruan et al., 2022). Alnaimat et al. (Alnaimat & Mathew, 2020) proposed a novel a magnetophoretic microfluidic device for separation of magnetic microparticles based on size. The schematic diagram of particle motion in the vertical plane is shown in Figure 10A I. And the corresponding trajectories of the microparticles in the microchannel were obtained by simulation analysis, as shown in Figure 9A II. For multi-magnetic particles separation, Ruan et al. (Ruan et al., 2022) proposed an arc comb structure microfluidic (ACSM) chip and carried out simulation-based parametric analysis of the new ACSM chip. By numerical simulations, the parameters (e.g. flow velocity, velocity ratio of the two inlets, the remnant flux density of the magnets, and dynamic viscosity of the buffer solution) affecting the particle trajectory and separation efficiencies were investigated and optimized. The distribution of magnetic field and velocity field in the separation chamber obtained by simulation is shown in Figure 10BI. And the particles trajectories (Figure 10B II) shows that the larger particles prefer to exit from bottom outlet because they experience bigger magnetic force and hydrodynamic forces. According to investigating and optimizing the separation conditions, six kinds of magnetic particles were separated based on size by the ACSM chip. This research reveals that the relevant simulation results have a significant effect on the fundamental understanding of particle dynamic behaviour in the chip chamber.

Figure 10. (A) I. Schematic diagram of particle motion in a vertical plane of a microfluidic device. II. Particles separation in microfluidic device. Red line represents 10 µm microparticle and blue line represents 1 µm microparticle. Reproduced from Ref. (Alnaimat & Mathew, 2020) with the permission of Taylor & Francis Group. (B) The numerical simulations result. I. Distribution of magnetic field and velocity field in the separation chamber at Br = 1.2 T, υ1 = 2.9 mm/s and υ2 = 8.7 mm/s. Br is Remanent flux density, υ1 and υ2 are average flow velocity at the inlet 1 and the inlet 2. II. Simulation results of the particles trajectories in the specific separation chamber at υ1 = 0.8. mm/s, Br = 0.8 T. Reproduced from Ref. (Ruan et al., 2022) with the permission of Taylor & Francis Group.

The above numerical simulation researches used to evaluate the feasibility of microfluidic design were mainly aimed at revealing the particle focusing analysis of specific particles under specific conditions, so as to verify the effectiveness of its structure design and Re setting on separating particles and cells. Recently, to meet the design requirements of inertial microfluidic systems for efficient separation of various types of actual particles, inertial microfluidic have been integrated with active separation methods, and the diversity and complexity of microfluidic structures have continued to increase accordingly. A single numerical simulation study of inertial microfluidics is no longer practical, so it is urgent to establish a multi-field coupled microfluidic numerical simulation study to accurately predict the focus position of particles affected by multiple physical fields. Numerical simulation research as a supplement and alternative method for experimental verification, is a low-cost and efficient way to evaluate and optimize the reliability of microfluidic device design (Issakhov et al., 2019). However, owing to the number of existing numerical simulation methods and the commercial CFD software (Yan et al., 2021) are large and the applications are chaotic, researchers cannot choose quickly and accurately. Therefore, it is urgent to establish an applicability standard of numerical simulation method based on specific sorting particles to improve its friendliness to researchers.

5. Conclusion

In this review, we have briefly summarized all computational techniques for inertial microfluidic modeling in theoretical perspective. In these techniques, the DNS methods based on the Navier-Stokes equations have been intensively studied and developed for inertial particulate flows, which are highlighted in this review. These methods are mainly used for assessing the effects of different parameters on particle migration within a microchannel, including ALE, DLM, IBM, and FSPP. Generally, three aspects affect the motion of particles in channel: particle properties, the channel geometry and the fluid properties. Admittedly, there has been significant progress in inertial modeling of particle motion, especially for spherical particle, simple microchannel and Newtonian fluid. But various complex cases (e.g. non-straight channels, non-spherical/deformable particles and viscoelastic fluids) are still in their infancy and require thorough investigation. That’s precisely what the researchers have been focusing on in recent five years, as the microchannels used for sorting real particles such as cell and bacteria are getting more complicated. Compared to the migration of simple rigid spherical particles in straight Newtonian fluids, the modeling of complex case in line with the real situation is complicated and time-consuming. Furthermore, the simulation results obtained by the developed codes and models are usually inconsistent with the reality. Lately, by combining diverse DNS methods, combining DNS methods with LBM and incorporating asymptotic solutions into the LPT method, several new modeling strategies have been proposed to obtain more realistic simulation results of complex situations. Those coupling methods are indeed advanced for researchers to go deeply into the particles behaviour in realistic complex micorchannels. To date, the broad applications of numerical simulation-assisted microfluidic technology have indicated its promising prospects for microfluidics, single-cell analysis, and early-stage disease diagnostics, but there are still gaps between the numerical simulation and clinical application. First, the codes and models for complex situations should be proposed and improved accordingly. Second, the existing numerical simulation research methods are numerous and chaotic. To improve its friendliness to researchers, it is urgent to establish the applicability standard of numerical simulation method for specific particle and particular research. Last, as increasing investigators has used the commercial CFD packages, new commercial and scientific computational packages need to be developed to accurately predict the migration of realistic particles in non-Newtonian fluids.

Author contributions

Ruiju Shi conceived and written the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by Youth Innovation Fund of Ji Hua Laboratory: [Grant Number X220231XD220], the Science and Technology Project of Ji Hua Laboratory: [Grant Number X190191TD190].