2 Problem statement

In this work, we study flow of a dilute suspension of neutrally buoyant microspheres over a confined rectangular cavity using both microfluidic experiments and lattice-Boltzmann simulations. The cavity is connected to a rectangular channel. Due to the effect of inertial lift force inside the channel, particles equilibrate in close proximity to the channel walls (Chun & Ladd Reference Chun and Ladd2006; Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Tonner2009). The flow separation at the leading cavity wall leads to formation of a vortex bubble which can entrap particles from the channel (later we explain that ‘vortical flow’ is a precise terminology for the flow inside the cavity). The schematic of the flow geometry and coordinate system are shown in figure 1. Fluid flow and dynamics of suspended particles depend on cavity length, width and height; or cavity size along the $X$ , $Y$ and $Z$ directions respectively. The flow direction in the channel is specified by $X$ . The cavity length is abbreviated by $L$ , the width by $W$ and the height by $H$ . The cavity aspect ratio is defined as the ratio of cavity length to width, $\unicode[STIX]{x1D706}=L/W$ . Dynamics of particles is also affected by flow inertia. For fluid flow inside a rectangular channel, inertia is represented by $Re$ , defined using hydraulic diameter $h_{d}$ , average flow velocity $\langle u\rangle$ and fluid kinematic viscosity $\unicode[STIX]{x1D708}$ as $Re=\langle u\rangle h_{d}/\unicode[STIX]{x1D708}$ . The hydraulic diameter is defined using channel cross-sectional area $A$ and perimeter $P$ as $h_{d}=2A/P$ . Other length scales such as the primary direction of the velocity gradient ( $Y$ ) or particle diameter ( $d$ ) can also be used to compute $Re$ . In a pressure-driven flow, the shear rate varies inside the channel. Therefore, $Re$ at particle length scale can be defined using the average fluid velocity, the channel width $C_{Y}$ and particle diameter $d$ as $Re_{p}=(d/C_{Y})^{2}Re_{y}$ , where $Re_{y}=\langle u\rangle C_{Y}/\unicode[STIX]{x1D708}$ (Kahkeshani et al. Reference Kahkeshani, Haddadi and Di Carlo2016). For the rest of this study, we proceed with the common definition of $Re$ for rectangular conduits which is based on $h_{d}$ . We emphasize that $Re$ characterizes fluid inertia and is distinguished from particle inertia, which is often represented by the Stokes number, $St$ . Using particle density $\unicode[STIX]{x1D70C}_{p}$ , the fluid density $\unicode[STIX]{x1D70C}_{f}$ and $Re_{p}$ , the Stokes number can be defined as $St=(\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f})Re_{p}$ . Due to equal densities, neutrally buoyant particles have similar inertia to the fluid, i.e. $St=Re_{p}$ . Regardless of geometrical simplicity, it will be discussed that particle entrapment is affected by various factors. Finite size particles deviate from flow streamlines due to flow disturbance. Here, we study the effect of particle size on entrapment and dynamics inside the cavities for finite size particles. A size ratio is defined using the diameter of the particle $d$ and the channel width $C_{Y}$ , as $\unicode[STIX]{x1D6FF}=d/C_{Y}$ . The channel width, or the size of the channel in the $Y$ direction, is held constant at $C_{Y}=40~\unicode[STIX]{x03BC}\text{m}$ for all experiments.

Figure 1.

(a) Microfluidic geometry for studying the interaction of dilute suspensions with the flow cell inside the cavity. Directions of the coordinate system are indicated in the figure. The flow direction is specified by $X$ . (b) Two inlets are implemented in the channel entrance in order to infuse the buffer solution and suspension separately: for particle entrapment, buffer solution is infused first to generate the flow cell and the suspension is injected after stoppage of the buffer flow. For depletion studies, the infusion process is reversed.

The cavity size affects the number of entrapped particles. The geometry of an open cavity is similar to a backward facing step followed by a forward facing step, as it is illustrated in figure 1(b), with a primary vortex bubble extending between the leading and trailing cavity walls. We study entry of hard spherical particles into the cavity and follow their dynamics inside the flow cell. Our results add more details to explanations of particle motion inside recirculating flow in a tubular expansion presented by Karino & Goldsmith (Reference Karino and Goldsmith1977), abbreviated to K&G hereafter. Considering that K&G examined the motion of red blood cells, the particle motion was influenced by the deformation and orientation of red blood cells. Significant precision of microfluidic platforms enables us to study the dynamics of isolated particles inside the cavities, as opposed to the experiments of K&G which were affected by particle aggregates. Here, we study motion of neutrally buoyant hard spheres in order to isolate the physics of non-deformable particles and present a more clear interpretation. Shape diversity, size distribution and deformability in a blood mixture prevent a rigorous exploration of particle motion in experiments. Therefore, we only varied the diameter of particles to make progress in studying the particle entry and dynamics in the vortical flow. Studying the flow of deformable particles and capsules has turned into an active research area (Kumar, Rivera & Graham Reference Kumar, Rivera and Graham2014; Kumar & Graham Reference Kumar and Graham2015; Sinha & Graham Reference Sinha and Graham2016) and the motion of deformable human cells in recirculating flows can be the subject of a separate study. We discuss the effect of cavity aspect ratio by changing cavity length and width. We point out in § 7.2 that confinement of the cavity in the $Z$ direction induces spiralling fluid flow which can affect particle entrapment. Besides, varying cavity size in the $Z$ direction leads to change of the particle distribution inside the channel. Therefore, and to keep the distribution of particles inside the channel identical for all experiments, we only present the effect of $\unicode[STIX]{x1D706}$ for devices of the same size in the $Z$ direction, i.e. $H=70$ ( $\unicode[STIX]{x03BC}\text{m}$ ). More discussion on the effect of confinement in the $Z$ direction requires details of the fluid flow and particle distribution inside the channel.

