Shape-design for stabilizing micro-particles in inertial microfluidic flows

TL;DR

This paper introduces a numerical framework for designing micro-particles with shapes optimized for stable central positioning in inertial microfluidic flows, using shape parametrization, stability metrics, and Bayesian optimization.

Contribution

It presents a novel shape-design methodology for micro-particles that maximizes stability at the channel centerline, incorporating NURBS parametrization and adaptive Bayesian optimization.

Findings

Characteristic 'fish', 'bottle', 'dumbbell' shapes enhance stability.

High k and Re particles show increased stabilization.

Modified dumbbell designs outperform previous rod-disk models.

Abstract

Design of microparticles which stabilize at the centerline of a channel flow when part of a dilute suspension is examined numerically for moderate Reynolds numbers ($10 \le Re \le 80$). Stability metrics for particles with arbitrary shapes are formulated based on linear-stability theory. Particle shape is parametrized by a compact, Non-Uniform Rational B-Spline (NURBS)-based representation. Shape-design is posed as an optimization problem and solved using adaptive Bayesian optimization. We focus on designing particles for maximal stability at the channel-centerline robust to perturbations. Our results indicate that centerline-focusing particles are families of characteristic "fish"/"bottle"/"dumbbell"-like shapes, exhibiting fore-aft asymmetry. A parametric exploration is then performed to identify stable particle-designs at different k's (particle chord-to-channel width ratio) and Re's ($0.1 \le k \le 0.4, 10 \le Re \le 80$). Particles at high-k's and Re's are highly stabilized when compared to those at low-k's and Re's. A comparison of the modified dumbbell designs from the current framework also shows better performance to perturbations in Fluid-Structure Interaction (FSI) when compared to the rod-disk model reported previously (Uspal & Doyle 2014) for low-Re Hele-Shaw flow. We identify basins of attraction around the centerline, which span larger release-angle-ranges and lateral locations (tending to the channel width) for narrower channels, which effectively standardizes the notion of global focusing in such configurations using the current stability-paradigm. The present framework is illustrated for 2D cases and is potentially generalizable to stability in 3D flow-fields. The current formulation is also agnostic to Re and particle/channel geometry which indicates substantial potential for integration with imaging flow-cytometry tools and microfluidic biosensing-assays.

Figures

Tables

Table 1. Table 1: Equilibrium values for R ​ e p = 15 𝑅 subscript 𝑒 𝑝 15 Re_{p}=15 ( r e ​ q , u p , u f , ω p , d ​ u d ​ y subscript 𝑟 𝑒 𝑞 subscript 𝑢 𝑝 subscript 𝑢 𝑓 subscript 𝜔 𝑝 𝑑 𝑢 𝑑 𝑦 r_{eq},u_{p},u_{f},\omega_{p},\frac{du}{dy} denote, resp., stable location, particle linear-velocity, undisturbed fluid velocity at centroid-height of particle, angular-velocity of particle, and velocity-gradient at particle-centroid)

	Yang et. al. (2006)	Present	Error (%)
Equilibrium Location, $\frac{2r_{eq}}{H}$	0.454	0.446	1.652
Linear Velocity, $\frac{u_p}{u_f}$	0.963	0.973	1.038
Angular Velocity, $\frac{{\omega }_p}{\left(\frac{1}{2}\frac{du}{dy}\right)}$	0.9289	0.9418	1.389

Table 2. Table 2: Non-dimensional lift-gradients, d ​ L ∗ d ​ y ∗ 𝑑 superscript 𝐿 𝑑 superscript 𝑦 \frac{dL^{*}}{dy^{*}} , evaluated for a circular particle at the stable equilibrium point, at select confinements, k 𝑘 k , and Reynolds numbers, R ​ e 𝑅 𝑒 Re (the negative sign indicates a restoring lift) - used for computing damping coefficients of non-circular shapes by means of a dynamic shape-factor [ 28 ]

	$Re=10$	$Re=20$	$Re=40$	$Re=60$	$Re=80$
$k=0.1$	-0.0246	-0.0237	-0.0201	-0.0159	-0.0103
$k=0.2$	-0.1522	-0.1444	-0.1297	-0.1041	-0.0731
$k=0.3$	-0.2597	-0.2561	-0.2523	-0.2378	-0.2156
$k=0.4$	-0.3068	-0.3831	-0.4325	-0.4490	-0.4478
$k=0.5$	-0.4660	-0.3917	-0.3810	-0.4881	-0.5064

Equations

$$\nabla \cdot \mathbf{u}=0$$

$$\mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p+\frac{1}{Re}{\nabla }^2\mathbf{u}$$

$$\displaystyle F_{x}=2\Big{[}\oint_{\Gamma_{p}}\big{\{}-p\overline{\overline{I}}+\frac{1}{Re}(\nabla{\bm{u}}+\nabla^{T}{\bm{u}})\big{\}}\cdot\hat{\bm{n}}\,d\Gamma\Big{]}\cdot\hat{\bm{i}}=0$$

$$u(\mathbf{x}\in \mathrm{walls})=-u_p$$

$$u(\mathbf{x}\in {\Gamma }_p)=0,v(\mathbf{x}\in {\Gamma }_p)=0$$

$$F_y(y_p=0,{\theta }_p=0)$$

$${\tau }_z(y_p=0,{\theta }_p=0)$$

$$m\frac{{\mathrm{d}}^2{y}_p}{\mathrm{d}{t}^2}=F_y(y_p,{\theta }_p)-\alpha \frac{\mathrm{d}{y}_p}{\mathrm{d}t}(\alpha >0)$$

$$I\frac{{\mathrm{d}}^2{\theta }_p}{\mathrm{d}{t}^2}={\tau }_z(y_p,{\theta }_p),$$

$$\frac{\mathrm{d}{y}_p}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{v}_{particle}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\theta }_p}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\omega }_{particle}}{\mathrm{d}t}$$

$$\displaystyle{A}=\bigg{(}\frac{\mathrm{\partial}\bm{f}}{\mathrm{\partial}\bm{X}}\bigg{)}_{\bm{X_{0}}}$$

