Utilization of direct forcing immersed boundary methods for the optimization of inertial focusing microfluidics Patrick Giolando1, Hui Ma2, Barrett Davis2, Tamara Kinzer-Ursem2, Steve Wereley1 1 School of Mechanical Engineering, Purdue University 2 Weldon School of Biomedical Engineering, Purdue University 1. Introduction Inertial focusing microfluidics have gained significant momentum in the last decade for their ability to separate and filter mixtures of particles and cells based on size [1-3]. However, the most important feature is that the separation is passive, without the need for external forces. At the heart of inertial focusing is the balance between counteracting lift forces: shear and wall- induced lift. Shear-induced lift is a product of the curvature of the fluid flow and the rotation of the particle in the flow, while wall-induced lift is generated by the disturbance of the fluid by the particle near a wall. This phenomenon was first observed by Segre and Silberberg for the focusing of particles in a pipe, and was later extended to the focusing of cells and particle in rectangular channels [4]. Taking advantage of inertial focusing we explore particle capture utilizing an expanded channel microfluidics chip design. By expanding a small region of the straight channel microvortices form in the well, which allows for size selective trapping of particles [1, 2]. Modeling the two-way coupled with traditional finite volume method (FVM) or finite element method (FEM) can prove costly as with each time step the body-fitted grid would have to be re- meshed [5, 6]. Immersed boundary methods (IBM) offer a cost effective solution, rather than solving the fluid equations on a body fitted grid they are solved on a regular cartesian grid and the boundary conditions are imposed on the fluid domain by the addition of a forcing term. The Navier-Stokes equations are solved for on this Eulerian grid, while the Newton-Euler equations that govern the motion of the particle are solved for on the Lagrangian grid that defines the surface of the particle. The direct forcing IBM accounts for the force and torque acting on the particle by the requirement of the predict fluid velocity to be the local velocity of the particle on the surface of particle [7, 8]. Interpolations are required to find the velocity of the fluid on the Lagrangian nodes, and then a spreading function to return the force on the Lagrangian nodes back to the Eulerian grid [9]. 2. Methods 2.1 Experimental Methods 2.1.1 Design and Fabrication of Microfluidics Chips For this study, two chip designs were produced to evaluate the inertial focusing properties of particles. The first chip design is a straight channel microfluidics chip fabricated in polydimethylsiloxane (PDMS, SYLGARD184 Silicone Elastomer Kit) using a master mold developed with standard photolithography. A silicon wafer with dry resist film (PerMX 3050 series, DuPont Electronic Technologies) was used to develop the mold and cast the PDMS, which was then bound to a glass slide (Sigma Aldrich) and baked for 45 minutes at 95 ֯C. A second series of chip designs were produced with an expanded channel to produce microvortices for particle retention. These microfluidics chips were designed by cutting pressure sensitive adhesive (PSA, ARseal 90880, Adhesive Research, cyclic olefin copolymer (COC, Zeon Zeonor zf14-188) with a laser cutter (Universal Laser Systems, VSL350). 2.1.2 Particle Imaging Experimental data for the motion of buoyant polystyrene microparticles of size 1 and 7.32 µm in diameter, and acrylic particles of size 20 µm in diameter (Fisher Scientific) were captured using a fluorescence microscope (Axio Observer, Zeiss). Shutter speeds were varied from 1/10 to 1/10000 s, allowing for both particle image velocimetry (PIV) data analysis and the observation of complete pathlines. The 1, 7.32, and 20 µm particles were mixed with deionized water at a concentration of 5x107 , 2.5x105 , and 6x104 particles per mL, respectively. Particles were collected in a syringe (BD syringe, 1mL) and pumped through PEEK tubing (IDEX, 1569) utilizing a syringe pump (KD Scientific Inc) to vary flow conditions (0.1 ≤ Re ≤ 500). 2.1.3 Image Analysis The expanded well design produced two separate fluid domains with dramatically different Reynolds numbers (Re), which required separate in-house algorithms for the quantification of particle velocities. For regions of slow fluid motion, particle motion was captured in image pairs that produced a strong cross-correlation, which is ideal for PIV analysis. However, in regions of high Reynolds flow the particles became faint streaks, and cross-correlation procedures became nonviable. Particle streak velocimetry (PSV) was employed to find and quantify the length of the streaks to determine the particles velocity. The first few image pairs were used to segment the images into regions of large and small Re flow to be quantified by either PIV or PSV, respectively. The first image pair from the stack was split into interrogation window and evaluated with a FFT cross-correlation algorithm to evaluate the strength of cross-correlation for each window. Windows that produced cross-correlation resembling a Dirac function were quantified with PIV for the entire image stack. Regions with a cross-correlation domains of lowered and wider peaks in the cross-correlation domain were segmented with a morphological operator to look for streaks. If streaks were identified of significant size PSV was used to quantify the particle motion in the window for the entire stack. If there were no significant streaks in the window, then the particle velocity of that window was set to zero to avoid extraneous velocity vectors. This allowed for rapid and accurate quantification of the particle motion in the entire fluid domain within the microfluidics chip. Figure 1. PIV/PSV Algorithm. (A) Image pair used to produce mask to separate windows for PIV or PSV analysis. (B) Cross-correlation of a window in the image pair. (C) Segmentation of images to identify streaks. (D) Final quantification of particle velocities in the entire system. Regions quantified with PIV/PSV are red/blue, respectively. 2.2 Numerical Methods 2.2.1 Governing Equations Particle motion was governed by Newton-Euler equations, while the fluid motion was governed by the Navier-Stokes equations. The incompressible, Newtonian Navier-Stokes equations: 𝜕𝑢𝑖 𝜕𝑡 + 𝜕 𝜕𝑥𝑗 𝑢𝑖𝑢𝑗 = − 𝜕𝑝𝑖 𝜕𝑥𝑖 + 1 𝑅𝑒 𝜕 𝜕𝑥𝑗 𝜕𝑢𝑖 𝜕𝑥𝑗 (1) 𝜕𝑢𝑖 𝜕𝑥𝑖 = 0 (2) where 𝑢𝑖 is the ith component of the fluid velocity, 𝑥𝑗 is the jth dimension, and 𝑝𝑖 is the ith component of pressure. Notable this system was 2D flow and gravity was negligible in this system for small buoyant particles. Re is defined as Re = 𝜌fUL 𝜇 (3) where 𝜌f is the density, 𝜇 is the dynamic viscosity, U is the characteristic velocity, and L is the characteristic length of the fluid. The particle was modeled by a string of interconnected nodes, the motion of which includes both translational and rotational velocity: 𝑢𝑝𝑛 = 𝑢𝑝 +𝜔𝑝 𝑥 𝑟 (4) where 𝑢𝑝𝑛 is the velocity of the particle node, 𝑢𝑝 is the velocity of the particle, 𝜔𝑝 is the angular velocity of the particle, and r is the radial arm from the position of the node to the center of mass of the particle. Assuming buoyancy the Newton-Euler equations became: 𝜌𝑝𝑉𝑝 𝜕𝑢𝑝 𝜕𝑡 = ∮ 𝜏 ∗ 𝑛𝑑𝑎 𝜕𝑉 + 𝐹𝑐 (5) 𝐼𝑝 𝜕𝜔𝑐 𝜕𝑡 = ∮ 𝑟 𝑥 (𝜏 ∗ 𝑛)𝑑𝑎 𝜕𝑉 + 𝑇𝑐 (6) where 𝜌𝑝 is the density of the particle 𝑉𝑝 is the volume of the particle, 𝐹𝑐 is the force introduced by a collision, 𝐼𝑝 is the moment of inertia of the particle, 𝑇𝑐 is the torque introduced by a collision, and 𝜏 is the total stress tensor acting on the particle. Finally, eq 1 is modified to include the forcing term: 𝜕𝑢𝑖 𝜕𝑡 + 𝜕 𝜕𝑥𝑗 𝑢𝑖𝑢𝑗 = − 𝜕𝑝𝑖 𝜕𝑥𝑖 + 1 𝑅𝑒 𝜕 𝜕𝑥𝑗 𝜕𝑢𝑖 𝜕𝑥𝑗 + 𝑓𝑖 (7) where the forcing term f is zero everywhere except in the vicinity of the particle. 2.2.2 Numerical Methods The Navier-Stokes equations are handled with pressure-correction scheme, which is highly compatible with the direct forcing IBM [5, 10]. The numerical method is semi-implicit second order difference method, where operator splitting was used to handle the linear terms with the 2nd order backward difference formula (BDF-2) and the nonlinear terms with the 2nd order Adams-Bashforth (AB-2) explicit method. Spatially the viscous terms were approximated with the 2nd order central difference (CD-2) method, an adaptive upwind-downwind scheme is used to approximate the convective terms, and the Laplacian in the Poisson problem is approximated with a 9-point scheme (∝= 1/3) [11]. The velocity and pressure components are solved for on a fully staggered grid, also known as a Marker And Cell (MAC) scheme, Fig 1 bc. The first step of Chorin’s projection method is to solve for the intermediate fluid velocity without the added forcing term: 𝑢𝑖 ∗−𝑢𝑖 𝑛 ∆𝑡 = 1 𝑅𝑒 𝜕2𝑢𝑖 𝑛 𝜕𝑥𝑗𝜕𝑥𝑗 − 𝜕(𝑢𝑖𝑢𝑗) 𝑛 𝜕𝑥𝑗 (8) where 𝑢𝑖 ∗ is the intermediate fluid velocity, and is notably not divergence free. As mentioned previously this is solved in two steps: linear and nonlinear terms. The second step is to solve for the forcing term, which done by first interpolating the intermediate fluid velocity onto the Lagrangian nodes: 𝑈𝑙 ∗ = ∑ 𝑢𝑖𝑗 ∗ 𝛿𝑑(𝑥𝑖𝑗 − 𝑋𝑙 𝑛)∆𝑥∆𝑦𝑖𝑗 (9) where capital letters represent values on the Lagrangian grid, and 𝛿𝑑 is the Dirac delta function, Fig 2a,b. The forcing term on the Lagrangian nodes is then computed: 𝐹𝑙 𝑛+1/2 = 𝑈𝑝(𝑋𝑙 𝑛)−𝑈𝑙 ∗ ∆𝑡 (10) where 𝑈𝑝(𝑋𝑙 𝑛) is the velocity of the particle node, which includes both translational and rotational velocity. Equ 10 handles both the no slip and no penetration boundary conditions. Finally, the forcing term is interpolated back onto the cartesian grid: 𝑓𝑖𝑗 𝑛+1/2 = ∑ 𝐹𝑙 𝑛+1/2 𝛿𝑑(𝑥𝑖𝑗 − 𝑋𝑙 𝑛)∆𝑉𝑙𝑙 (11) where ∆𝑉𝑙 is the volume of the Lagrangian grid cells, Fig 2c. The third step is to update the intermediate velocity to account for the forcing term: 𝑢𝑖 ∗∗−𝑢𝑖 ∗ ∆𝑡 = 𝑓𝑛+1/2 (12) Figure 2. Interpolation Between Eulerian and Lagrangian Grid. The immersed boundary method utilized two independent grids to resolve the fluid-surface interface. A method for interpolating between the fully staggered Eulerian grid and the Lagrangian grid is achieved utilizing IBM. (A) The 4-point Dirac delta function the kernel is used to interpolate between the grids. (B) The interpolation of the intermediate fluid velocity onto the Lagrangian nodes. Dark blue triangles are the intermediate U-velocity nodes, which are known and used to find the intermediate fluid U- velocity at the Lagrangian node (large black dot), Equ 9. The Lagrangian forcing term is then found to impose the no- slip boundary conditions, Equ 10. (C) The spreading of the forcing term from the Lagrangian nodes onto the Eulerian nodes. The forcing term from each of the neighboring Lagrangian nodes (large black nodes) are interpolated onto the Eulerian node (dark blue triangle), Equ 11. Each node has a discrete volume, ΔV, associated with it such that the collection of nodes form a thin shell around the particle. In the fourth step, the projection function is solved: − 𝜕2𝑝𝑛+1 𝜕𝑥𝑗𝜕𝑥𝑗 = ( 𝜕𝑢𝑗 ∗ 𝜕𝑥𝑗 ) 1 ∆𝑡 (13) In the final step, the approximated pressure is used to update the divergence free velocity at the n+1 step: 𝑢𝑖 𝑛+1−𝑢𝑖 ∗∗ ∆𝑡 = − 𝜕𝑝𝑛+1 𝜕𝑥𝑖 (14) After solving for the velocity and pressure of the fluid domain at the next time point, the position of the particle needs to be updated. Equs 5 and 6 are modified: (𝜌𝑝-𝜌f)𝑉𝑝 𝑑𝑢𝑝 𝑑𝑡 = −𝜌𝑓 ∑ 𝐹𝑙 𝑛+1/2 ∆𝑉𝑙 𝑁𝑙 𝑙=1 + 𝐹𝑐 𝑛+1/2 (15) 𝐼𝑝 𝑑𝜔𝑐 𝑑𝑡 = −𝜌𝑓 ∑ 𝑟𝑙 𝑛 𝑥 𝐹𝑙 𝑛+1/2 ∆𝑉𝑙 𝑁𝑙 𝑙=1 + 𝑇𝑐 𝑛+1/2 (16) Collisions were modeled using a short ranged repulsive force, based on the work of Glowinski et. al [12, 13]. 𝐹𝑐,𝑃1 𝑛+1/2 { 0 𝑖𝑓 |�⃗�𝑃1,2| > 2𝑅𝑃 + ∆𝑟𝑐 , 𝜅𝑐 ( 2𝑅+∆𝑟𝑐−�⃗⃗⃗�𝑃1,2 ∆𝑟𝑐 ) 2 �⃗⃗⃗�𝑃1,2 |�⃗⃗⃗�𝑃1,2| 𝑖𝑓 2𝑅𝑃 < |�⃗�𝑃1,2| < 2𝑅𝑃 + ∆𝑟𝑐 , 𝜅𝑐 ( 2𝑅+∆𝑟𝑐−�⃗⃗⃗�𝑃1,2 ∆𝑟𝑐 ) 4 �⃗⃗⃗�𝑃1,2 |�⃗⃗⃗�𝑃1,2| 𝑖𝑓 |�⃗�𝑃1,2| < 2𝑅𝑃 (17) where �⃗�𝑃1,2 = �⃗�𝑃1 − �⃗�𝑃2, ∆𝑟𝑐 is the threshold distance for the repulsive force, and 𝜅𝑐 = 𝜌𝑝𝑉𝑝|𝑔|. 3. Results 3.1 Model Verification To verify the accuracy of the computational model, the model was compared to experimental data of increasing complexity. The accuracy of the numerical scheme solving for the Navier- Stokes equation was first verified by comparing the model to experimental data for the lid driven flow in a square cavity for Reynolds number of 100, 400, and 1000, Fig 3a. Cross-sections of fluid velocity predicted by the model were compared to Ghia et. al [14] resulting in an average relative error of 0.802%. The accuracy of the IBM was then evaluated by comparing the model to another classical problem in fluid mechanics, flow past a stationary cylinder modeled with IBM mesh, Fig 3b. For mid-range Re, from 10 to 100, steady vortices downstream of the cylinder were formed. Experimental data for the structure of the steady vortices were compared to model predictions with an average relative error of 4.139%. These steady vortices become unsteady at higher Re which transitioned to unsteady von Karmann vortices for higher Re. The vortex shedding creates oscillations in lift and drag forces on the cylinder, Fig 3c. Finally, the coefficient of drag for the cylinder was compared to experimental data over varying Re. The model predicted the coefficient of drag on a cylinder for Re = 0.1, 1, 10, 30, 120, and 1000 with an average relative error of 3.016%, Fig 3d (red dots). Figure 3. Model Verification. (A) Results from Ghia et. al [14] (black dots) are compared to the model predicted for fluid velocity across mid-sections (grey lines) for the lid driven cavity flow. Normalized velocity magnitude is plotted as a contour plot. (B) Pseudo streamlines for the flow over a cylinder with Re=30 and definitions of structural parameters of the steady vortices. (C) Comparison of experimental data [15] for the flow over a cylinder with predicted values from the model. (D) Coefficient of lift and drag for the flow over a cylinder at a Re of 100. (E) Comparison of experimental data to model predicted values for the coefficient of drag for varying Re. 3.2 Inertial Focusing With both the accuracy of numerical methods for solving the acceleration of the fluid and solid phase verified, the model was used to predict the inertial focusing of particles. Experimental data was collected for the migration of 7.32 µm particles, in a straight 50 µm wide PDMS channel (blockage ratio of 0.146), into two beams of particles off the side of the walls. Images were collected every 500 µm downstream of the inlet with short and long exposure to monitor the distribution of particles along the width of the channel (W), as well as their velocities, Fig 4a,b. Cross sections of the channel were analyzed to quantify the radial distribution of particles across the width of the channel, Fig 4c. Model predictions for the distribution of particles along the width of the channel were compared to the experimental data. The model predicted the final equilibrium position of the particles at 0.120 y/W from the wall, Fig 4d. Experimental data found the particles migrating to 0.1255 y/W, which is in agreement with previously measured equilibrium position of 0.125 y/W [16]. The relative error in the model prediction of 4%, or 0.025 µm which is well within in the standard deviation of the experimental data, 2.5 µm. Fig 4. Particle Migration Across Streamlines. (A) Light microscope image of 50 µm channel with 7.32 µm beads. (B) Fluorescents of the 7.32 µm beads at a Re of 20 and an exposure of 1/10 s. (C) Averaged cross section of fluorescence normalized to FWHM. (D) Experimental data (black dots) for the average radial position of particles downstream of the inlet compared to predicted results from IBM simulations (navy line). 3.3 Particle Capture Experimental data for motion of 1, 7.32, and 20 µm beads was collected for varying Re to quantify the size selective capture of particles as a function of well geometry and Re. The 1 µm fluorescent beads were used to evaluate the Re dependent development of the microvortices within the wells (not included). Fig 5. Particle Capture in Microvortices. (A) Long exposure images of 7.32 µm beads for Re = 1, and (B) Re = 100. (C) Short exposure images were then used to collect PIV Data for Re =1 and (D) Re = 100. (E) Comparison of u- velocity of mid-slice of PIV/PSV data (black dots) and model predictions (blue line) for Re = 1 and (F) Re = 100. Notably the difference between the data point circled in red and the model predictions for the fluid is the inertia of the particle that keeps the particle moving quickly as it enters the well before it decelerates. 4.1 Error in PIV/PSV algorithm The algorithm for quantifying the velocity of particles within the fluid performed well for the motion of particles in small and mid-range Re flow. The algorithm accurately quantified the motion of particles in the entire system for Re = 1, Fig 5. However, for Re = 100 and 300 the algorithm could only accurately quantify the velocity of particles within and entering the well, Fig 5. Inside the channel (Re = 100 and 300) the particle velocity was to great for the particle streaks to be captured over the background noise, Fig 5. 4.2 Microfluidics device fabrication The current method of fabrication for the expanded chip design relies on a laser cutter to etch the design. This not only produces rough boundaries that will likely affect the inertial focusing of the particles, but also limits the channel size to >150µm which is too large to observe inertial focusing of 7.32µm which limits the chip design to either capturing 7.32µm beads with a larger well or only capturing 20µm beads. To increase the resolution of the well microfluidics chip, the design will be fabricated using standard photolithography in the same way the straight chip design was fabricated. 4.3 Microfluidics device optimization With the model able to predict the inertial focusing of the particles, as well as the development of microvortices at varying Re, the model will be used to aid in the design and optimization of the well microfluidics devices. Utilizing IBM any type of well geometry can be modeled as well as allow for the modeling of elastic, deformable cells. These simulations can be used to evaluate the rate of particle drift, which is directly related to their size and deformability. The opening of the well can then be modified to increase or decrease the amount of time the particles spend moving past the opening, thus allowing or preventing a size specific particle from drifting into the well and becoming captured. The shape of the well can also be modified to modulate the strength of the vortices to aid in particle retention. 4.4 Biomedical engineering applications Size selective capturing has numerous applications for filtering biological samples. The filtration of cancerous cells out of the blood stream for analysis is a practical use of these chip designs as the larger cancerous cells will experience a larger shear force and become captured in the wells much more readily than red blood cells [2]. Inertial focusing is also greatly influenced by deformability of the particle; increased deformability shifts the focusing of the particle towards the centerline. This could be taken advantage of by filtering out malaria infected red blood cells, which become rigid during the incubation of the viral infection. 5. Conclusion Utilizing continuous and direct forcing immersed boundary methods we developed a model that accurately simulates the inertial focusing of particles. The model was then expanded to model the expansion of the channel with IBM and simulate the capture of 7.32 µm particles within the wells. Future work will aim at utilizing the developed IBM model to optimize the geometry to enhance size selectivity and efficiency of the particle capture. The model will also be expanded to account for deformable and rigid cells. References [1] J. Zhou, S. Kasper, and I. Papautsky, "Enhanced size-dependent trapping of particles using microvortices," Microfluid Nanofluidics, vol. 15, no. 5, pp. 611-623, 2013, doi: 10.1007/s10404-013-1176-y. [2] S. C. Hur, A. J. Mach, and D. Di Carlo, "High-throughput size-based rare cell enrichment using microscale vortices," Biomicrofluidics, vol. 5, no. 2, pp. 022206-022206-10, 2011, doi: 10.1063/1.3576780. [3] A. J. Mach and D. Di Carlo, "Continuous scalable blood filtration device using inertial microfluidics," Biotechnol. Bioeng, vol. 107, no. 2, pp. 302-311, 2010, doi: 10.1002/bit.22833. [4] G. SegrÉ and A. Silberberg, "Radial Particle Displacements in Poiseuille Flow of Suspensions," Nature (London), vol. 189, no. 4760, pp. 209-210, 1961, doi: 10.1038/189209a0. [5] W.-P. Breugem, "A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows," Journal of computational physics, vol. 231, no. 13, pp. 4469-4498, 2012, doi: 10.1016/j.jcp.2012.02.026. [6] J. Lee and D. You, "An implicit ghost-cell immersed boundary method for simulations of moving body problems with control of spurious force oscillations," Journal of computational physics, vol. 233, pp. 295-314, 2013, doi: 10.1016/j.jcp.2012.08.044. [7] A. Gronskis and G. Artana, "A simple and efficient direct forcing immersed boundary method combined with a high order compact scheme for simulating flows with moving rigid boundaries," Computers & fluids, vol. 124, pp. 86-104, 2016, doi: 10.1016/j.compfluid.2015.10.016. [8] J. Yang and F. Stern, "A simple and efficient direct forcing immersed boundary framework for fluid–structure interactions," Journal of computational physics, vol. 231, no. 15, pp. 5029-5061, 2012, doi: 10.1016/j.jcp.2012.04.012. [9] A. M. Roma, C. S. Peskin, and M. J. Berger, "An Adaptive Version of the Immersed Boundary Method," Journal of computational physics, vol. 153, no. 2, pp. 509-534, 1999, doi: 10.1006/jcph.1999.6293. [10] M. Uhlmann, "An immersed boundary method with direct forcing for the simulation of particulate flows," Journal of computational physics, vol. 209, no. 2, pp. 448-476, 2005, doi: 10.1016/j.jcp.2005.03.017. [11] R. Lynch, "Fundamental Solutions of 9-point Discrete Laplacians; Derivationand Tables," ed. Purdue e-Pubs, 1992. [12] R. Glowinski, T. W. Pan, T. I. Hesla, D. D. Joseph, and J. Périaux, "A Fictitious Domain Approach to the Direct Numerical Simulation of Incompressible Viscous Flow past Moving Rigid Bodies: Application to Particulate Flow," Journal of computational physics, vol. 169, no. 2, pp. 363-426, 2001, doi: 10.1006/jcph.2000.6542. [13] M. H. Abdol Azis, F. Evrard, and B. van Wachem, "An immersed boundary method for flows with dense particle suspensions," Acta mechanica, vol. 230, no. 2, pp. 485-515, 2019, doi: 10.1007/s00707-018-2296-y. [14] U. Ghia, K. N. Ghia, and C. T. Shin, "High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method," Journal of computational physics, vol. 48, no. 3, pp. 387-411, 1982, doi: 10.1016/0021- 9991(82)90058-4. [15] M. Coutanceau and R. Bouard, "Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 2. Unsteady flow," J. Fluid Mech, vol. 79, no. 2, pp. 257-272, 1977, doi: 10.1017/S0022112077000147. [16] J. Zhou and I. Papautsky, "Fundamentals of inertial focusing in microchannels," Lab Chip, vol. 13, no. 6, pp. 1121-1132, 2013, doi: 10.1039/c2lc41248a.