We study flow of dilute suspensions of $d=10$ , $15$ and $20~\unicode[STIX]{x03BC}\text{m}$ particles, corresponding to $\unicode[STIX]{x1D6FF}=0.25,0.375$ and $0.5$ . The values of $\unicode[STIX]{x1D6FF}$ indicate that particle diameters are comparable to the channel width; therefore, disturbance generated by particles is no longer negligible and particles are considered ‘finite size’. The motion of a finite size particle is influenced by hydrodynamic forces of the disturbance flow and can be drastically different from tracers which primarily follow streamlines. For practical applications, such as blood flow in arterial stenosis or aneurysms, or fluidic devices for separation of tumour cells from blood samples, particle size ranges between the size of red blood cells with diameter $d\simeq 6{-}8~\unicode[STIX]{x03BC}\text{m}$ , white blood cells with $d\simeq 8{-}12~\unicode[STIX]{x03BC}\text{m}$ and cancer cells, which are normally larger than other blood components with $d\simeq 10{-}25~\unicode[STIX]{x03BC}\text{m}$ . Our range of particle diameter is close to the size of components in blood flow and microfluidic cell sorting and separation applications; although we present motion patterns which can be extended to any particle size. In this work, we vary flow rate in the range of $100{-}1000~\unicode[STIX]{x03BC}\text{L}~\text{min}^{-1}$ corresponding to $Re=30{-}308$ . We have scaled all velocities by $\langle u\rangle$ which is obtained as the ratio of volumetric flow rate $Q$ to the channel cross-sectional area $\langle u\rangle =Q/A$ . Further increase of the flow rate leads to large deformation of the polydimethylsiloxane (PDMS) microchannels (Sollier et al. Reference Sollier, Go, Che, Gossett, O’Byrne, Weaver, Kummer, Rettig, Goldman and Nickols2014) and leakage from inlets. Based on the definition of $Re_{p}$ and for $d=20~\unicode[STIX]{x03BC}\text{m}$ particles, $Re_{p}=5.6{-}56.9$ in the range of flow rates examined in this study. Considering the range of $Re_{p}$ , inertial effects should be fully considered in studying particle dynamics. To eliminate the effect of collision induced fluctuations on particle entry (Haddadi et al. Reference Haddadi, Shojaei-Zadeh and Morris2016), the volumetric fraction of solid phase in suspensions is $\unicode[STIX]{x1D719}=0.0001$ . The results remain invariant of $\unicode[STIX]{x1D719}$ by testing $\unicode[STIX]{x1D719}=0.001$ suspensions.

The manuscript is organized as follows: the experimental procedure is presented in § 3, followed by a brief explanation of the governing equations and computational parameters for discrete particle simulations with the lattice-Boltzmann method in § 4. We begin our discussion by presenting our flow visualization experiments in § 5. Dynamics of an isolated particle inside the cavity and the limit cycle trajectory will be presented in § 6. Then we study flow of a dilute suspension over the cavities in § 7 by discussing particle entry and exit from the cavities and the influence of the limit cycle trajectory on particle entrapment. We conclude the manuscript by our final remarks in § 8.

3 Experimental procedure

The microfluidic device in this study is composed of a single rectangular channel connected to a cavity, which has been fabricated by a standard PDMS replica moulding process (McDonald et al. Reference McDonald, Duffy, Anderson, Chiu, Wu, Schueller and Whitesides2000). A master mould was fabricated by spin coating SU-8 photoresist (MicroChem) on a silicon substrate followed by exposure through computer-aided design masks. Using the master moulds, the PDMS microchannels were fabricated by casting a Sylgard 184 Elastomer Kit. The inlet and outlets are punched into the PDMS using a pin vise. In order to enclose the channels, the PDMS microfluidic channel is bonded to the glass slide by surface activation of PDMS and glass using oxygen plasma. The final device is annealed in an oven at $65\,^{\circ }\text{C}$ for 24 h before use. The error of the fabrication procedure in the $X$ and $Y$ directions was determined by microscopy. In the $Z$ direction, the surface profile of the mould was measured using a Veeco Dektak surface profiler. In all directions, the error of the fabrication process is below $1~\unicode[STIX]{x03BC}\text{m}$ .

To address the effect of cavity width and length separately, we have fabricated cavities with constant width ( $W=310~\unicode[STIX]{x03BC}\text{m}$ ) and various lengths ( $L=620{-}1550~\unicode[STIX]{x03BC}\text{m}$ ); similarly to study the effect of cavity width we kept the length constant at $L=930~\unicode[STIX]{x03BC}\text{m}$ and changed the width of cavities ( $W=310{-}660~\unicode[STIX]{x03BC}\text{m}$ ). Table 1 summarizes device dimensions which are discussed in this work. Herein, we refer to the device aspect ratio $\unicode[STIX]{x1D706}$ indicated in this table. All of our results are presented for $H=70~\unicode[STIX]{x03BC}\text{m}$ cavities, although we have examined $H=45~\unicode[STIX]{x03BC}\text{m}$ devices to study the stable motions inside the cavities. We pointed out before that wall confinement in the $Z$ direction or height of the microfluidic device affects particle distribution inside the channel before entrapment into the cavities; therefore, by keeping the device height constant, we only examine the effect of $\unicode[STIX]{x1D706}$ on particle entrapment and dynamics.