$$A={\left[\begin{array}{cccc}\frac{\partial f_1}{\partial y_p}& \frac{\partial f_1}{\partial v_{particle}}& \frac{\partial f_1}{\partial {\theta }_p}& \frac{\partial f_1}{\partial {\omega }_{particle}}\\ \frac{\partial f_2}{\partial y_p}& \frac{\partial f_2}{\partial v_{particle}}& \frac{\partial f_2}{\partial {\theta }_p}& \frac{\partial f_2}{\partial {\omega }_{particle}}\\ \frac{\partial f_3}{\partial y_p}& \frac{\partial f_3}{\partial v_{particle}}& \frac{\partial f_3}{\partial {\theta }_p}& \frac{\partial f_3}{\partial {\omega }_{particle}}\\ \frac{\partial f_4}{\partial y_p}& \frac{\partial f_4}{\partial v_{particle}}& \frac{\partial f_4}{\partial {\theta }_p}& \frac{\partial f_4}{\partial {\omega }_{particle}}\end{array}\right]}_{{\mathbf{X}}_0}$$

$$={\left[\begin{array}{cccc}0& 1& 0& 0\\ \frac{1}{m}\frac{\partial F_y(y_p,{\theta }_p)}{\partial y_p}& -\frac{\alpha }{m}& \frac{1}{m}\frac{\partial F_y(y_p,{\theta }_p)}{\partial {\theta }_p}& 0\\ 0& 0& 0& 1\\ \frac{1}{I}\frac{\partial {\tau }_z(y_p,{\theta }_p)}{\partial y_p}& 0& \frac{1}{I}\frac{\partial {\tau }_z(y_p,{\theta }_p)}{\partial {\theta }_p}& 0\end{array}\right]}_{{\mathbf{X}}_0}$$

$$\mathbf{P}(\xi )=\frac{{\sum }_{i=0}^n{w}_i{\mathbf{P}}_i{N}_{i,m}(\xi )}{{\sum }_{i=0}^n{w}_i{N}_{i,m}(\xi )}$$

$$N_{i,m}(\xi )=\frac{\xi -t_i}{t_{i+m-1}-t_i}{N}_{i,m-1}(\xi )+\frac{t_{i+m}-\xi }{t_{i+m}-t_{i+1}}{N}_{i+1,m-1}(\xi )$$

$$N_{i,1}(\xi )=\{\begin{array}{ll}1,& \xi \in [t_i,t_{i+1})\\ 0,& \xi \notin [t_i,t_{i+1})\end{array}$$

$$\mathbf{T}=\{t_0=t_1=\cdots =t_{m-1}<t_m\le t_{m+1}\le \cdots \le t_n<t_{n+1}=\cdots =t_{n+m}\}$$

$$\underset{Y\equiv Y_1,Y_2,...,Y_{nv}}{\mathrm{argmin}}\max (\mathrm{Real}({\lambda }_i(\mathbf{P})))+\mathrm{Regularization}$$

$$m_c\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{t}^2}=L_c(y)-{\alpha }_c\frac{\mathrm{d}y}{\mathrm{d}t}({\alpha }_c>0)$$

$$\displaystyle{\Big{(}\frac{\alpha_{c}}{m_{c}}\Big{)}}^{2}+\frac{4}{m_{c}}\frac{dL}{dy}<0$$

$$\displaystyle\frac{dL}{dy}<0\implies-2\sqrt{m_{c}\Bigl{|}{\frac{dL}{dy}}\Bigr{|}}<\alpha_{c}<2\sqrt{m_{c}\Bigl{|}{\frac{dL}{dy}}\Bigr{|}}$$

$$\displaystyle\alpha_{c}>0\implies 0<\alpha_{c}<2\sqrt{m_{c}\Bigl{|}{\frac{dL}{dy}}\Bigr{|}}$$

$$\alpha =\frac{K_s{d}_{v,s}}{K_c{d}_{v,c}}{\alpha }_c$$

$$\displaystyle\mathrm{where,}\;\alpha_{c}=\sqrt{m_{c}\Bigl{|}{\frac{dL}{dy}}\Bigr{|}}$$

$$MLE(\theta )=\log ({\mathbf{y}}^T\mathbf{c})+\frac{1}{n}\sum _{i=1}^n\log ({\lambda }_i(K))$$

$$C^{\ast }=C+\alpha {\int }_S{\kappa }^{2}ds$$

$$C^{*}=C+\log_{10}\bigg{(}\int_{S}\kappa^{2}ds\bigg{)}$$

Full text

1 Introduction

A particle released in flow undergoes a time-evolution of position and velocity as governed by the net forces and torques acting on it, and its long-term behavior is a strong function of its shape. It is widely known that for Stokes number (Stk, the ratio of the characteristic time of the particle to that of the flow) $\ll$ 1, particles flow with negligible inertia along fluid streamlines [1]. However, for larger values of Stk, particle inertia becomes prominent, introducing a variety of nonlinear behavior in flow which can critically affect multiphase-flow applications across scales, e.g., massive oil drills with slurry transport, or lab-on-a-chip devices dealing with hemodynamics, red-blood cell Fluid-Structure Interaction (FSI) [2], and micro/nano-biosensors [3]. In the context of fluid inertia, early studies of the passive manipulation of non-inertial particles date to the mid-twentieth century when randomly-dispersed, rigid, spherical particles released in laminar pipe flow were observed to concentrate to annuli at $\approx 0.6R$, where $R$ is the pipe radius [4, 5]. Spheroids, bodies-of-revolution, ellipsoids at low-$Re$ in uniform, linear-shear, or unbounded paraboloidal flow, resp., were examined analytically by [6, 7, 8]. Subsequent studies extended these results for bounded flows with finite fluid inertia numerically for focusing of arbitrary shapes [9, 10, 11]. Experiments by [12] using disks, cylinders, and h-shapes further generalized the focusing characteristics of arbitrary particles. More recent studies have advanced insight into the low-$Re$ dynamics of irregular particles. [13] reported on the dynamics of torus-shaped particles modelled to represent polymers. They obtained a closed-form solution to the translation velocity of a rotating, force-free torus particle as a function of its slenderness ratio and angular velocity, and numerically studied its translation along a cylindrical rail. [14] studied the motion of thin axisymmetric particles placed in a low-$Re$ linear shear flow, where they calculate power-laws for effective aspect-ratios for families of particle-shapes to control tumbling orbits, and explore shape-tunability. [15] experimentally studied trajectories of micron-scale glass rods in a microfluidic channel. They examine the deviation from traditional Jeffery orbits based on the degree to which axisymmetry is broken. Steady, and non-tumbling motion of particles is crucial to the performance of high-throughput optical scanning devices such as sheath-flow cytometers. The success for capture of key physical, biochemical, and morphological characteristics of the investigated cells is a direct consequence of their spatial orientation at the detection-region [16]. These studies represent a “forward problem” of the fully-coupled fluid-particle system in that they detail the motion characteristics of predefined shaped particles in flow, for a given channel/flow-rate. However, as our capabilities to manufacture complex 3D shaped objects has advanced  [17], it is also of interest to study the “inverse problem”, which would then become an interesting engineering question of identifying particle geometries that satisfy certain desirable motion characteristics long after release.