Dilute suspensions of rigid particles are prepared by dispersing polystyrene spherical particles (Thermo Scientific) $(\unicode[STIX]{x1D70C}=1.02~\text{g}~\text{cm}^{-3})$ in a buffer solution composed of deionized water, glycerol ( $0.08v/v$ ) and a minute amount of triton X100 surfactant (Sigma-Aldrich). The density and kinematic viscosity of the buffer solution are $\unicode[STIX]{x1D70C}=1.023~\text{g}~\text{cm}^{-3}$ and $\unicode[STIX]{x1D708}=1.15\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ respectively. The initial fraction of the solid particles in the suspension supplied by the manufacturer is $0.01v/v$ . To prepare dilute suspensions, $20~\unicode[STIX]{x03BC}\text{L}$ of suspension is dispersed in 40 mL of buffer solution using vortex mixing. In order to conduct flow visualization experiments, $40~\unicode[STIX]{x03BC}\text{L}$ of a dilute suspension of $1~\unicode[STIX]{x03BC}\text{m}$ fluorescent polystyrene tracers, with initial solid fraction of $0.01v/v$ is dispersed in 40 mL of buffer solution. Solutions are infused into microchannels using a Harvard PHD 2000 syringe pump. We have utilized high-speed and fluorescent microscopy imaging to study interaction of particles with the vortex inside the cavity. The imaging details will be explained in the following sections.

Table 1.

Microfluidic devices used for investigation of particle interaction with the flow cell. We refer to the device aspect ratio $\unicode[STIX]{x1D706}=L/W$ in this work. Length, width and height are cavity dimensions along the $X$ , $Y$ and $Z$ directions respectively.

4 Governing equations and computational method

We present a brief description of governing equations and the lattice-Boltzmann method (LBM) for hard sphere suspensions. We refer to the background articles stated in the following for more details of the computational method. The LBM is developed for hard spherical particles dispersed in a Newtonian fluid. The governing equations of the fluid phase in dimensionless form are

where length can be made dimensionless by $h_{d}$ , velocity by $\langle u\rangle$ and pressure by $\unicode[STIX]{x1D702}\langle u\rangle /h_{d}$ , where $\unicode[STIX]{x1D702}$ is dynamic viscosity of the fluid. Newtonian dynamics governs the motion of particles. The translational ( $\boldsymbol{V}_{i}$ ) and rotational ( $\unicode[STIX]{x1D734}_{i}$ ) velocities of particles are calculated from where $\boldsymbol{F}_{i}$ and $\boldsymbol{T}_{i}$ are force and torque on particle $i$ , which has mass $m_{i}$ and moment of inertia $I_{i}$ .

The LBM model for hard sphere suspensions was developed by Ladd (Reference Ladd1994a ,Reference Ladd b ), with further modifications by Aidun, Lu & Ding (Reference Aidun, Lu and Ding1998) for dynamics of moving particles with arbitrary densities. The hard sphere LBM is based on a discrete Boltzmann equation on Cartesian grids coupled with Newtonian dynamics for the motion of particles. Theoretical details of the LBM model used in this work have been discussed in Ladd (Reference Ladd1994a ,Reference Ladd b ) and in a review article by Aidun & Clausen (Reference Aidun and Clausen2010). The method has been successfully utilized for computing rheological properties of suspensions (Haddadi & Morris Reference Haddadi and Morris2014) and studying the effect of inertia at bulk flow, such as flow of suspensions around square and circular cylinders (Haddadi et al. Reference Haddadi, Shojaei-Zadeh, Connington and Morris2014, Reference Haddadi, Shojaei-Zadeh and Morris2016).

In flow of a dilute suspension over and inside an open cavity, particles may collide with one another and with cavity walls. Therefore, particle motion at small separations from solid walls should be resolved using hydrodynamic lubrication and collision forces. Lubrication force at small gaps between solid surfaces has been implemented by construction of the grand resistance matrix for lubrication interaction between a spherical particle and a plane wall, based on the assumption that fluid flow at small gaps is inertialess (Nguyen & Ladd Reference Nguyen and Ladd2002). However, the hydrodynamic lubrication force does not prevent overlap between solid surfaces, especially for collision of the spherical particle with the walls and particle–corner collisions at high inertia. Various methods have been employed to resolve motion of particles at small surface separations (Ardekani & Rangel Reference Ardekani and Rangel2008). In Haddadi et al. (Reference Haddadi, Shojaei-Zadeh and Morris2016) a model based on a combination of lubrication force and elastic collision has been discussed for computing the particle motion at small gaps between solid surfaces, especially collision of a spherical particle with a sharp corner. According to this model the gap between solid surfaces is defined as $h$ and lubrication force governs the motion of particles for separations smaller than a lubrication cutoff, defined as $h_{N}=0.0015{-}0.05d$ , depending on the kinematic viscosity. For particle overlap defined as $h\leqslant 7.5\times 10^{-7}$ , an elastic collision force governs motion of particles.

Further details of the computational method are discussed in Haddadi et al. (Reference Haddadi, Shojaei-Zadeh and Morris2016). We have utilized the same tool in the present study and we refer to the mentioned article for further explanations. The confined cavity geometry is implemented by placement of two rectangular objects inside the computational box. The axial flow is along the $X$ direction and the box is bounded by walls in the $Y$ and $Z$ directions. The computational box is periodic along the $X$ direction. The majority of numerical simulations in our work are presented for $\unicode[STIX]{x1D706}=3$ cavities, with sizes $930\times 310\times 70~\unicode[STIX]{x03BC}\text{m}$ along the $X$ , $Y$ and $Z$ directions, respectively. The numerical simulations for the $\unicode[STIX]{x1D706}=3$ cavity are conducted in a computational box of $768\times 176\times 35$ lattice units (l.u.) along the $X$ , $Y$ and $Z$ directions. The results are invariant of resolution by repeating simulations in a larger $1536\times 352\times 70$ l.u. domain. By placing rectangular blocks in a $768\times 176\times 35$ l.u. domain, the leading cavity wall is located 150 l.u. from the inlet. Similarly, the trailing wall is symmetrically placed 150 l.u. away from the outlet. To simulate the motion of $20~\unicode[STIX]{x03BC}\text{m}$ particles, the particle resolution is chosen to be $10.08$ lattice units per diameter of particle (l.u. $/d$ ) in a $768\times 176\times 35$ l.u. domain.

Figure 2.

Fluid streamlines inside $\unicode[STIX]{x1D706}=3$ and $\unicode[STIX]{x1D706}=5$ cavities, with lengths $L=930$ and $1550~\unicode[STIX]{x03BC}\text{m}$ , at various $Re$ . The arrows specify direction of the flow ( $X$ ). Other dimensions of the cavities, i.e. width ( $Y$ ) and height ( $Z$ ), are identical ( $W=310$ and $H=70~\unicode[STIX]{x03BC}\text{m}$ ). Size details of the cavities are indicated in table 1. Forming a complete flow cell inside a longer cavity requires higher $Re$ .

5 Fluid flow over confined cavities

In this section we present images of fluorescent tracers in order to discuss fluid streamlines in the cavities. We demonstrate that vortex morphology depends on cavity size and $Re$ . The flow of tracer suspension is monitored using a Nikon Ti eclipse microscope and a Phantom v2010 camera. Depending on the flow rate, the exposure time changes from $1$ to 5 s, with higher exposure for higher flow rates. Due to their small size compared to the cavity width or length ( ${\sim}O(0.001)$ ), fluorescent tracers generate negligible flow disturbance and follow the fluid streamlines. For the cavities in the present study, distribution of tracers are visualized in great detail by taking magnified snapshots from multiple subregions and joining images to construct the whole picture of flow in the cavity. Flow of Newtonian fluids inside the cavity has always been a question of great scientific and technological interest. Despite its simple geometry, the flow inside a cavity can manifest various phenomena; such as corner eddies in Stokes flow (Moffatt Reference Moffatt1964) and vortical flows in turbulence (Citro et al. Reference Citro, Giannetti, Brandt and Luchini2015). Here we present the longitudinal extent of the flow cell inside the cavity. It is pointed out that the fluorescent imaging illustrates the flow morphology within a finite depth of field inside the cavity. The depth of field depends on the objective lens magnification and is within $\pm 5.8~\unicode[STIX]{x03BC}\text{m}$ from the object plane. We explain more details of the streamlines in § 7 and discuss the effect of wall confinement on formation of spiralling fluid flow. As was stated earlier, ‘vortical flow’, rather than recirculation, is a more precise terminology for the flow cell inside the wall-confined cavity.

We compare streamlines for different size cavities for a wide range of $Re$ . The influence of length on morphology of the flow cell is discussed first. Figure 2 displays streamlines inside $\unicode[STIX]{x1D706}=3$ and $5$ cavities. It is reminded that the width and height of cavities are identical, i.e. $W=310$ and $H=70~\unicode[STIX]{x03BC}\text{m}$ (table 1). The streamlines are compared for $Re=62$ , $123$ and $216$ for both cavities. At $Re=62$ , an incomplete flow cell forms in both $\unicode[STIX]{x1D706}=3$ and $5$ cavities. The incomplete flow cell denotes a triangular vortex bubble extending to a portion of the cavity length. With increasing $Re$ , the length of the vortex progressively increases until the vortex reaches the trailing wall. A vortex covering the full cavity length is termed a complete flow cell. Higher $Re$ is required to generate a complete flow cell for a longer cavity. By further increase of $Re$ to $216$ , the vortex core in the $\unicode[STIX]{x1D706}=3$ cavity approaches closer to the trailing cavity wall. The vortex core is a qualitative definition for the central zone of the vortex where tracers move in small orbits. In figure 3, we present streamlines for $\unicode[STIX]{x1D706}=3$ and $\unicode[STIX]{x1D706}=2.02$ cavities to compare the effect of cavity width on fluid flow at identical lengths and heights of $L=930$ and $H=70~\unicode[STIX]{x03BC}\text{m}$ . At each $Re$ , the location of the vortex core along the flow direction $X$ is relatively the same for cavities with different width.

Figure 3.

Fluid streamlines inside $\unicode[STIX]{x1D706}=3$ and $\unicode[STIX]{x1D706}=2.02$ devices, showing the fluid flow inside $W=310$ and $460~\unicode[STIX]{x03BC}\text{m}$ cavities. The length and height of both cavities is identical, i.e. $L=930~\unicode[STIX]{x03BC}\text{m}$ and $70~\unicode[STIX]{x03BC}\text{m}$ respectively. The streamlines are compared for $Re=62$ , $123$ and $216$ . The arrows indicate the flow direction inside the channel.

6 Dynamics of an isolated finite size particle inside the cavity

In this section, we present microfluidic measurements of the isolated particle dynamics inside the cavity. The LBM simulations are also utilized to study aspects of particle motion. We demonstrate that an isolated finite size particle forms a stable ‘limit cycle’ trajectory inside the cavity and discuss the topology of the limit cycle for different size cavities for various $Re$ . By referring to the equation of motion of a finite size particle, it is discussed that formation of the limit cycle is a distinct feature of finite size particles; considering that a finite size particle is affected by hydrodynamic contributions to the particle force that are negligible for tracers.

Figure 4.

Experimental trajectory of an isolated particle inside $\unicode[STIX]{x1D706}=3$ (a,b), $2.02$ (c,d) and $5$ (e,f) cavities at $Re=123$ and $216$ . The cavity size is made dimensionless using the diameter of the large, i.e. $d=20~\unicode[STIX]{x03BC}\text{m}$ , particle as the length scale. The colour bar specifies the particle velocity scaled by the average flow velocity inside the channel $\langle u\rangle$ .