[18] identified ring-shaped particles which do not tumble in shear flow, in contrast to the tumbling behavior of axisymmetric particles reported by numerous earlier works; the particle shapes were derived as perturbations to a circular shape to obtain zero torques on the particle. Recent advances [19] have established self-aligning and centreline-focusing characteristics of asymmetric particles in Hele-Shaw flow - specifically, “dumbbell” and “trumbbell” shapes. Furthermore, modern fabrication techniques [20, 21, 22] allow for scalable fabrication of arbitrary-shaped microparticles, presenting an abundant landscape for the design of customized application-specific microparticles. Understanding and leveraging the behavior of arbitrarily-shaped rigid particles in flow is an increasingly relevant area of study, but the current state-of-the-art is confined to the pursuit of two different thrusts: the forward problem, which seeks to understand the time-evolution of trajectories and potential focusing locations; and the design of particles under assumptions of zero inertia and unbounded-flows for desirable characteristics such as rotational stability, and centerline-focusing. While much of the previously described work is useful in limited context, there are no generalized methodologies that can be utilized for particle design in more complex scenarios, such as flows with finite fluid inertia. Here, nonlinear behavior in the fluid-structure interaction gives rise to competing forces such as shear-gradient and wall-lifts, making design difficult for different flow fields and channel geometries. This motivates the work undertaken in this paper, where we formulate a framework for shape design over arbitrary flow speeds, channel and particle geometries. We introduce an approach for design of stabilizing particles in confined flows at the channel-centerline. The two main applications which form the basis for the current work are sheath-flow cytometry [23], and raft-particles which act as stable microcarriers for biological cell-specimens [24] for imaging flow cytometry. Sheath-flow techniques use streams of co-flows to focus sample-fluid to a narrow region around the channel-centerline for optical interrogation. However, it is well-known that the centerline and its neighborhood is an unstable region [25, 16] for spherical- and disk-shaped particles; such particles also tumble at the off-centre stable points. We propose to improve these systems by designing dynamically-stable particles that focus to the centerline for such applications which guide particles to neighborhoods around the centerline using co-flows upstream of the scanning region, thereby transforming a purely-active particle manipulation technique to a hybrid active-passive manipulation technique. The inherent particle stability could reduce the sheath-flow volumes required to constrict conventional particles to the channel centerline. We illustrate a methodology for particle-designs in 2D channel-flows, for centerline-focusing to perturbations about the centerline. The proposed methods are easily generalizable to 3D flow-fields.

The rest of the paper is organized as follows: section 2 details the numerical methods employed for the optimization problem, section 3 details the design for centerline-stability with corresponding validation using FSI, and parametric design across a range of $k$’s and $Re$’s with relative comparison of performance of designed-particles using full fluid-structure interaction simulations for centerline-stability. We conclude in section 4 by framing several open questions and subsequent avenues of work. The appendices contain mathematical details of several steps in the framework: A provides validation-cases for the quasi-dynamic approach; B details the damping-coefficient calculation for an arbitrary particle; C describes test-cases for the fixed-budget optimization problem for convergence in total number of required iterations; and D details the convergence studies of the shape-parametrization employed.

2 Methods

2.1 Evaluating stability of an arbitrary particle

We seek a particle design $\bm{P}$ that is stable at the channel centerline ($y_{p}=0$), with its longitudinal axis aligned with the flow direction, $\theta_{p}=0$ (see FIG. 1 for a schematic). The particle should satisfy the following two requirements:

1. (A)

Equilibrium: The nominal configuration $(y_{p}=0,\theta_{p}=0)$ should be a force free (and torque free) configuration. Furthermore, we expect no lateral particle velocity ($v_{p}=0$) as well as zero angular velocity ($\omega_{p}=0$) at this configuration. This will ensure that the particle will remain in this orientation. 2. (B)

Stable equilibrium: The nominal configuration $(y_{p}=0,\theta_{p}=0)$ should be stable to $\delta y$, $\delta\theta$ perturbations. This will ensure that the particle restores back to its nominal configuration after any perturbation.

We evaluate requirement (A) using a quasi-dynamic approach [26, 27], where the steady-state Navier-Stokes equations are solved in a translating reference frame moving at $u_{p}$, the particle streamwise velocity. The approach assumes that the particle instantaneously achieves a streamwise velocity, $u_{p}$, that results in the net axial force, $F_{x}$ becoming zero. This assumption has been previously validated for spherical particles in inertial flow [26, 27]. In this scheme, a particle with no lateral and angular velocity will also have no linear and angular velocity in the moving reference frame. The Navier-Stokes equations around the particle is solved to compute the fluid flow field $\bm{u}$, pressure field $p$, and the unknown streamwise velocity $u_{p}$:

[TABLE]

where $\Gamma_{p}$ is the surface of the particle, $\bm{P}$. The inlet and outlet boundary conditions are chosen to have fully developed parabolic velocity profiles. The particle is far enough away from inlet and outlet for the local disturbance around the particle to not affect the boundary velocity profiles - this lets us use parabolic velocity profiles at the boundaries equivalent to a particle flowing in a long channel (where the current domain is a “section” of that channel such that fully-developed profiles can be imposed with reasonable confidence; typically the channel length is $\geq 30$ times the characteristic length of the particle). Once the steady state Navier-Stokes equations are solved, we compute the lateral force and torque acting on the particle:

[TABLE]

For the first requirement (equilibrium) to be satisfied, the lift $F_{y}$ and torque $\tau_{z}$ must vanish. This is a quantitative measure to ensure that the centerline position ($y_{p}=0,\theta_{p}=0$) is an equilibrium location for a given particle shape. This requirement is trivially satisfied for (top-down) symmetric particles located at the centerline. We use an in-house finite element method framework to solve the Navier-Stokes equations with prescribed boundary conditions. A schematic of the approach is shown in Fig. 2.