Motion of an isolated particle inside the cavity is studied by analysis of high-speed camera recordings of a Nikon eclipse Ti7 microscope and a Phantom v711 high-speed camera. To compute the trajectory of an isolated particle in the cavity, a dilute suspension ( $\unicode[STIX]{x1D719}\sim O(10^{-4})$ ) is infused into the microchannel. Subsequently after entrapment of one particle, which can be monitored in real time imaging, infusion of the suspension is stopped. The vortex inside the cavity is maintained by flow of the buffer solution from the second inlet after stoppage of the suspension flow. Motion of the particle is monitored over a long time span of 10 min and is recorded every 1 min. Each video contains a stack of images recorded in a time span of 4 s, at rates of $4000$ , $5000$ or $6000$ frames per second (f.p.s.), depending on $Re$ (higher frame rate for higher $Re$ ). The particle trajectory is later computed by image analysis using a Hough transformation algorithm to track centre of mass location of the particle.

We present the topology of particle trajectories for $\unicode[STIX]{x1D706}=3$ , $\unicode[STIX]{x1D706}=2.02$ and $\unicode[STIX]{x1D706}=5$ cavities in figure 4. The cavity size is made dimensionless by the diameter of the largest particle, i.e. $d=20~\unicode[STIX]{x03BC}\text{m}$ . The colour bar specifies the particle velocity, which is scaled by the average flow velocity in the channel $\langle u\rangle$ . An isolated finite size particle orbits on a stable trajectory without deviation. The particle accelerates and decelerates as it traces its stable orbit. Due to the large magnitude of axial flow velocity inside the channel, the particle velocity increases as it moves forward in close proximity of the channel stream. The particle velocity is lower near the bottom wall of the cavity. Changing cavity size or $Re$ alters the fluid flow and consequently results in change of the limit cycle topology.

Stable motion of an isolated particle on a closed orbit implies formation of a limit cycle trajectory. The time scale of one complete orbit on the stable orbit $t^{\ast }$ depends on $\unicode[STIX]{x1D706}$ and $Re$ . For instance, in a $\unicode[STIX]{x1D706}=3$ cavity, the time scale of one complete orbit at $Re=154$ is $t^{\ast }=0.009$  s. Therefore, in a time span of 10 min, the particle completes at least $66\,000$ orbits without deviation. However, it can also be argued that the stable motion is observed within the limited time span of the experiment and in sufficiently long times, the particle may continue its outward motion and exit the cavity; or move towards other orbits inside the cavity. We have utilized lattice-Boltzmann simulations to corroborate the presence of a unique limit cycle orbit inside the cavities. The limit cycle is a trajectory attractor with a basin of attraction extending in both interior and exterior regions of the limit cycle. Trajectories in the interior zone spiral outward and exterior trajectories spiral inward to approach the limit cycle. Therefore, the limit cycle defines a bound for outward and inward spiralling motions. Further spiralling motion of trajectories beyond the limit cycle orbit leads to overlap of inward and outward spirals, which is not permitted by the mass conservation principle.

Figure 5.

Simulations of spiralling motion of (a) interior and (b) exterior trajectories of a $d=10.08$ (l.u.) particle towards the limit cycle orbit at $Re=123$ inside a $\unicode[STIX]{x1D706}=3$ cavity. (c) The basin of attraction extends to the channel, considering that a particle that enters the cavity spirals toward the limit cycle trajectory. (d) The LBM simulation of the limit cycle trajectory at $Re=123$ . In the background, the fluid streamlines at $z=1.73$ plane are exhibited. The fluorescent image of flow for $Re=123$ is displayed in figure 2(c). The symbols in (a,b) show the location of the initial point.

Numerical simulations are conducted by releasing an isolated particle from various initial locations inside the cavity after reaching a steady flow. We present the LBM simulations of particle trajectories at $Re=123$ . The fluid streamlines and the limit cycle topology of microchannel experiments at $Re=123$ are shown in figures 2(c) and 4(a) respectively. As shown in figure 5(a), a particle released from the interior region (closer to the core of the vortex) spirals outward to reach the limit cycle motion. Exterior trajectories spiral inward towards the limit cycle, as exhibited in figure 5(b). Considering that particles can also enter the cavity from the channel and approach the stable orbit, the basin of attraction of a limit cycle extends outside the cavity and into the channel. As seen in figure 5(c), a particle moving inside the channel enters the cavity and reaches a stable motion on the limit cycle trajectory. We later discuss conditions of particle entry into the cavities.

We state the difference between streamlines and limit cycle orbits obtained by LBM simulations and microfluidic experiments. Figure 6 exhibits LBM simulation of fluid flow computed on the $z=1.73$ symmetry plane along with fluorescent imaging of tracer distributions. It is stated that fluorescent imaging of tracers provides a picture of streamlines in a finite depth of field. Although the vortex morphology is relatively the same, the location of the vortex core and topology of the limit cycle is slightly different in the LBM simulations and microfluidic experiments. The difference can be due to the deformability of PDMS leading to slight change in the size of the microchannels, microgravity effects due to small density differences between polystyrene microspheres and the buffer fluid and the finite resolution of LBM simulations.

Figure 6.

Comparison between LBM simulation and microfluidic experiments of fluid flow and single particle dynamics inside the $\unicode[STIX]{x1D706}=3$ cavity at $Re=123$ . In (a) we demonstrate LBM simulations of fluid streamlines, which can be compared to distribution of tracers in (b). Difference between topology of the limit cycle orbit in LBM simulations and microfluidic experiments is shown in (c). Slight difference is observed between location of the vortex core and topology of the limit cycle in LBM simulations and microfluidic experiments.