After assessment of requirement (A), we next quantify stability in terms of restoring forces and torques when the particle is perturbed from its stable location, $(y_{p}=0,\theta_{p}=0)$. We again utilize the assumption of the quasi-dynamic approach, i.e. the particle instantaneously achieves a streamwise velocity, $u_{p}$, that results in the streamwise force $F_{x}$ becoming zero 111A common simplification employed for creeping flows is that of the force-free particle. A scaling analysis on the equation of motion reveals that the forces vanish in the limit of $Re\rightarrow 0$, which implies instantaneous equilibrium of the particle in all directions every point along its trajectory. In the case of inertial flows, however, this assumption is valid only in the stream-wise direction along which the particle exhibits relatively faster responses to the underlying flow-field, especially for localized perturbations around the centerline. Typical channel-lengths reported in literature for inertial migration [26] indicate smaller time-scales for stream-wise motion than lateral motion. This assumption has been often used, and well validated for spherical particles  [26, 27]. . We also assume that when the particle is released at its perturbed location, its non-streamwise linear ($v_{p}$) and angular velocity ($\omega_{p}$) components are zero. This assumption rests on the idea that the particle is perturbed impulsively from the centerline and has yet to respond to forces arising from the the now-asymmetric flow field. Due to this, the particle can be assumed to have zero lateral and angular velocities when perturbed. The classic example of a simple pendulum is motivation for this assumption. Whether we manually draw the bob to a certain angle away from the mean and then release it from rest, or we give it a slight push from its mean, the bob will tend back to the mean position, thus indicating that trends in force (and torque) about the mean position play a more prominent role than the initial conditions. We do not consider any time dependent effects, computing only the restoring force $F_{y}$ and torque $\tau_{z}$ after an impulsive perturbation 222This is a rather strong assumption, but provides a consistent estimate of the instantaneous response after an impulsive perturbation. Alternatively, the time-dependant response can be computed after an impulsive perturbation. This turns out to be extremely compute intensive, requiring full scale FSI simulations.. We solve the set of equations defined in Eqn. 1 for a perturbed particle configuration $(y_{p}=\delta y,\theta_{p}=\delta\theta)$, and compute the resulting lift force $F_{y}(y_{p}=\delta y,\theta_{p}=\delta\theta)$ and torque $\tau_{z}(y_{p}=\delta y,\theta_{p}=\delta\theta)$

A simplistic approach to gauging stability is to consider only the sign of the resulting lift and torque responses to perturbations. As long as the responses ensure a restoring motion towards the nominal configuration – i.e., $F_{y}<0$ if $\delta y>0$, and $\tau_{z}<0$ if $\delta\theta>0$ – the nominal configuration can be considered a stable equilibrium. The magnitude of the restorative response can be a quantitative measure of the stability, and can be used to rank order various particle shapes. However, this approach has two disadvantages: (1) it does not consider the coupled effects of a restoring torque and force on linear/angular velocities or particle displacement (i.e., it is non-intuitive whether $F_{y}$ should be $>0$ or $<0$ for $\delta\theta>0$, and similarly whether $\tau_{z}$ should be $>0$ or $<0$ for $\delta y>0$); and (2) it does not account for over-damped scenarios where the particle could oscillate about the nominal configuration.

We instead define simplified equations of motion for the particle based on $F_{y}$ and $\tau_{z}$, which will be used to analyze stability in response to small perturbations:

[TABLE]

\justify

where $\alpha$ is the damping coefficient for the particle in the $y$-direction, $m$ is the mass of the particle, and $I$ its moment of inertia about the z-axis. The damping-coefficient for a general shape is approximated using restoring-lifts on a circular particle in plane-Poiseuille flow at the same $(k,{Re})$ (see TAB. 2). The damping coefficient is arrived at by first computing an approximate damping-coefficient for a hydrodynamically equivalent circular particle assuming an under-damped motion, and then adjusting this value using a dynamic-shape factor [28]. This calculation is detailed in appendix B333Although it is not known apriori whether a hydrodynamically-focussed circular particle exhibits rapid decay (critical damping), we assume the worst-case scenario of under-damped motion to allow for oscillations, which is the slowest compared to critical/over-damped decay, and thus design for the same. Another important note in this regard is that we assume the damping coefficient is independent of the particle-location in the channel, drawing from the Stokes drag-analogy due to low transverse-speed. This is again evident from the typical microfluidic channel-lengths it takes for spherical beads to focus [26]..

The above system of second-order equations (4) is converted into four first-order equations, in terms of \Big{(}y_{p},\frac{dy_{p}}{dt},\theta_{p},\frac{d\theta_{p}}{dt}\Big{)}:

[TABLE]

This represents a first order dynamical system, which we use established stability theory to rigorously quantify stability [29]. Specifically, given a first order equation, $\frac{\mathrm{d}\bm{X}}{\mathrm{d}t}=\bm{f(X)}$, we can evaluate the stability of the system about a (hyperbolic) equilibrium point $\bm{X_{0}}$ by linearizing the system about $\bm{X_{0}}$. The linearized system about the equilibrium point is given by $\frac{\mathrm{d}\bm{X}}{\mathrm{d}t}={A}\bm{X}$, where ${A}$ is the Jacobian of the system expressed as:

[TABLE]

Stability is quantified in terms of the eigenvalues of ${A}$. If the real parts of all eigenvalues $\lambda_{i}$ of ${A}$ are negative, the equilibrium point is considered stable. Increasingly negative eigenvalues indicate faster transit to the equilibrium after perturbations [29].

Comparing equations (6) and (5), with $\bm{X}=[y_{p},v_{particle},\theta_{p},\omega_{particle}]^{T}$ the resulting Jacobian is:

[TABLE]

[TABLE]

\justify