In order to explain the stable motion of a finite size particle on the limit cycle trajectory, we discuss the hydrodynamic force generated on a finite size particle inside the cavity. Small particles follow streamlines closely, however, the hydrodynamic stress of the disturbance flow of a finite size particle leads to deviation from streamlines. We first use the equation of motion of the particle in order to explain various contributions on the particle force, leading to deviation from fluid streamlines. We later address the fact that inertial contributions are implicitly stated by this equation. For a finite size particle suspended in a non-uniform flow, the equation of motion has been presented by Maxey & Riley (Reference Maxey and Riley1983) as a form of Newton’s second law. The motion of a spherical particle with radius $a$ is described by

where $\boldsymbol{V}$ is the centre of mass velocity of the particle located at $\boldsymbol{Y}(t)$ at time $t$ and $\boldsymbol{u}$ is the undisturbed fluid velocity. The mass of the fluid and the particle are specified by $m_{f}$ and $m_{p}$ respectively. The dynamic viscosity of the suspending fluid is denoted by $\unicode[STIX]{x1D702}$ . The undisturbed fluid flow generates a force equivalent to $m_{f}(\text{D}\boldsymbol{u}/\text{D}t)|_{\boldsymbol{Y}(t)}$ on the particle at the particle centre of mass $\boldsymbol{Y}(t)$ ; the net force as would act on a deformable fluid element of the same size. The force of the undisturbed flow is specified by the first term on the right-hand side of equation (6.1). An accelerating immersed object must also accelerate the surrounding fluid, generating a force on the object denoted as the added mass which is represented by the second term on the right-hand side of the equation of motion. The third term describes the body forces on the particle; which, for instance, can be a gravitational force $\boldsymbol{f}=(m_{p}-m_{f})\boldsymbol{g}$ . The hydrodynamic stress on the surface of a non-deformable particle generates the hydrodynamic drag which appears in all examples of particle motion in a fluid. The last term on the right-hand side corresponds to this hydrodynamic drag force. It should be pointed out that we have not included the Basset history force in (6.1), considering that the history force affects the object for a motion starting from rest. The terms containing the Laplacian of the fluid velocity field $\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}$ are the Faxén corrections to take the curvature of the velocity field into account. The curvature of the fluid velocity cannot be neglected across the length of a large particle. Except for the forces of the undisturbed fluid velocity, i.e. $m_{f}(\text{D}\boldsymbol{u}/\text{D}t)|_{\boldsymbol{Y}(t)}$ , other terms of the equation of motion are obtained by assuming that the flow is inertialess, i.e. Stokes flow. Therefore, the particle acceleration appears only implicitly in time derivatives of the particle velocity $\text{d}\boldsymbol{V}/\text{d}t$ . Considering the range of $Re_{p}$ for flow of finite size particles in microchannels ( $Re_{p}=5.9{-}59.5$ for $d=20~\unicode[STIX]{x03BC}\text{m}$ particles in a $C_{Y}=40~\unicode[STIX]{x03BC}\text{m}$ channel), the inertial contributions to the particle dynamics are significant compared to viscous drag forces. Besides, the effect of wall and inertial lift force are not considered by (6.1). We point out that both inertial lift and wall forces can substantially influence the motion of particles inside the cavity.

It is pointed out that the acceleration of a hard spherical particle is represented by the centre of mass acceleration. The acceleration of each material point inside the rigid body with distance $\boldsymbol{r}$ from the particle centre of mass $\boldsymbol{a}=\unicode[STIX]{x1D736}+\unicode[STIX]{x1D734}\times \boldsymbol{r}+\unicode[STIX]{x1D74E}\times \unicode[STIX]{x1D74E}\times \boldsymbol{r}$ , where $\unicode[STIX]{x1D736}$ denotes the centre of mass acceleration, $\unicode[STIX]{x1D734}\times \boldsymbol{r}$ the azimuthal acceleration due to the angular acceleration $\unicode[STIX]{x1D734}$ and $\unicode[STIX]{x1D74E}\times \unicode[STIX]{x1D74E}\times \boldsymbol{r}$ the centrifugal acceleration due to the angular velocity $\unicode[STIX]{x1D74E}$ , respectively. For a spherical particle with a symmetric distribution of material points around the centre of mass, both azimuthal and centrifugal accelerations are zero and the centre of mass acceleration $\unicode[STIX]{x1D736}$ represents acceleration of the object. By referring to (6.1), we intend to emphasize that contributions to the hydrodynamic force generated by finite size and non-deformability result in significant deviation from fluid streamlines. Therefore, as opposed to tracers which follow fluid streamlines and spiral inside the cavity, a finite size particle approaches a stable orbit.

The presence of a limit cycle trajectory for a finite size particle in recirculating flow has been discussed in a number of works to result from the balance of different forces on the particle. For a particle orbiting inside the recirculating flow behind a bluff body, increasing particle size and density leads to outward displacement of the particle in each orbit (Haddadi et al. Reference Haddadi, Shojaei-Zadeh, Connington and Morris2014). The outward motion of the finite size particle is due to particle inertia, which increases with particle size and density (Haddadi et al. Reference Haddadi, Shojaei-Zadeh and Morris2016). Consequently, the particle does not follow closed streamlines of the recirculating wake and migrates outward. The limit cycle orbit can only be explained if there are opposing migration mechanisms. In K&G, the large residence time of red cell clusters was argued to result from the inertial lift force inside the tube, which generates an inward migration towards the vortex core, opposing outward motion of the cell clusters in the vortical flow region. The interplay between migration due to centrifugal force and outward motion of high density particles in the vortices of a Taylor–Couette device leads to formation of stable limit cycle trajectories and particle segregation (Wereley & Lueptow Reference Wereley and Lueptow1999). It has been discussed in Haddadi et al. (Reference Haddadi, Shojaei-Zadeh, Connington and Morris2014) that the outward inertial motion of a finite size particle inside the wake behind a bluff body is opposed by hydrodynamic lubrication force between a finite size particle and the obstacle. The hydrodynamic lubrication between surfaces of the obstacle and the particle prevents further outward displacement. Therefore, a stable limit cycle orbit forms close to wake boundaries.

Figure 7.

The inward migration generated by inertial lift force opposes outward migration of the particle towards the boundaries of the flow cell. The combined effects of outward and inward migrations lead to formation of a limit cycle trajectory inside the cavity. The inertial lift force on the particle is generated by velocity gradients, which are exhibited for three sections, $x=10$ , $20$ and $30$ inside the cavity at $Re=123$ . The direction of the lateral lift force is indicated in the figure.

Here we present a qualitative explanation for the motion of isolated particles on the stable limit cycle trajectory that is based on the effect of inertial lift force balancing the outward inertial motion of particles. However, due to the absence of a mechanistic description for lift force, our proposed mechanism is a baseline description. In figure 7 we present numerical simulations of the limit cycle trajectory at $Re=123$ . It has been shown that inertial particles migrate outward in the wake behind bluff bodies (Haddadi et al. Reference Haddadi, Shojaei-Zadeh and Morris2016). The rate of outward motion increases with particle inertia. Similarly, inertial particles released in the core of the cavity vortex migrate outward towards the vortex boundaries (figure 5 a). In each cycle of motion on the limit cycle orbit, the particle is subjected to inertial lift force due to the spatially varying velocity, which repels the particle from zones of steep velocity gradient. Profiles of axial velocity at different sections in the $X$ direction are shown in figure 7. The velocity gradient inside the channel exerts a lateral force while the particle traces the orbit in the forward direction, leading to inward migration towards the vortex core. Upon reversing direction, the particle is subjected to a velocity gradient in the opposite direction. The velocity gradient in the opposite direction is weaker, however it can also generate inward migration away from the vortex boundaries towards the vortex core. Therefore, outward displacement of a particle in the flow cell is opposed by inward migration generated by inertial lift force and particle reaches a stable orbit.

The limit cycle topology varies with particle size. We use scaling arguments of the lift in a straight channel to qualitatively explain variations of the limit cycle topology. In all scalings, there is a power-law dependence between $F_{L}$ and $d$ . Numerous scaling arguments delineate the relationship between $F_{L}$ and $d$ in tube or slit geometries (Saffman Reference Saffman1965; Cox & Brenner Reference Cox and Brenner1968; Ho & Leal Reference Ho and Leal1974; Asmolov Reference Asmolov1999). However, the migration pattern is different in a microchannel. As opposed to inertial pipe flow, which generates migration towards an annulus (Matas et al. Reference Matas, Glezer, Guazzelli and Morris2004), particles migrate towards discrete inertial equilibrium positions in a rectangular channel (Chun & Ladd Reference Chun and Ladd2006; Di Carlo et al. Reference Di Carlo, Edd, Humphry, Stone and Tonner2009; Miura, Itano & Sugihara-Seki Reference Miura, Itano and Sugihara-Seki2014). By solving the fluid flow around a force free particle coupled with particle dynamics, Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Tonner2009) presents $F_{L}\sim d^{3}$ for an experimentally relevant range of $1\leqslant Re_{p}\leqslant 20$ and $0.05\leqslant \unicode[STIX]{x1D6FF}\leqslant 0.2$ . Using asymptotic expansions in a three-dimensional channel geometry, Hood, Lee & Roper (Reference Hood, Lee and Roper2015) have presented $F_{L}\sim d^{4}$ scalings for $\unicode[STIX]{x1D6FF}\rightarrow 0$ . They show that with increasing $\unicode[STIX]{x1D6FF}$ , $F_{L}\sim d^{5}$ provides more accurate relationship between magnitude of lateral force and particle diameter. In figure 8, we present the limit cycle trajectories for $d=10$ and $20~\unicode[STIX]{x03BC}\text{m}$ particles, corresponding to $\unicode[STIX]{x1D6FF}=0.25$ and $0.5$ , inside $\unicode[STIX]{x1D706}=3$ cavity. At $Re=123$ , the limit cycle orbit of the larger $d=20~\unicode[STIX]{x03BC}\text{m}$ particle forms closer to the vortex core and $d=10~\unicode[STIX]{x03BC}\text{m}$ particle orbits closer to vortex boundaries. Although for intermediate values of $Re=185$ and $216$ , the difference between orbit topologies of $d=10$ and $20~\unicode[STIX]{x03BC}\text{m}$ particles is not pronounced, there is a clear distinction between large and small particle limit cycle orbits at $Re=308$ .

Figure 8.

Experimental measurements of the limit cycle trajectory for $d=10$ (○) (red) and $d=20$ (○) (blue) particles inside $\unicode[STIX]{x1D706}=3$ cavity at (a)  $Re=123$ , (b)  $Re=185$ , (c)  $Re=216$ and (d)  $Re=308$ . Diameter of the large particle $d=20~\unicode[STIX]{x03BC}\text{m}$ is used to make the lengths dimensionless.

We mention that the limit cycle orbit is a unique trajectory attractor which forms on the symmetry plane along the $Z$ direction, i.e. $z=H/2$ plane. All trajectories spiral towards the middle plane due to the flow symmetry. In figure 9, we exhibit LBM simulations of spiralling motion of a $\unicode[STIX]{x1D6FF}=0.5$ particle towards the midplane along the $Z$ direction ( $z=1.73$ ). The particle spirals towards the limit cycle on the midplane and maintains its planar motion on a stable limit cycle orbit.