This representation (equations (7)) is a generalization to higher dimensions of the 1D-case of a circular particle focusing within a straight 2D channel, where stability is interpreted in terms of the slope of the lift-versus-transverse coordinate curve at equilibrium locations [30]. For fully-3D flow fields, the Jacobian, $A$, would be of size $10\times 10$, in contrast to the 2D case where $A$ is a $4\times 4$ matrix. The additional 6 rows (and columns) appear due to the remaining three directions; one linear ($z$), and two rotational (about the $x$, and $y$ axis, based on Fig. 1). However, for highly-confined geometries depth-wise (along Z-direction) in 3D, a quasi-2D approximation would permit us to formulate stability of the particle using a $4\times 4$ system (equations (7)), while still maintaining fully-3D flow-fields. We construct the Jacobian using finite-difference based gradients, and assess particle stability using the real-parts of its eigenvalues, $\lambda_{i}$ ($1\leq i\leq 4$):

• •

if $\mathrm{Real}(\lambda_{i})<0\hskip 1.99997pt(i=1,2,3,4)$ - the given equilibrium point is stable, otherwise

• •

the given equilibrium point is unstable

Thus, we quantitatively evaluate the stability of an arbitrarily-shaped particle, without the aforementioned pitfalls regarding dynamic force couplings and overdamping behavior. We next turn to a compact parametrization of particle shape.

2.2 Shape-parametrization

The shape of the particle is represented using Non-Uniform Rational B-Spline (NURBS) curves with B-spline basis functions [31]. This compact representation enables generation of a variety of smooth shapes, with local control on curve-shape using control-point weights. Any point on the curve, $\bm{P}=\bm{P}(\xi)$, is given as a function of the parameter, $\xi\hskip 1.99997pt(0\leq\xi\leq 1)$, as in [32]:

[TABLE]

\justify

where, $\bm{P}_{i}=(X_{i},Y_{i})$ are the control-points, $w_{i}$ are the control-points weights, and $N_{i,m}(\xi)$ are the piecewise polynomial B-spline basis functions of order, $m\,(=4)$, given by:

[TABLE]

[TABLE]

[TABLE]

\justify

where $\bm{T}$ defines the knot vector and knot-spans that govern the continuity of the curve and its derivatives. We use a uniform knot-vector, equal weights ($w_{i}=1$) and constrain a given number of control points at predefined $X$-locations so that they are free to move only along the Y-direction for the design problem. The NURBS-curve defines the top-half of the particle, which is mirrored about the $XZ$ plane to complete the shape. For any given shape, the $X_{i}$’s denote the interior, equally-spaced $X$-coordinates of the variable control-points for the shape which have the same values for all shapes, and $Y_{i}$’s denote the interior $Y$-coordinates. The two ends are fixed on the $X$-axis, so that $X_{nv+1}-X_{0}=1$ (non-dimensional), $Y_{0}=\frac{H}{2}$, and $Y_{nv+1}=\frac{H}{2}$. So, the shape is parametrized as $\bm{P}(\xi;\{Y_{0},Y_{1}...Y_{nv+2}\})$, for a given set of $Y_{i}$’s such that $0.1+\frac{H}{2}\leq Y_{i}\leq 0.5+\frac{H}{2}(1\leq i\leq nv)$. Additionally, before any perturbations, we place the shape such that its centroid coincides with the center of the channel, in both X and Y directions.

2.3 Design problem

We are interested in designing particles that exhibit stability to perturbations. We have quantified stability in terms of the eigenvalues of the Jacobian of the associated dynamical system. We frame the design problem as an optimization problem, i.e, find parameters $\{Y_{1},Y_{2},...,Y_{nv}\}$ that define a particle, $\bm{P}$, that minimizes a cost functional. We choose the cost functional to be the maximum (real) eigenvalue of the system. The eigenvalues represent decay-rates of the solution trajectories of the system. For the current system, each shape would have 4 such eigenvalues. Thus for finding stable shapes, we require that all real-parts of eigenvalues be negative. We identify the least-negative eigenvalue, and desire this eigenvalue to be as negative as possible. This automatically minimizes all other eigenvalues to improve overall restoration-rates. Thus formally the optimization is:

[TABLE]

\justify

Regularization terms are added to ensure that wiggly shapes with large curvature changes are penalized (the two types of regularization terms used are detailed in appendix D). The computational effort in solving the Navier-Stokes equations and the subsequent eigenvalue problem leads to a costly objective function. Therefore, we use a Bayesian strategy to optimize particle shape (the Bayesian frame work is described in detail in appendix C). Convergence with respect to the number of control-points for shape-representation is analyzed first, because it is of interest to weigh significant improvements (if any) in the stability of the particle against the complexity in resolving the actual fabrication process. Specifically fabrication processess like 3D printing [33], stop-flow lithography [34], optofluidic fabrication [35], continuous-flow lithography [21] may not be able to capture fine features in the particle shape (wiggles/nooks) introduced by a larger number of control-points, since effects such as the diffusion of a crosslinking photoinitiator might act to smoothen the shape. We discuss the regularization and convergence aspects in appendix D. Furthermore, the number of original cost-function evaluations required for a reasonable approximation of the response surface is typically far lesser compared to that required with traditional optimization techniques, such as evolutionary algorithms. In the present context, as opposed to the usual approach of optimizing one infill-criterion per update of the Radial-Basis surrogate, we perform asynchronous optimization using multiple infill-criteria by giving a range of weights for exploration (high-variance) vs. exploitation (high-mean), which updates the response surface at multiple points after each iteration. This enables efficient utilization of High-Performance Computing (HPC) resources, a useful method for accelerating convergence in computationally difficulty problems (the reader is directed to [36] for additional details).

3 Results and Discussion

3.1 Centerline design: $k=0.3$, $Re=20$

Particles are designed for local stability around the centerline at zero orientation. This means that once particles have been guided close to the centerline upstream of the test-section using existing techniques such as pinched-flow fractionation [37], they will locally focus and remain rotationally stable at the centerline. Families of fish-like shapes were gathered from the optimization run for the case of $k=0.3$, $Re=20$ using 8 control-points ($nv=6$) as shown in FIG. 3. While the best shapes (FIG. 3a) for this configuration singularly appear to be variants of fishes, the entire range of stable shapes (FIG. 3b) appears to be more varied in terms of the local curvature, area, aspect-ratio, etc., with the fore-aft asymmetry for some designs not being as prominent as with the set of best shapes. This allows some leeway in the fabrication process while retaining essential self-stabilizing characteristics. The pressure fields around stable particles from QD-snapshots at perturbed locations (FIG. 3c, 3d) indicate prominent transverse gradients across the length of the particle, especially at the extremities (sections A-A’, B-B’). Moreover, the asymmetry in these gradients in the fore- compared to the aft-segments acts to stabilize such particles so that a $y_{p}$-perturbation leads to a negative-lift, and positive torque, but a $\theta_{p}$-perturbation gives rise to a positive lift, but negative torque444It should be noted that an examination consisting of lift-vs.-$`y^{\prime}$ or torque-vs.-$`\theta^{\prime}$ trends alone would be misleading due to the inherent coupling of lift and torque as functions of $`y^{\prime}$ and $`\theta^{\prime}$.. It is also interesting to note the presence of features such as an intermediate-“lobe” at the mid-section of stable particles.

3.2 Global stability

The stability metric discussed previously is constructed using perturbations local to the channel centerline. However, it is also of interest to study a global notion of the particle’s tendency to focus to any stable locations in the $y_{p}-\theta_{p}$ space, much like the cross-sectional force-maps used to study particle focusing in inertial migration [38]. For a given ($k$, $Re$) configuration, a particle is placed at different lateral and angular locations throughout the channel, and at each location, that streamwise velocity is solved which yields zero net force (see FIG. 2). For the current configuration ($k=0.3$,$Re=20$), we pick three shapes: highly-stable, and weakly-stable, and unstable, and construct force-torque maps and corresponding $\omega$-limit sets [29] as shown in FIG. 4. For any two stable particles, one is called more stable if the largest real-part of its eigenvalues is more negative than that of the other. For the strongly- and weakly-stable shapes, we see that there are finite basins of attraction at $(y_{p},\theta_{p})=(0,0)$, which is at the channel centerline with zero inclination; for the unstable particle, however, there is no such basin. The basins of attraction for the centerline for both the stable shapes are approximately spanned by $0\leq\frac{y_{p}}{a}\leq 0.5$, and $-0.3\leq\frac{\theta_{p}}{\pi}\leq 0.15$. Practically, these basins serve as a guiding estimate of the feasible release-locations of particles that lead to their focusing to the centerline. Although the current work is formulated for maximal stability of particles to perturbations after they have been focused to the centerline, the basins of attraction provide design-bounds for release-locations before focusing.

We demonstrate the behavior of the three types of shapes in transient, full-coupled FSI simulations for a release-location of $(y_{p},\theta_{p})=(0,-10^{\circ})$ (simulation details in appendix E). The trajectories for these simulations are shown in FIG. 5 (visuals in FIG. 6). We see that there is a conclusive demarcation between the time-evolution of lateral and angular positions for the stable and unstable particles. Specifically, the unstable particle tends to rapidly destabilize and tumble over time, drifting away from the centerline monotonically. The stable particles on the other hand, display strongly contained and restorative trajectories, with a displacement from the centerline that is only a fraction of the height reached by the unstable particle. More importantly, it is seen that the angular displacements asymptotically reach $0^{\circ}$ without oscillations, confirming validity of the damping-coefficient estimation in relation to the rank-ordering of shapes in terms of their individual stabilities. Additionally, from the lateral trajectories for the stable particles, we notice that the stability metrics are reflected well in the full-physics simulations. The highly-stable particle has the smallest initial overshoot from to the angular perturbation, whereas the weakly-stable particle overshoots to twice as much height. However, these overshoots tend to gradually decay over time as shown in FIG. 5c. In FIG. 6, the unstable particle was inspected with the QD-method first to check that it is unstable for both $0^{\circ}$ and $180^{\circ}$ orientations, to ensure its suitability for validating the designed stable shapes against. We picked this particular unstable shape among many others to illustrate the fact that although the stable and unstable shapes in this case both look visually similar (as “fish” shapes), stability cannot be guaranteed on qualitative arguments alone. The erratic nature of the unstable particle in both headfirst and tailfirst orientations is well-captured in the FSI-trajectories, where the particle shows no tendency to stabilize in either orientation, and the rate of tumbling builds over time. Lastly, although we provide an initial angular perturbation, the difference between the stabilizing behavior of the highly-stable vs. weakly-stable particles is reflected in the lateral positions over time, but not the angular alignment. This is illustrative of the fact that the motion in the $y_{p}$ and $\theta_{p}$ directions is fundamentally coupled and cannot be decoupled by examining the gradients in lift or torque alone. The eigenvalues of the stability matrix essentially achieve this by giving information about the stability of the particle, and this is consequently reflected in the trajectories which can be expressed in terms of $e^{-\lambda_{i}t}$.

3.3 Centerline design: $k=0.1-0.4$, $Re=10-80$

We next design families of stable shapes (FIG. 7, 8) for a parametric range of the flow parameters, $0.1\leq k\leq 0.4$, and $10\leq Re\leq 80$. These families indicate that the most-stable particles are classes of “fishes”/“bottles”/“dumbbells”. It is observed that low-confinements, and low-$Re$ tend to give high aspect ratio, rod-like shapes, whereas high-confinements, and high-$Re$ tend to favor low-to-moderate aspect ratio shapes. The asymmetric make-up of the rod-like particles would seem to be a consequence of smaller velocity-gradients across the particle, which the lobe and longer “lever-arm” of the particle would leverage to realign the particle after being perturbed. The trends indicate that in general, there is a good amount of variability between all possible stable shapes for any given ($k$, $Re$) configuration (see FIG. 8), suggesting multiple local minima in the global landscape of the chosen cost-function. However, among optimally stable shapes, a large variance is only seen in families designed for low confinements (see FIG. 7). Additionally, it is seen that shapes which are optimal in the orientation reported here ($\theta_{p}=0^{\circ}$) are highly unstable in the $\theta_{p}=180^{\circ}$ orientation, which suggests uniqueness of the preferred stable orientation in flow.