Our results for motion of isolated particles inside the vortex show distinctive behaviour of finite size particles in the vortical flow inside a wall-confined cavity. The limit cycle forms for both $H=45$ and $70~\unicode[STIX]{x03BC}\text{m}$ cavities, although we presented our results for a $H=70~\unicode[STIX]{x03BC}\text{m}$ cavity. The stable motion of the finite size particle significantly increases the particle residence time inside the cavity and generates a particle retention zone. In addition to biological applications such as blood flow in aneurysms or stenosis geometries, formation of the limit cycle in a flow cell has important consequences in some industrial processes. For example, in cleansing indented surfaces from particulate pollutants, the presence of retention zones can lead to imperfect release of particles. In addition to particle size, the magnitude of the lift force depends on $Re$ . Therefore, there is a critical $Re$ , denoted by $Re_{cr}$ , below which the limit cycle trajectory no longer forms. Large particles are more influenced by the lift force, resulting in formation of the limit cycle at a smaller $Re_{cr}$ . For instance in a $\unicode[STIX]{x1D706}=3$ cavity, the value of $Re_{cr}$ for $d=20~\unicode[STIX]{x03BC}\text{m}$ is $84$ . For smaller $d=15$ and $10~\unicode[STIX]{x03BC}\text{m}$ particles, the values of $Re_{cr}$ are $200$ and $216$ respectively.

Figure 9.

The limit cycle trajectory forms on the symmetry plane, i.e. the midplane along the $Z$ direction. Particles on other planes spiral towards the symmetry plane. The LBM simulation of spiralling motion of a $\unicode[STIX]{x1D6FF}=0.5$ particle released from $z=-d$ towards the limit cycle on the symmetry plane, projected on (a) $XY$ and (b) $XZ$ planes at $Re=123$ .

7 Dilute suspension flow over open cavities

In the following, we study flow of a dilute suspension over the cavity and discuss particle entrapment. The particle exchange between the cavity and channel is explained by investigation of particle entry and exit. Formation of the limit cycle, which has been discussed in previous sections, is crucial for maintaining particles inside the cavity after entrapment. We begin by discussing the effect of the cavity aspect ratio $\unicode[STIX]{x1D706}$ on the number of entrapped particles at various $Re$ . Then we explain different factors affecting the particle entrapment in cavities. We demonstrate that the number of entrapped particles inside the cavity reaches a maximum at a specific $Re$ . It will be discussed that the combination of vortex morphology and particle dynamics leads to maximum entrapment. Finally, depletion of the cavities generated by hydrodynamic interaction between particles will be presented.

8 Concluding remarks

In this work we have shown that finite size particles form stable limit cycle trajectories inside confined cavities. We discussed the effect of cavity aspect ratio ( $\unicode[STIX]{x1D706}=L/W$ ) and $Re$ on topology of the limit cycle orbit. We qualitatively explained that a finite size particle migrates outward in the flow cell. The inertial lift force $F_{L}$ which scales with particle size ( $F_{L}\sim d^{3}$ Di Carlo et al. (Reference Di Carlo, Edd, Humphry, Stone and Tonner2009), $F_{L}\sim d^{4}$ Hood et al. (Reference Hood, Lee and Roper2015)), generates inward migration opposing the outward migration of particles towards vortex boundaries. The limit cycle trajectory of a small particle forms closer to the vortex outer boundary and with increasing diameter, the particle orbits closer to the vortex core. In the absence of interaction between multiple particles, the particle motion on the limit cycle is stable. Our detailed discussions of the limit cycle trajectory, clarify the reason for significantly longer residence time of large particles inside the vortex, observed for the first time by Karino & Goldsmith (Reference Karino and Goldsmith1977). Formation of the stable limit cycle trajectory is important in biological and industrial applications, such as retention zones inside aneurysms and cleansing of indented surfaces from particulate pollutants.

We discussed that due to large $Re_{p}$ , inertial contributions on particle acceleration should be fully considered. We described the particle entry as the collective effect of flow morphology, particle acceleration and fluid exchange between the cavity and the channel, induced by wall confinement. However, distinct effect of each contribution is not clearly known. Although inertial lift forces generate lateral migration, leading to closer proximity of the particle and vortex, their effect is combined with other hydrodynamic contributions on particle acceleration such as added mass, proximity to wall and varying drag. At low $Re$ , small number of particles are entrapped due to combined effects of insufficient lateral migration of particles and small vortex size. With increasing $Re$ , entrapment of particles inside the flow cell is due to inertial lift force, deviation of particles from flow streamlines due to other hydrodynamic contributions and exchange of fluid between the cavity and channel. An inertial particle with $\unicode[STIX]{x1D70F}_{p}>0$ , does not instantaneously respond to steep changes of flow direction and continues its motion towards the flow cell. At high $Re$ , particles tend to follow the axial flow inside the channel, without manifesting sufficient deviation towards the vortex.

We have examined the effect of cavity aspect ratio and $Re$ on the number of entrapped particles. It has been observed that for each $\unicode[STIX]{x1D706}$ , the number of entrapped particles is maximized at a specific $Re$ , denoted by $Re_{max}$ . With increasing cavity length the number of particles inside the cavity increases and $Re_{max}$ shifts to higher values, i.e. a horizontal shift. With increasing the width, the number of entrapped particles increases, without a significant change of $Re_{max}$ (vertical shift).

Collisions between multiple entrapped particles generates deviation from stable limit cycle orbit and particles exit the cavity. We have examined particle exit and discussed that the depletion depends on the number of particles inside cavities. Finally, we have shown that the particle entrapment inside the cavities depends on the size of particles and by adjusting $Re$ , it is possible to entrap particles based on their size. This result can be used for separation of large human cells from a background of blood cells.