It is interesting to note that stability of particles increases with an increase in either confinement or $Re$. In the case of spherical particles of a fixed size, the wall-lift force $W_{L}\propto k^{6}Re^{2}$, whereas the shear-gradient lift $W_{SG}\propto k^{3}Re^{2}$ (neglecting slip-shear, and rotation-induced lift). However, for the particles designed in this work, we observe that stability increases with $k$ at constant $Re$ and vice-versa, with increase in stability much larger for an increase in $k$, than that in $Re$. If we assume that the forces scale in a qualitatively similar manner as spherical particles, we can conclude that stability is strongly governed by wall-lift forces. Additionally, we see that there exist shapes which are stable across multiple ($k$, $Re$) configurations. In this context, we note that a number of flow configurations contain what may be termed “modified-dumbbells” as optimal designs for local focusing, which we regard as improvements to the 3D asymmetric dumbbell shapes (rod-disk model) proposed by [19] in Hele-Shaw flow, although the current work is based in 2D. In order to test the performance of the designed shape with the ones reported earlier by [19], we choose the configuration of $k=0.3$, and $Re=60$. The conventional dumbbell shape was reconstructed using values reported by [39] ($\tilde{s}=3.3$ in their work). The QD cost-function on this particle revealed that it is unstable in both orientations, $0^{\circ}$, and $180^{\circ}$, in 2D as well as 3D. When the modified and conventional dumbbell particles were simulated using full-FSI (FIG. 10), it was found that the proposed designs performed better in comparison. Specifically, the lateral and angular trajectories reveal that the initial phase seems qualitatively similar, where both particles tend to restore to the centerline after the initial overshoot from the centerline. Over time, however, the unstable particle rapidly destabilizes away from the centerline, in contrast to the stable design. This suggests that there is much scope yet for improvement to centerline-focusing shapes (either local criteria for near-centerline release as in sheath-flow cytometry, or for global criteria for arbitrary release-locations in the channel as a passive manipulation technique) that have been previously reported for assumptions such as low-$Re$, unbounded flows, and so on.

3.4 Stability under high confinement

We designed particles for $k=0.5$ and $Re=10$ to examine basins of attraction (FIG. 9). From the basins for the wider channel cases (i.e., smaller $k$, see FIG. 4), we see that the basin for $(y_{p},\theta_{p})=(0,0)$ for the present configuration includes nearly the entire lateral range of release-locations, and $-0.35\leq\frac{\theta_{p}}{\pi}\leq 0.1$ for release-angles. This is a significant coverage compared to the wider channels, where the lateral range was restricted to a notably smaller span. This suggests applicability of the design optimization approach in this work toward global focusing, with a significant range of particle release-angles and locations in narrower channels. Such basins of attraction would be even more significant for confinement ratios ($\approx 0.72$) akin to those used in recent works [24]. From the parametric sweep for centerline-stable particles (FIG. 7), we also test the performance of flowing particles for additional configurations (using FSI), varying in $k$, or $Re$, or both. Specifically, we test for response to angular perturbations - $(k,Re)=(0.4,20),(0.4,80)$ - as well as transverse perturbations - $(k,Re)=(0.7,20)$, and trajectories are shown in FIGS. 11, 12, 13, respectively, along with candidate shapes. It is apparent that particles restore monotonically in the direction of the initial perturbation, while the displacement in the other direction is non-monotonic, tending to shift the particle away from equilibrium at first, until stabilizing stresses begin to re-position it to the mean location. As with earlier cases, particles deemed to be unstable by our cost-function drift away significantly from the centerline, accompanied by a tumbling motion.

4 Conclusions

We have demonstrated a computational framework for designing self-stabilizing particles in 2D inertial laminar flow, geared towards microfluidic cell-scanning devices such as sheath-flow cytometers. Stable particle designs group into families of “fish”/“bottles”/

“dumbbells” shapes depending on channel confinement and flow conditions, suggesting the existence of multiple optimal designs per configuration. Designed particles have been conclusively shown to exhibit stability to perturbations in contrast to particles deemed unstable, as verified computationally using two-way coupled FSI simulations. The basins of attraction for wide channels (low $k$) reveal a finite region of release-locations around the centerline for local focusing (particle release near the channel centerline), whereas those for narrow channels tend to cover all lateral locations in the channel, suggesting a higher possibility of global focusing (far-release) to the centerline. The design methodology in the present work has been demonstrated for purely 2D scenarios, but is easily extensible to 3D channels which introduce blunted flow-profiles. The methods discussed herein are envisioned to lay a basis for future work including design for global stability and robustness in system parameters, which will eventually result in optimized designs of particles for high-throughput performance – conceivably with non-Newtonian, complex bio-fluids. In addition, we also see scope for exploring modified shape-descriptors, including non-monotonic shapes, Bezier-PARSEC parameterization, and Elliptic Fourier Descriptors, which may yet reveal a richer phase-space of stable designs. Finally, more in-depth sensitivity analysis on the control-points and the use of low-dimensional models could aid in computational efficiency, and make the framework more readily applied to unconventional channel geometries.

Appendix A Validation

The quasi-dynamic (QD) framework is validated by simulating the Segre-Silberberg effect [4] in 2D [30]. 15 particle-locations are chosen in the radial direction, and particle Reynolds numbers of 1.67, 15, and 37.5 are used, where $k=0.15$. The results are shown in FIG. 1, and are in excellent agreement with previous reports, including particle equilibrium location, and velocities (TAB. 1).

Appendix B Damping-coefficient calculation

If $\alpha_{c}$ is the damping-coefficient, $L_{c}$ is the restoring lift at the stable location, and $m_{c}$ is the mass of the circular cylinder, then the equation of motion in the $y$-direction is:

[TABLE]

\justify

Assuming a linear variation of the position-dependent lift close to the stable-point, $L_{c}(y)=\frac{dL_{c}}{dy}y$, and drawing from the traditional spring-mass-damper system, the condition for under-damped motion for equation (12) becomes,

[TABLE]

\justify

Since we would like to use the damping coefficient as an estimate in the order-of-magnitude sense, any value in the above range should appropriately capture variations in damping with varying particle-shapes and sizes, true to a given configuration and flow parameters. The damping coefficient of any arbitrary shape, $\alpha$, is then computed as:

[TABLE]

\justify

The ${}^{\prime}s^{\prime}$ subscripts in equation (13) refer to the non-circular particle, and the ${}^{\prime}c^{\prime}$ refers to quantities pertinent to the circular particle. $K$, and $d_{v}$ stand for the dynamic shape-factor, and volume-equivalent diameter, as defined in [28]. The $\frac{dL}{dy}$ term is computed using the non-dimensional counterparts from TAB. 2.

Appendix C Fixed-budget problem

The optimization is performed using an in-house radial-basis function surrogate framework based on the ‘matern5/2’ kernel covariance function, along with a lower-confidence bound infill strategy for adaptive sampling [40], [41]. A maximum likelihood estimator is minimized [42] for optimizing the length-scale parameter in the kernel function after each infill as follows:

[TABLE]

\justify

where, $\theta$, is the characteristic length-scale parameter for the basis functions, $\bm{y}$ is the function-vector after the latest infill, $\bm{c}$ is the weight-vector, ${}^{\prime}n^{\prime}$ is the number of sampled locations, and $\lambda_{i}(K)$ are the eigenvalues of the covariance matrix, $K$, where $K^{-1}\bm{y}=\bm{c}$. The initial dataset is generated using a Latin Hypercube Sampling (LHS), and the total number of infills performed is restricted due to the fixed-budget nature of the problem. Metamodelling of the response-surface allows for expensive cost-functions to be conveniently approximated using strategically sampled points, wherein the approximations and their derivatives are trivial to compute in terms of computational cost, from a minimization point-of-view. We formulate the minimization as a fixed-budget problem as limited by the computational resources available. For each optimization run, we add as many infills (locations where the true function is evaluated each iteration as deemed “optimal” according to the infill-criterion) as it takes to arrive at a converged minimum, as also is informed by previous optimizations performed on the original cost-function using evolutionary algorithms. The first test is performed on a two-variable problem, where every candidate-shape is defined by 8 control-points in total where all but two of the control points are fixed. Specifically, the variables $`y_{4}^{\prime}$, and $`y_{7}^{\prime}$ from FIG. 1 are chosen to vary. FIG. 1c illustrates the original function-surface which is created using a $20\times 20$ mesh-grid of the variables, and FIG. 1d depicts all the 400 shapes used for generating it. The optimization is performed using twice, each with a different realization of the initial-sampling for the surrogate. The metamodel response-surfaces from these two runs are shown in FIGS. 1a, & 1b. It is seen that the reconstructions, which required a total of 60 function evaluations, are remarkably similar to the original function-surface, which required 400 points. The metamodel is able to demarcate the stable and unstable regions of the original function reasonably well. Additionally, the infill-criterion appears well-“calibrated” due to the fact that we see a clustering of infill-locations towards the minima-regions of the function and a sparse set farther away, and this resolution of the optimal regions is conducive to the design problem, where we are interested in capturing optima well. Additionally, we use the 2:1 infill-to-initial sampling, and 10 initial-points-per-dimension, thumb-rules [40] for this test, which seems reasonable from findings. Our second test comprises a higher-dimensional optimization ($nv=6$). Contrary to the two-variable case, we run infill-to-initial ratios of 2, 5, and 10, and examine the range of fitnesses for the stable shapes (FIG. 1e). We find that higher stabilities are only achieved for infill-to-initial ratio of about 10, and this more closely matches our benchmarks from the Genetic Algorithm (GA). Although the absolute difference in the best-fitness values between ratios of 5 and 10 is about 16%, we choose a ratio of 10 for all further runs due to the scope of improvement seen here.

Appendix D Convergence in control-points: $k=0.2$, $Re=20$

\justify

We run design-optimization for the following tests:

• •

with 4, 6, 10, and, 14 variable-parametrization - without curvature-penalty regularization

• •

with 6 (without regularization) and, 10 variable-parametrization (with scaled curvat-

ure-penalty)

• •

with 6, and, 10 variable-parametrization - both with log-curvature penalty

For the first case (FIG. 1), it is seen that the best-stability values increase with the number of control-points due to the fact that the shape-representation now allows for introduction of crucial features which act to stabilize the particle more. However, it is seen that the average trend of the shapes remains similar - fore-aft asymmetry characterized by a major aft-segment with one or more smaller fore-segments. But for ${}^{\prime}nv^{\prime}=10,14$ (FIG. 1c, 1d), the shapes contain a significant number of small-scale features/bends which is undesirable. For the second case (FIG. 2), the cost-function ($C^{*}$) for a candidate-shape is evaluated using the stability ($C$) and a scaled curvature-integral (regularization) over the shape:

[TABLE]

\justify

where, $\kappa$ is the curvature at a segment of length, $ds$, on the shape; $\alpha=0.1$ (configuration-dependent) is chosen to offset the large variations in curvature with moderate changes in shapes, so as to not discard potentially-stable shapes. The penalized-shapes have fitness values which are quite similar to the shapes with the unregularized cost-function, and also similar profiles to a certain extent. This essentially suggests that we can do away with the higher-$`nv^{\prime}$ parametrization with regularization by using $`nv^{\prime}=6$, without the added dimensionality. However, it appears that the regularization used for the second case penalizes candidate-shapes drastically, due to which we also investigate a log-curvature-integral regularization to reduce the effect of the penalty:

[TABLE]

\justify

For the third case, we again see that although both $`nv^{\prime}$s produce different profiles, the stability values are close which confirms that lower $`nv^{\prime}$ without penalty would produce shapes of stability similar to those of a higher $`nv^{\prime}$ with regularization. Thus, $`nv^{\prime}=6$ is used for all further optimization runs.

Appendix E Fluid-Structure Interaction (FSI) problem setup

The transient simulations for designed particles are performed using ANSYS Fluent 18.1. The Dynamic-Mesh feature is utilized to account for changing particle-boundary ($\Gamma_{p}$) location every time-step, using a boundary-fitted mesh. The fluid-flow equations are solved dimensionally using a finite-volume approach, and the particle positions and velocities are computed using the 6-DOF rigid body solver, and updated using a multi-step predictor-corrector scheme. The particle chord-length (${}^{\prime}a^{\prime}$) is always taken to be $0.1\,m$, and the fluid density, and viscosity are taken to be, $1000\,\frac{kg}{m^{3}}$, and $25\,Pa-s$. We use no-slip walls (zero-velocity), fully-developed inlet and zero-pressure outlet boundary conditions. The initial conditions are set to be fully-developed velocity and pressure fields throughout the channel, and the particle is released at rest, at a prescribed perturbation. The mesh close to the particle-surface is refined sufficiently to preserve the particle-shape as it deforms, and the time-step imposed is such that the particle does not move more than half the characteristic length of the surface-elements. The particle is taken to be neutrally-buoyant and the mass and moment properties are externally supplied through a User-Defined Function (UDF), which is invoked by the solver concurrently during runtime. A neighborhood around the deforming zone is locally remeshed every few time-steps to ensure a minimum quality of elements based on skewness, for accurate interpolation of flow-fields from previous time-steps.